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

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.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) С++, последовательная реализация