Участник:Бротиковская Данута/Алгоритм k-means

Алгоритм k средних (k means)
Последовательный алгоритм
Последовательная сложность	[math]O(\Tau*k*n*d)[/math]
Объём входных данных	[math] n*d [/math]
Объём выходных данных	[math] n [/math]

Авторы страницы Данута Бротиковская и Денис Зобнин

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

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

Алгоритм k средних (k means) – неиерархический, итерационный метод кластеризации, получивший большую популярность благодаря своей простоте, наглядности реализации и достаточно высокому качеству работы. Был изобретен в 1950-х годах математиком Гуго Штейнгаузом и почти одновременно Стюартом Ллойдом. Особую популярность приобрел после публикации работы МакКуина в 1967.

Алгоритм представляет собой версию EM-алгоритма, применяемого также для разделения смеси гауссиан. Основная идея k means заключается в том, что данные произвольно разбиваются на кластеры, после чего итеративно перевычисляется центр масс для каждого кластера, полученного на предыдущем шаге, затем векторы разбиваются на кластеры вновь в соответствии с тем, какой из новых центров оказался ближе по выбранной метрике.

Цель алгоритма заключается в разделении [math]n[/math] наблюдений на [math]k[/math] кластеров таким образом, чтобы каждое наблюдение придележало ровно одному кластеру, расположенному на наименьшем расстоянии от наблюдения.

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

Дан набор из [math]n[/math] d-мерных векторов [math]X=\{\mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_n\}[/math]. Алгоритм k means разбивает набор [math]X[/math] на [math]k, k\lt =n[/math] наборов [math]S=\{S_1, S_2, ..., S_k\}, S_i \cap S_j= \varnothing, i \ne j, [/math] таким образом, чтобы минимизировать сумму квадратов расстояний от каждой точки кластера до его центра. Другими словами:

[math]\arg\min_{S} \sum\limits_{i=1}^k \sum\limits_{x \in S_i} \lVert \mathbf{x}- \mathbf{\mu}_i \rVert^2,[/math]

[math](1)[/math]

где [math]\mathbf{\mu}_i[/math] – центры кластеров, [math]i=\overline{1,k}[/math]

Шаги алгоритма:

1. Начальный шаг: инициализация кластеров

   Выбирается произвольное множество точек [math]\mu_i, i=\overline{1,k}[/math], рассматриваемых как начальные центры кластеров.

2. Распределение векторов по кластерам

   [math]\forall \mathbf{x}_i \in X, i=\overline{1,n}: \mathbf{x}_i \in S_j \iff j=\arg\min_{k}||\mathbf{x}_i-\mathbf{\mu}_k||^2[/math]

3. Пересчет центров кластеров

   [math] \forall i=\overline{1,k}: \widetilde{\mu_i} = \cfrac{1}{|S_i|}\sum_{x\in S_i}x[/math]

4. Проверка условия останова:

      if [math]\exist i\in \overline{1,k}: \mu_i \ne \widetilde{\mu_i}[/math] then
          goto 2;
      else
          stop
   

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

Вычислительным ядром являются шаги 2 и 3 приведенного выше алгоритма: распределение векторов по кластерам и пересчет центров кластеров.

Распределение векторов по кластерам предполагает вычисление расстояний между каждым вектором [math]\mathbf{x}_i \in X, i= \overline{1,n}[/math] и центрами кластера [math]\mathbf{\mu}_j, j= \overline{1,k}[/math]. Таким образом, данный шаг предполагает [math]k*n[/math] вычислений расстояний между d-мерными векторами.

Пересчет центров кластеров предполагает [math]k[/math] вычислений центров масс [math]\mathbf{\mu}_i[/math] множеств [math]S_i, i=\overline{1,k}[/math] представленных выражением в шаге 3 представленного выше алгоритма.

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

Инициализация центров масс [math]\mu_1, ..., \mu_k[/math].

Наиболее распространенными являются следующие стратегии:

• Метод Forgy
  В качестве начальных значений [math]\mu_1, ..., \mu_k[/math] берутся случайно выбранные векторы.
