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

Эта работа ждет рассмотрения преподавателем Дата последней правки страницы: 07.02.2017 Авторы этой статьи считают, что задание выполнено.

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

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

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

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

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

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

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]-го кластера.

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

Шаг 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] 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

Суммирование всех элементов матрицы с заранее определёнными индексами [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)
clusterization = zeros(N) # sparse version of M

while not stop_criteria(J_history):
    J_curr = 0

    # update M, J
    for i in range(N):
        M[i, clusterization[i]] = 0
        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)
    
    # update C
    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 Свойства алгоритма

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

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

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

Недостатки:

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

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

Для демонстрации работы алгоритма была создана собственная реализация алгоритма с использованием фреймворка Intel® Threading Blocks

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

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

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