Участник:Urandon/Строгий алгоритм С средних (Hard C-Means, HCM)
Содержание
- 1 Свойства и структура алгоритмов
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритмов
1.1 Общее описание алгоритма
Строгий алгоритм C средних (Hard C-Means; Jang, Sun and Mizutani, 1997[1]) пытается определить центры кластеров в многомерном признаковом пространстве[2]. Цель заключается в том, что бы сопоставить каждой точке признакового пространства соответствующий кластер.
Назначение алгоритма заключается в кластеризации больших наборов числовых данных представленных признаковыми описаниями объектов. Достоинства алгоритма в классе алгоритмов, решающих данную задачу, -- лёгкость реализации и вычислительная простота. Недостатки же -- явное задание количества кластеров, отсутствие гарантии нахождения оптимального решения.
На вход алгоритму задаются точки конечномерного признакового пространства и количество кластеров. Так же могут задаваться пороговые значения, использующиеся в критерии останова алгоритма: величина целевого функционала, скорость изменения функционала при выполнении основного цикла алгоритма. Целевой функционал алгоритма -- среднее внутрикластерное расстояние.
На первом шаге алгоритма производится инициализация центров кластеров. Это может быть осуществлено путём выбора случайных значений. Затем выполняется основной цикл алгоритма:
- Расчёт рядовой матрицы, сопоставляющей каждую точку признакового пространства с ближайшим к ней кластерному центру. Иначе это можно рассматривать как оптимизацию целевого функционала по принадлежности точек к кластерам
- Расчёт целевого функционала. Проверка условий критерия останова и выход из цикла в случае его выполнения.
- Пересчёт кластерных центров как центров масс точек кластера. Представляет собой оптимизацию целевого функционала по координатам центров кластеров.
Алгоритм гарантирует монотонное невозрастание целевой функционала, следовательно, сходимость к локальному минимуму.
Целевая функция алгоритма полностью аналогична целевой функции 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] 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)
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 по точкам. Так же заводятся счётчики: сколько точек принадлежит кластеру с данным центроидом. Затем покоординатно производится нормировка.
1.8 Ресурс параллелизма алгоритма
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 Программная реализация алгоритма
Для демонстрации работы алгоритма была создана собственная реализация алгоритма[5] с использованием фреймворка Intel® Threading Blocks [6]
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
1) Matlab, последовательная реализация
2) С++, последовательная реализация
3 Литература
- ↑ 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.
- ↑ Jan Jantzen (1998). Neurofuzzy Modelling.
- ↑ Нейский И. М. Классификация и сравнение методов кластеризации. // Интеллектуальные технологии и системы. Сборник учебно-методических работ и статей аспирантов и студентов. – М.: НОК «CLAIM», 2006. – Выпуск 8. – С. 130-142.
- ↑ Reduce paralle pattern [1]
- ↑ https://github.com/urandon/hcm_tbb
- ↑ https://www.threadingbuildingblocks.org/