• Метод случайно разделения (Random Partitioning)
  Для каждого вектора [math]x_i \in X, i=\overline{1, n}[/math] выбирается случайным образом кластер [math]S_1, ..., S_k[/math], после чего для каждого полученного кластера вычисляются значения [math]\mu_1, ..., \mu_k[/math].

Распределение векторов по кластерам

Для этого шага алгоритма между векторами [math]x_i \in X, i=\overline{1, n}[/math] и центрами кластеров [math]\mu_1, ..., \mu_k[/math] вычисляется расстояние по формуле:

[math]\rho(\mathbf{v_1}, \mathbf{v_2}) = \lVert \mathbf{v_1}- \mathbf{v_2} \rVert^2[/math]

[math](2)[/math]

Пересчет центров кластеров

Для этого шага алгоритма производится пересчет центров кластера по формуле:

[math] \mu = \cfrac{1}{|S|}\sum_{x\in S}x[/math]

[math](3)[/math]

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

algorithm k-means is
    1. Инициализировать центры кластеров [math]\mathbf{\mu}_i^{(1)}, i=\overline{1,k}[/math]
    2. [math]t \leftarrow 1[/math]
    3. Распределение по кластерам

       [math]S_i^{(t)}=\{\mathbf{x}_p: \lVert\mathbf{x}_p-\mathbf{\mu}_i^{(t)}\rVert^2 \leq \lVert\mathbf{x}_p-\mathbf{\mu}_j^{(t)}\rVert^2 \quad \forall j=\overline{1,k}\},[/math]

       где каждый вектор [math]\mathbf{x}_p[/math] соотносится единственному кластеру [math]S^{(t)}[/math]
    4. Обновление центров кластеров

       [math]\mathbf{\mu}_i^{(t+1)} = \frac{1}{|S^{(t)}_i|} \sum_{\mathbf{x}_j \in S^{(t)}_i} \mathbf{x}_j [/math]

    5. if [math]\exist i \in \overline{1,k}: \mathbf{\mu}_i^{(t)} \ne \mathbf{\mu}_i^{(t+1)}[/math] then
           goto 3;
       else
           stop

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

Обозначим [math]\Theta_{centroid}^{d, m}[/math] временную сложность вычисления центорида кластера, число элементов которого равна [math]m[/math], в d-мерном пространстве.

Аналогично [math]\Theta_{distance}^d[/math] – временная сложность вычисления расстояния между двумя d-мерными векторами.

Сложность шага инициализации [math]k[/math] кластеров мощности [math]m[/math] в d-мерном пространстве – [math]\Theta_{init}^{k, d, m}[/math]

• Стратерия Forgy: вычисления не требуются, [math]\Theta_{init}^{k, d, m} = 0[/math]
• Стратегия случайного разбиения: вычисление центров [math]k[/math] кластеров, [math]\Theta_{init}^{k, d, m} = k*\Theta_{centroid}^{d, m}, m \le n[/math]

Cложность шага распределения d мерных векторов по [math]k[/math] кластерам – [math]\Theta_{distribute}^{k, d}[/math]

На этом шаге для каждого вектора [math]x_i \in X, i=\overline{1, n}[/math] вычисляется k расстояний до центров кластеров [math]\mu_1, ...\mu_k[/math]

[math]\Theta_{distribute}^{k, d} = n*k*\Theta_{distance}^d[/math]

Сложность шага пересчета центров [math]k[/math] кластеров размера [math]m[/math] в d-мерном пространстве – [math]\Theta_{recenter}^{k, d, m}[/math]

На этом шаге вычисляется [math]k[/math] центров кластеров [math]\mu_1, ...\mu_k[/math]

[math]\Theta_{recenter}^{k, d, m} = k*\Theta_{centroid}^{d, m}[/math]

Рассчитаем [math]\Theta_{centroid}^{d, m}[/math] для кластера, число элементов которого равно [math]m[/math]

[math]\Theta_{centroid}^{d, m}[/math] = [math]m * d[/math] сложений + [math]d[/math] делений

Рассчитаем [math]\Theta_{distance}^d[/math] в соответствие с формулой [math](2)[/math]

