Участник:Urandon/Строгий алгоритм С средних (Hard C-Means, HCM)

Материал из Алговики
Перейти к навигации Перейти к поиску
Symbol wait.svgЭта работа прошла предварительную проверку
Дата последней правки страницы: 09.02.2017
Данная работа соответствует формальным критериям.
Проверено ASA.


Автор описания Хомутов Никита

1 Свойства и структура алгоритмов

1.1 Общее описание алгоритма

Строгий алгоритм C средних (Hard C-Means; Jang, Sun and Mizutani, 1997[1]) пытается определить центры кластеров в многомерном признаковом пространстве[2]. Цель заключается в том, что бы сопоставить каждой точке признакового пространства соответствующий кластер.

Назначение алгоритма заключается в кластеризации больших наборов числовых данных представленных признаковыми описаниями объектов. Достоинства алгоритма в классе алгоритмов, решающих данную задачу, -- лёгкость реализации и вычислительная простота. Недостатки же -- явное задание количества кластеров, отсутствие гарантии нахождения оптимального решения.

На вход алгоритму задаются точки конечномерного признакового пространства и количество кластеров. Так же могут задаваться пороговые значения, использующиеся в критерии останова алгоритма: величина целевого функционала, скорость изменения функционала при выполнении основного цикла алгоритма. Целевой функционал алгоритма -- среднее внутрикластерное расстояние.

На первом шаге алгоритма производится инициализация центров кластеров. Это может быть осуществлено путём выбора случайных значений. Затем выполняется основной цикл алгоритма:

  1. Расчёт рядовой матрицы, сопоставляющей каждую точку признакового пространства с ближайшим к ней кластерному центру. Иначе это можно рассматривать как оптимизацию целевого функционала по принадлежности точек к кластерам
  2. Расчёт целевого функционала. Проверка условий критерия останова и выход из цикла в случае его выполнения.
  3. Пересчёт кластерных центров как центров масс точек кластера. Представляет собой оптимизацию целевого функционала по координатам центров кластеров.

Алгоритм гарантирует монотонное невозрастание целевой функционала, следовательно, сходимость к локальному минимуму.

Целевая функция алгоритма полностью аналогична целевой функции K-Means[2].

1.2 Математическое описание алгоритма

На вход алгоритм принимает точки [math]\{u_k\}_{k=1}^{N}[/math] признакового пространства [math]\mathbb{R}^D[/math] и количество кластеров [math]C[/math].

На выходе алгоритм возвращает рядовую матрицу[math]M[/math], содержащую информацию о принадлежности элементов кластерам: [math] m_{ij} = \begin{cases} 1, & u_j \in C_i, \\ 0, & \text{otherwise}, \end{cases} [/math] где [math]C_i[/math] — это множество точек, принажлежащих [math]i[/math]-му кластеру, и а [math]c_i[/math] — центр [math]i[/math]-го кластера.

Алгоритм состоит из пяти шагов[3], включающих инициализацию и цикл оптимизации.

Шаг 1 – инициализация центров кластеров

[math] c_i \sim Uniform \{u_1, \dots, u_N\} [/math] -- инициализация путём случайного выбора центров кластеров из точек входного набора. Начальные центры кластером так же могут быть заданы пользователем явно. [math]test[/math]

Шаг 2 – вычисление рядовой матрицы [math]M[/math]

[math]M = \|m_{ik}\|[/math].

Здесь:

[math] m_{ik} = \begin{cases} 1, & \|u_k-c_i\|^2 \leqslant \|u_k-c_j\|^2 \qquad \forall j \neq i \\ 0, & \text{otherwise}. \end{cases} [/math]

где [math]N[/math] — количество элементов во входном наборе данных.

Представляет собой шаг [math]J \to \min_{M}[/math] при условии, что каждый объект попадает хотя бы в один кластер.

Шаг 3 – расчёт целевого функционала

[math] J = \sum\limits_{i=1}^C J_i = \sum\limits_{i=1}^C \left( \sum\limits_{k:\, u_k \in C_i} \|u_k -c_i \|^2 \right) .[/math]

На этом шаге, после расчёта целевого функционала проверяем критерий останова (зависит от реализации) и принимаем решение, следует ли остановить цикл оптимизации.


Шаг 4 – пересчёт кластерных центров

[math] c_i = \frac{1}{|C_i|} \sum\limits_{k:\, u_k \in C_i} u_k, [/math]