[math]\Theta_{distance}^d[/math] = [math]d[/math] вычитаний + [math]d[/math] умножений + [math](d-1)[/math] сложение

Предположим, что алгоритм сошелся за [math]\Tau[/math] итераций, тогда временная сложность алгоритма [math]\Theta_{k-means}^{d, n}[/math]

[math]\Theta_{k-means}^{d, n} \le \Theta_{init}^{k, d, n} + \Tau*(\Theta_{distribute}^{k, d} + \Theta_{recenter}^{k, d, n})[/math]

Операции сложения/вычитания:

[math]\Theta_{k-means}^{d, n} \le k*n*d+ \Tau*(k*n*(2d-1) + k*n*d) = k*n*d+ \Tau*(k*n*(3d-1)) \thicksim O(\Tau*k*d*n)[/math]

Операции умножения/деления:

[math]\Theta_{k-means}^{d, n} \le k*d + \Tau*(k*n*d + k*d) = k*d + \Tau * k*d*(n+1) \thicksim O(\Tau*k*d*n) [/math]

Получаем, что временная сложность алгоритма k-means кластеризации [math]n[/math] d-мерных векторов на [math]k[/math] кластеров за [math]\Tau[/math] итераций:

[math] \Theta_{k-means}^{d, n} \thicksim O(\Tau*k*d*n) [/math]

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

Рис. 1. Схема распределения векторов по кластерам

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

В последовательной реализации алгоритма наиболее вычислительно сложным шагом является распределение данных по кластерам – он требует [math]n*k[/math] вычислений расстояний до центров кластеров (см. раздел «Последовательная сложность алгоритма»). При этом определение соответствующего кластера для двух разных объектов [math]\mathbf{x}_i[/math] и [math]\mathbf{x}_j[/math] может производиться независимо.

При наличии [math]P[/math] вычислительных потоков этап распределения данных по кластерам может быть выполнен параллельно – по [math]\frac{n}{P}[/math] объектов в каждом потоке (поток [math]thread_i, i = \overline{1..P}[/math] обрабатывает подмножество объектов [math]\{\mathbf{x}_{\frac{(i-1)*n}{P}}, ..., \mathbf{x}_{\frac{i*n}{P}}\}[/math]). При этом каждому вычислительному потоку передаются координаты центров всех кластеров [math]\mu_1, ..., \mu_k[/math].

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

Работа алгоритма состоит из [math]\Tau[/math] итераций, в каждой из которых происходит вычисление расстояний между всеми векторами и центрами кластеров + перерасчет центров масс. Вычислим параллельную сложность [math]\Psi_*[/math] каждого из шагов, а также параллельную сложность всего алгоритма, [math]\Psi_{k-means}[/math]. Будем исходить из предположения, что может быть использовано любое необходимое число потоков.

Вычисление расстояний между всеми векторами и центрами кластеров.

Поскольку на данном шаге для каждой пары векторов [math]x_i, i=\overline{1,n}[/math] и [math]\mu_j, j=\overline{1,k}[/math] операции вычисления расстояния не зависят друг от друга, они могут выполняться параллельно. Тогда, разделив все вычисление расстояний на n потоков, получим, что в каждом потоке будет выполняться только одна операция вычисления расстояния между векторами размерности d. Таким образом, параллельная сложность данного шага определяется сложностью параллельной операции вычисления расстояния между d-мерными векторами, [math]\Psi_{distance}^d[/math] и сложностью определения наиболее близкого кластера (паралельное взятие минимума по расстояниям), [math]\Psi_{min}^k[/math]. Для оценки [math]\Psi_{distance}^d[/math] воспользуемся параллельной реализацией Нахождение частных сумм элементов массива сдваиванием. Аналогично, [math]\Psi_{min}^k = log(k)[/math]. В результате, [math]\Psi_{distance}^d = O(log(d))[/math]. Таким образом,

[math]\Psi_{distribute} = \Psi_{distance}^d + \Psi_{min}^k = O(log(d))+O(log(k)) = O(log(kd))[/math]

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

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

• Целое положительное число [math]k[/math] – количество кластеров.
• Матрица из [math]n*d[/math] элементов – координат векторов (наблюдений).