Представляет собой шаг [math]J \to \min_{c_1, \dots, c_N}[/math]


Шаг 5 – зацикливание

Переход по циклу к шагу 2.

1.3 Вычислительное ядро алгоритма

Для каждой пары [math]\{u_i, c_k\}, i \in \mathbb{N}_N, k \in \mathbb{N}_C[/math] согласно алгоритму необходимо вычислить расстояние [math]\|u_i - c_k\|^2[/math] в [math]D[/math]-мерном пространстве. Такие расстояние вычисляются для [math]NC[/math] пар. Таким образом, происходит вычисление [math]O(NCD)[/math] арифметичеких операций.

1.4 Макроструктура алгоритма

Алгоритм представляет собой итерационный процесс, в котором поочерёдно происходит покоординатная оптимизация функционала [math]J[/math].

[math] J(M, C) \to \max_{M} [/math]

[math] J(M, C) \to \max_{C} [/math]

Оптимизация функционала по [math]M[/math] приводит к вычислению расстояний между объектами и центроидами, поиском ближайших.

Оптимизация по [math]C[/math] приводит к пересчёту центроидов как центров масс новых конфигураций кластеров.


Алгоритм распадается на макрооперации, примитивами которых выступают вектора линейных пространств.

Вычисление разностей векторов с последующим вычислением евклидовой нормы на втором шаге алгоритма. [math] d_{ik} = \|u_k - c_i\|^2 [/math]

Нахождение аргминимума в векторе при фиксированном объекте [math]k[/math]: [math] i^* = \arg\!\min\limits_{i = 1 \dots C} d_{ik} [/math]. После чего [math]m_{ik} := [i = i^*][/math]. Данная операция сводится к использованию ассоциативной [math]\min[/math] и паттерна Reduse[4]


Суммирование всех элементов матрицы с заранее определёнными индексами [math] J = \sum\limits_{i=1}^C J_i = \sum\limits_{i=1}^C \left( \sum\limits_{k:\, u_k \in C_i} d_{ki} \right) = \sum\limits_{k=1}^N d_{k, C(u_k)} [/math]

Умножение вектора на скаляр и сложение векторов [math] c_i = \frac{1}{|C_i|} \sum\limits_{k:\, u_k \in C_i} u_k, [/math]

1.5 Схема реализации последовательного алгоритма

Последовательная схема реализована в виде python-подобного псевдокода. Для оптимизации вычислений хранится структура, хранящая номер кластера для каждой точки. Это оптимизация с учётом разреженности матрицы [math]M[/math] 

# Инициализация
J_history = []
# центроиды выбираются случайно среди объектов
c = u[random.choice(N, C), :] 
M = zeros(N, C)
# хранилище индексов кластеров для каждого объекта
# используется для обновления ранговой матрицы M
clusterization = zeros(N) 

# Основной цикл оптимизации
# повторяем его пока не сработает критерий останова
# предполагаем, что критерий использует значения
# функционала качества кластеризации
while not stop_criteria(J_history):
    J_curr = 0

    # Оптимизация матрицы M
    # при фиксированных центроидах
    # пообъектный цикл
    for i in range(N):
        # вычисляем расстояния до кластеров
        # находим минимальный
        d_min, j_min = +inf, None
        for j in range(C):
            d_ij = 0
            for k in range(D):
                d_ij += (u[i,k] - c[j,k])**2
            if d_ij < d_min:
                d_min, j_min = d_ij, j
        # добавляем расстояние до ближайшего 
        # кластера к фунционалу качества кластеризации
        J_curr += d_min

        # обнуление старого соответствия объекта кластеру
        # передача соответствия новому кластеру
        M[i, clusterization[i]] = 0
        clusterization[i] = j_min
        M[i, clusterization[i]] = 1
 
    J_history.append(J_curr)
    

    # Оптимизация центроидов при 
    # фиксированной ранговой матрице
    # покластерный цикл
    c = zeros(N, C)
    for j in range(C):
        # заводим счётчик, сколько объектов в кластере
        # он далее будет использоваться для
        # нормализации центроида
        n_pts_in_cluster = 0

        for i in range(N):
            # складываем центроид покоординатно только
            # с точками, принадлежащими кластеру
            if not M[i, j]:
                continue
            n_pts_in_cluster += 1
            for k in range(D):
                c[j, k] += u[i, k]

        # нормализация центроида
        for k in range(D):
            c[j, k] /= n_pts_in_cluster

1.6 Последовательная сложность алгоритма

На каждой итерации основного цикла алгоритма:

  • Шаг 2 [math]NC[/math] операций вычисления квадратов расстояний до центров кластеров: [math]NCD[/math] разностей, возведеий в квадрат и [math]NC(D-1)[/math] сумм
  • Шаг 2 [math]NC[/math] сравнений
  • Шаг 3 [math]N[/math] векторных сложений: [math]NCD[/math] сложений
  • Шаг 4 [math]N[/math] векторных сложений, [math]C[/math] векторных умножений на скаляр, [math]C[/math] сложений: [math]ND+C[/math] сложений и [math]CD[/math] умножений.

Итого, одна итерация основного цикла стоит [math]O(3NCD)[/math] элементарных операций. Количество итераций зависит от данных, от начального приближения кластеризации, от критерия останова.

1.7 Информационный граф

В информационном графе выделяются различные структурные группы вершин.

Группа 1 Обозначена [math]d(p,c)[/math], выделена зелёным цветом. Отвечает за вычисление расстояний между точками [math]p, u_i[/math] и центроидами [math]c[/math]. Расстояния вычисляются независимо как по точкам, так и по центроидам. Несмотря на независимость вычислений, их следует группировать по точкам, так как в дальнейшем результатов (расстояния до центров) будут аггрерироваться по точкам (операция [math]\arg\min[/math]).

Группа 2 Обозначена [math]\arg\min[/math], выделена жёлтым цветом. Отвечает за поиск ближайшего центроида до данной точки (при уже вычисленных расстояниях). Фактически, применяет паттерн Reduce с участием ассоциативной операцией [math]\arg\min[/math]. Находит как индекс ближайшего центроида (используется для обновления матрицы [math]M[/math]), так и расстояние до него (входит в сумму [math]J[/math]).

Группа 3 Обозначена [math]M[/math], выделена красным цветом. Обновляет значения матрицы [math]M[/math] для каждой точки. Индикатор, соответствующий старому центроиду обнуляется, а новому - выставляется 1.

Группа 4 Обозначена [math]J sum[/math], выделена сиреневым цветом. Осуществляет суммирование расстояний до ближайших центроидов. Применение паттерна Reduce с участием операции [math]+[/math].

Группа 5 Обозначена [math]upd c[/math], выделена синим цветом. Обновляет координаты центроидов. Для каждого центроида [math]c[/math] перебираются все точки [math]p[/math], где [math]M[c, p] = 1[/math]. Координаты таких точек суммируются (независимо покоординатно), возможно, с применением паттерна Reduce по точкам. Так же заводятся счётчики: сколько точек принадлежит кластеру с данным центроидом. Затем покоординатно производится нормировка.


Информационный граф алгоритма. Здесь [math]p[/math] - точки [math]D[/math]-мерного пространства, соответствующие исследуемым объектам (входные данные); [math]d(p,c)[/math] вычисление расстояний от точки то центров; [math]\arg\,\min[/math] вычисление минимума и аргминимума среди расстояний до центров; [math]M[/math] обновление элементов матрицы [math]M[/math]; [math]J, J_{sum}[/math] вычисление и аггрерирование функционала оптимизации; [math]upd c[/math] обновление центров кластеров

1.8 Ресурс параллелизма алгоритма

Так как алгоритм имеет итерационную природу, то оценим ресурс параллелизма одной его итерации.

Из информационного графа четко видна связь последовательная между группами вычислений, поэтому разберём группы по отдельности.

В первой (зелёной) группе необходимо вычислить [math]NC[/math] расстояний. Каждое вычисление расстояния сопоставляется вычислению [math]D[/math] покоординатных расстояний с последующим суммированием. Ширина ЯПФ [math]O(NCD)[/math], высота [math]O(\log D)[/math]. Отметим, что на практике суммрование сдваиванием может не быть самой оптимальной стратегией из-за оверхеда на обеспечение синхронизации: выходом может быть тоже иерархический паттерн Reduce, но арргерирующий больше двух результатов на одной вычислительной единице. Асимптотическая высота такого паттерна так же [math]O(\log NC)[/math]. При достаточно малых значениях [math]D[/math] оптимальная стратения может выродится в цикл на одной вычислительной единице.

Во второй группе (желтой) по каждой из [math]N[/math] точек проводится [math]\min[/math]-аггрегация [math]C[/math] результатов. Таким образом, ширина ЯПФ [math]O(NC)[/math], высота [math]O(\log C)[/math].