Объем входных данных
[math]1[/math] целое число + [math]n*d[/math] вещественных чисел (при условии, что координаты – вещественные числа).

Выходные данные
Целые положительные числа – номера кластеров, соотвествующие каждому вектору (при условии, что нумерация кластеров начинается с 1).

Объем выходных данных
[math]n[/math] целых положительных чисел.

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

Вычислительная мощность

Вычислительная мощность алгоритма k means равна [math]\frac{k*d*n*\Tau}{n*d} = k*\Tau [/math], где [math]k[/math] – число кластеров, [math]\Tau[/math] – число итераций алгоритма.

Детерминированность и Устойчивость

Алгоритм k means является итерационным. Количество итераций алгоритма в общем случае не фиксируется и зависит от начального расположения объектов в пространстве, параметра [math]k[/math], а также от начального приближения центров кластеров, [math]\mu_1, ..., \mu_k[/math]. В результате этого может варьироваться результат работы алгоритма. По этим причинам алгоритм не является ни детермирированным, ни устойчивым.

Сильные стороны алгоритма:

• Сравнительно высокая эффективность при простоте реализации
• Высокое качество кластеризации
• Возможность распараллеливания
• Существование множества модификаций

Недостатки алгоритма:

• Количество кластеров является параметром алгоритма

• Чувствительность к начальным условиям

  Инициализация центров кластеров в значительной степени влияет на результат кластеризации.

• Чувствительность к выбросам и шумам

  Выбросы, далекие от центров настоящих кластеров, все равно учитываются при вычислении их центров.

• Возможность сходимости к локальному оптимуму

  Итеративный подход не дает гарантии сходимости к оптимальному решению.

• Использование понятия "среднего"

  Алгоритм неприменим к данным, для которых не определено понятие "среднего", например, категориальным данным.

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

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

2.7.1 Открытое программное обеспечение

1. CrimeStat
Программное обеспечение, созданное для операционных систем Windows, предоставляющее инструменты статистического и пространственного анализа для решения задачи картирования преступности.
2. Julia
Высокоуровневый высокопроизводительный свободный язык программирования с динамической типизацией, созданный для математических вычислений, содержит реализацию k-means.
3. Mahout
Apache Mahout - Java библиотека для работы с алгоритмами машинного обучения с использованием MapReduce. Содержит реализацию k-means.
4. Octave
Написанная на C++ свободная система для математических вычислений, использующая совместимый с MATLAB язык высокого уровня, содержит реализацию k-means.
5. Spark
Распределенная реализация k-means содержится в библиотеке Mlib для работы с алгоритмами машинного обучения, взаимодействующая с Python библиотекой NumPy и библиотека R.
6. Torch
MATLAB-подобная библиотека для языка программирования Lua с открытым исходным кодом, предоставляет большое количество алгоритмов для глубинного обучения и научных расчётов. Ядро написано на Си, прикладная часть выполняется на LuaJIT, поддерживается распараллеливание вычислений средствами CUDA и OpenMP. Существуют реализации k-means.
7. Weka
Cвободное программное обеспечение для анализа данных, написанное на Java. Содержит k-means и x-means.
8. Accord.NET
C# реализация алгоритмов k-means, k-means++, k-modes.
9. OpenCV
Написанная на С++ библиотека, направленная в основном на решение задач компьютерного зрения. Содержит реализацию k-means.
10. MLPACK
Масштабируемая С++ библиотека для работы с алгоритмами машинного обучения, содержит реализацию k-means.
11. SciPy
Библиотека Python, содержит множество реализаций k-means.
12. scikit-learn
Библиотека Python, содержит множество реализаций k-means.
13. R
Язык программирования для статистической обработки данных и работы с графикой, а также свободная программная среда вычислений с открытым исходным кодом в рамках проекта GNU, содержит три реализации k-means.
14. ELKI
Java фреймворк, содержащий реализацию k-means, а также множество других алгоритмов кластеризации.

2.7.2 Проприетарное программное обеспечение

1. Ayasdi
2. Stata
3. Mathematica
4. MATLAB
5. SAS
6. RapidMiner
7. SAP HANA