В третьей (красной) группе всё примитивно. [math]N[/math] событий по опусканию одного и поднятию другого флага. Ширина [math]O(N)[/math], высота [math]O(1)[/math].

В четвёртой (сиреневой) группе: SUM-аггрерация [math]N[/math] результатов. Ширина ЯПФ [math]O(N)[/math], высота [math]O(\log N)[/math].

В пятой (синей) группе: по каждому из [math]C[/math] центроидов независимо покоординатно ([math]D[/math] координат) надо найти сумму произведений [math]u_{ik} \cdot M_{ij}[/math]. Так же надо завести счётчик. Ширина ЯПФ [math]O(C(D+1)N)[/math], высота [math]O(\log N + 1)[/math].

Итого, ширина ЯПФ алгоритма [math]O(NCD)[/math], высота [math]O(\log N + \log C)[/math]

1.9 Входные и выходные данные алгоритма

Входные данные:

  • [math]\{ u_k \}_{k=1}^N \subset \mathbb{R}^D[/math] — исходный набор точек признакового пространства.
  • [math]C[/math] — количество кластеров.
  • [math]criteria\, params[/math] — параметры критерия останова.

Объём входных данных: [math]ND + const[/math]

Выходные данные:

  • [math]M \in \mathbb{R}^{N \times C} [/math] — рядовая матрица
  • [math]\{c_k\}_{k=1}^C[/math] — координаты кластерных центров

Объём выходных данных: [math]C(N+D)[/math]

1.10 Свойства алгоритма

Назначение: кластеризация больших наборов числовых данных представленных признаковыми описаниями объектов.

Достоинства:

  • лёгкость реализации
  • вычислительная простота

Недостатки:

  • явное задание количества кластеров
  • отсутствие гарантии нахождения оптимального решения
  • чувствительность к начальному приближению координат кластерных центров

Соотношение последовательной [math]O(NCD)[/math] и параллельной [math](\log N + \log C)[/math] сложности в случае неограниченных ресурсов сублинейна по [math]N[/math] и [math]C[/math] (по количеству объектов и кластеров), линейна по [math]D[/math] (по размерности признакового пространства).

Вычислительная мощность, как отношение числа операций к суммарному объёму входных и выходных данных, линейна.

Алгоритм детерменирован при явном задании начального приближения для центроидов. Алгоритм чувствителен к выбору начального приближения.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

Для проверки масштабируемости алгоритма была создана собственная реализация алгоритма[5] с использованием фреймворка Intel® Threading Building Blocks [6].

Хотя реализация имеет ряд ограничений для полноценной оценки масштабируемости: предполагается запуск на одной машине используя мультитрединг -- всё же она может дать интерпретируемые результаты.

Используемый компилятор: оптимизирующий компилятор Microsoft® С/C++ версии 19.00.24215.1. Процессор: Intel Core i7-3157U (2 ядра, 4 потока).

Решалась задача в 128-мерном пространстве с поиском 10 кластеров. Результаты тестов:

№ объектов \ № потоков 1 2 4
100 165ms 89ms 49ms
1000 1612ms 886ms 485ms

Так как ресурсы тестовой системы сильно лимитированы, то получить сложность [math]O(\log N + \log C)[/math] и не удалось бы: распараллелены только покоординатные операции. Видно, что покоординатное распараллеливание даёт реальное ускорение в рамках столь ограниченной системы.

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

1) Matlab, последовательная реализация

2) С++, последовательная реализация

3 Литература

  1. ↑ Jang, J.-S. R., Sun, C.-T. and Mizutani, E. (1997). Neuro-Fuzzy and Soft Computing, number isbn 0-13-261066-3 in Matlab Curriculum Series, Prentice Hall, Upper Saddle River, NJ, USA.
  2. ↑ Jan Jantzen (1998). Neurofuzzy Modelling.
  3. ↑ Нейский И. М. Классификация и сравнение методов кластеризации. // Интеллектуальные технологии и системы. Сборник учебно-методических работ и статей аспирантов и студентов. – М.: НОК «CLAIM», 2006. – Выпуск 8. – С. 130-142.
  4. ↑ Reduce paralle pattern [1]
  5. ↑ https://github.com/urandon/hcm_tbb
  6. ↑ https://www.threadingbuildingblocks.org/
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Urandon/Строгий_алгоритм_С_средних_(Hard_C-Means,_HCM)&oldid=22489