Просмотр кода страницы Участник:Светлана Лукьяненко/Строгий алгоритм С средних (Hard C-Means, HCM)

{{assignment|Coctic|Chalker}}

{{algorithm
| name              = Строгий алгоритм С средних (Hard C-Means,HCM)
| serial_complexity = <math>O(nNCi)</math>
| pf_height         = <math>O(ilog_2(NC))</math>
| pf_width          = <math>O(nNC)</math>
| input_data        = <math>nN+1</math>
| output_data       = <math>C(N+n)</math>
}}


Авторы описания: [[Участник:Светлана_Лукьяненко|Лукьяненко Светлана]](1.1, 1.3, 1.4, 1.8, 1.10, 2.4, 2.7), [[Участник:Юрий_Комаров|Комаров Юрий]](1.2, 1.5, 1.6, 1.7, 1.8, 1.9, 2.4)

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

=== Общее описание алгоритма ===
'''Строгий алгоритм C средних''' ('''''Hard C-Means'''''; ''Jang, Sun and Mizutani, 1997''<ref>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.</ref>) позволяет распределить большие наборы числовых данных по кластерам в многомерном пространстве<ref>Jan Jantzen (1998). [http://faculty.petra.ac.id/resmana/private/fuzzy/nfmod.pdf Neurofuzzy Modelling.]</ref>. Алгоритм является управляемым, так как для работы с ним требуется задание одного параметра — числа кластеров, на которое необходимо разделить данные. 
Входными данными для алгоритма являются набор векторов и количество кластеров, также может задаваться пороговое значение величины объектной функции, по достижении которого расчёт прекращается. Выходными данными для алгоритма является рядовая матрица, отображающая принадлежность исходных элементов построенным кластерам, и центры построенных кластеров.

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

# Расчёт рядовой матрицы, сопоставляющей каждый вектор с ближайшим кластерным центром.
# Расчёт объектной функции; проверка условий остановки и завершение расчёта в случае их выполнения.
# Пересчёт кластерных центров.

Так как алгоритм является итеративным и зависит от начальных значений кластерных центров, он может сходиться к локальному минимуму.

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

==== Входные и выходные данные ====
На вход подаётся набор векторов <math>\{u_k\}_{k=1}^{N}</math> пространства <math>\mathbb{R}^n</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\}_{i=1}^{C}</math> центров этих кластеров.

==== Описание алгоритма ====
Алгоритм состоит из 5 последовательных шагов.<ref>Нейский И. М. [http://it-claim.ru/Persons/Neyskiy/Article2_Neiskiy.pdf Классификация и сравнение методов кластеризации.] // Интеллектуальные технологии и системы. Сборник учебно-методических работ и статей аспирантов и студентов. – М.: НОК «CLAIM», 2006. – Выпуск 8. – С. 130-142.</ref>

'''<u>Шаг 1</u>'''

'''Инициализация кластерных центров''' <math>c_i \, (i=1,2,\ldots,C)</math> может происходить за счёт случайного выбора необходимого числа точек из входного набора либо путём явного указания.

'''<u>Шаг 2</u>'''

'''Вычисление рядовой матрицы''' <math>M = \|m_{ik}\|</math>, элементы которой задаются следующими выражениями (<math>N</math> — количество элементов во входном наборе данных):
:<math>
m_{ik} = \begin{cases}
1, & \|u_k-c_i\|^2 \leqslant \|u_k-c_j\|^2 \;\; \text{for all } i \neq j; \;\; (i,j=1,2,\ldots,C; \; k = 1,2,\ldots,N)\\
0, & \text{otherwise}.
\end{cases}
</math>
<!-- При этом, как легко видеть, матрица <math>M</math> обладает следующими свойствами:
:<math>
\sum\limits_{i=1}^C m_{ij} = 1, \;\; \sum\limits_{i=1}^C\sum\limits_{j=1}^K m_{ij} = K.
</math> -->

Будем говорить, что <math>u_k \in C_i</math>, если <math>m_{ik} = 1</math>.

'''<u>Шаг 3</u>'''

'''Расчёт объектной функции:'''
:<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>

На этом шаге происходит '''<u>остановка и выход из цикла</u>''', если полученное значение ниже пороговой величины <math>\varepsilon_{tr}</math> или полученное значение не сильно отличается от значений, полученных на предыдущих циклах.

'''<u>Шаг 4</u>'''

'''Пересчёт кластерных центров''' выполняется в соответствии со следующим уравнением:
:<math>
c_i = \frac{1}{|C_i|} \sum\limits_{k:\, u_k \in C_i} u_k,
</math>
где под <math>|C_i|</math> подразумевается количество элементов в <math>i</math>-ом кластере.

'''<u>Шаг 5</u>'''

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

=== Вычислительное ядро алгоритма ===
Вычислительное ядро алгоритма составляют вычисления квадратов нормы расстояния от каждого из векторов исходной выборки до промежуточных кластерных центров <math>\|u_k-c_i\|^2</math>. На каждой итерации алгоритма происходит расчёт <math>CN</math> таких расстояний.

=== Макроструктура алгоритма ===
Основной операцией, выполняемой в алгоритме, как сказано выше, является операция нахождения квадратов расстояний между исходными векторами и кластерными центрами. 

<math>
d_{ki} = \|u_k-c_i\|^2
</math>

Также на шаге вычисления объектной функции производится множественное суммирование полученных расстояний:

<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>

а на шаге расчета кластерных центров происходит суммирование векторов:

<math>
c_i = \frac{1}{|C_i|} \sum\limits_{k:\, u_k \in C_i} u_k.
</math>

=== Схема реализации последовательного алгоритма ===
Реализация последовательного алгоритма состоит в итеративном повторении блока, описанного далее в форме псевдокода, до тех пор, пока не будет достигнута желаемая точность.

<source lang="pascal">
J = 0;
M = zeros(N,C);
for i = 1 to N //go through every element
   d_min = MAXINT;
   j_min = 1;
   for j = 1 to C //deside which claster does element belong to
      d_ij  = sq_norm(u[i] - c[j]);
      if d_ij < d_min
         d_min = d_ij;
         j_min = j;
   //generation of output data
   J += d_min; 
   M[i,j_min] = 1;
   c_new[j] += u[i];
   c_newCount[k]++;

</source>

=== Последовательная сложность алгоритма ===
На каждой итерации алгоритма требуется выполнение:
* <math>NC</math> операций вычисления квадрата нормы разности элементов;
* <math>NC</math> операций сравнения (вычитаний);
* <math>3N</math> операций сложения.

Каждая операция вычисления нормы разности векторов предполагает выполнение <math>2n-1</math> арифметических операций сложения и вычитания и <math>n</math> операций возведения в квадрат (умножений), поэтому получаем для каждой итерации:
* <math>(2nC + 3)N</math> арифметических операций сложения и вычитания;
* <math>nNC</math> операций умножения.

Ожидая, что желаемая точность достигается за <math>i</math> итераций, получаем итоговую сложность последовательной реализации алгоритма <math>O(nNCi)</math>.

Отметим, что ввиду сделанных замечаний, построить сколько-нибудь точную оценку числа итераций <math>i</math> не представляется возможным.

=== Информационный граф ===
Приведём аналитическое и графическое описание графа алгоритма.

Граф алгоритма состоит из 5 различных видов вершин.

'''Первая''' группа вершин отмечена символом '''''d''''', ей соответствует вычисление минимального расстояния до кластерных центров, включающее в себя следующую последовательность действий:
* независимое вычисление расстояний от заданного элемента <math>u_i</math> до каждого из кластерных центров;
* последовательное сравнение полученных значений с целью отыскания минимума;
* передача найденного минимального расстояния и соответствующего ему индекса кластера.

'''Вторая''' группа вершин отмечена символом '''''M''''', ей соответствует операция записи единички в столбец, отвечающий текущей ветке графа, в позицию, полученную на предыдущем шаге и соответствующую ближайшему кластерному центру.

'''Третья''' вершина отмечена символом '''''J''''', подразумевается, что на этом этапе происходит расчёт объектной функции.

'''Четвертая''' вершина обозначена треугольным маркером и помечена '''''if''''', на этом этапе происходит проверка условия достижения необходимой точности (условия остановки итерационного процесса).

'''Пятая''' соответствует операции пересчёта кластерных центров, она отмечена как '''''{c}'''''.

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

[[File:HardCMeans1.png|thumb|800px|center|Рисунок 1. Информационный граф алгоритма. Здесь '''''in''''' — ввод входных данных (<math>C</math> стартовых значений кластерных центров и 1 значение <math>u_i</math> на каждую вершину графа), '''''d''''' — вычисление минимального расстояния до кластерных центров (подробнее см. рисунок 2), '''''M''''' — запись единички в соответствующий столбец матрицы <math>M</math>, '''''J''''' — вычисление объектной функции по результатам полученных вычислений, '''''if''''' — проверка условия достижения заданной точности, '''''{c}''''' — пересчёт кластерных центров и '''''out''''' — вывод полученных данных (матрица <math>M</math> и <math>C</math> значений полученных кластерных центров).]]
[[File:HardCMeans2.PNG|thumb|400px|center|Рисунок 2. Подробное описание работы блока  '''''d'''''. Здесь '''''in''''' — ввод входных данных (по 1 значению центра кластера на каждую из <math>C</math> соответствующих вершин графа и значение <math>u_i</math>, единое для этого '''''d'''''-блока), '''''d1,...,dC''''' — прямое вычисление  расстояния до кластерных центров, '''''min''''' — вычисление минимума двух чисел (с сохранением соответствующего минимуму кластерного индекса) и '''''out''''' — вывод полученных данных (индекс кластерного центра, расположенного ближе всего к точке <math>u_i</math>, и соответствующее расстояние).]]

=== Ресурс параллелизма алгоритма ===
Для выполнения одной итерации алгоритма C средних необходимо выполнить параллельно ряд вычислений для каждого из <math>N</math> элементов:

# Вычисление <math>C</math> расстояний между элементом и кластерными центрами.
# Выполнение <math>(C - 1)</math> сравнений для определения минимального.
# Определение значений соответствующего рассматриваемому элементу <math>u_i</math> столбца матрицы <math>M</math>.

Эта последовательность операций может быть выполнена для каждого элемента <math>u_i</math> входной выборки независимо и поэтому составляет основной ресурс параллелизма алгоритма.

При этом, как можно видеть из информационного графа алгоритма, часть вычислений проводится в последовательном режиме, например, расчёт значения объектной функции <math>J</math>. Аналогичная ситуация имеет место и с пересчётом кластерных центров — разделив эту операцию на <math>C</math> независимых процессов, мы потеряем в производительности за счёт увеличения накладных расходов, связанных с пересылкой данных и координацией вычислений.

Стоит отметить, что при малом объёме выборки и размерности системы попытка использования параллельной реализации алгоритма приведёт к понижению производительности, поскольку проведение последовательных расчётов потребует меньше времени, нежели пересылка данных (хотя бы матриц размера <math>nC</math>, определяющих кластерные центры).

Ширина [[глоссарий#Ярусно-параллельная форма графа алгоритма|ЯПФ]] составляет <math>O(nNC)</math> (<math>N</math> блоков <math>d</math>, в каждом из которых независимо выполняются <math>C</math> операций вычислений квадратов норм, каждая из которых предполагает ширину не более <math>n</math>). 

Высота одной итерации алгоритма оценивается через <math>log_2(NC)</math> (<math>N</math> сложений для расчёта <math>J</math> можно провести за <math>log_2(N)</math> операций, сложность блока <math>d</math> равна <math>O(log_2(C))</math>). Полагая, что число итераций равно <math>i</math>, получаем итоговую высоту ярусно-параллельной формы алгоритма: <math>ilog_2(NC)</math>

=== Входные и выходные данные алгоритма ===
'''Входные данные:''' 
*<math>\{u_k\}_{k=1}^{N} \subset \mathbb{R}^n</math> — исходный набор векторов;
*<math>C</math> — количество кластеров (<math>2 \leqslant C \leqslant N</math>);
* ''(опционально)'' <math>\varepsilon_{tr}</math> — пороговая точность работы алгоритма.

'''Объём входных данных:''' <math>nN+1</math>.

'''Выходные данные:'''
*<math>M</math> — рядовая матрица из пространства <math>\mathbb{R}^{C \times N}</math>.
*<math>\{c_k\}_{k=1}^{C}</math> — набор построенных кластерных центров.

'''Объём выходных данных:''' <math>C(N+n)</math>.

=== Свойства алгоритма ===
'''Назначение:''' кластеризация больших наборов числовых данных.

'''Достоинства:'''
* лёгкость реализации;
* вычислительная простота.

'''Недостатки:'''
* нет гарантии, что итерационный процесс сойдётся к оптимальному решению;
* производительность существенно зависит от выбора начальных кластерных центров.

== Программная реализация алгоритма ==
=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
==== Локальность реализации алгоритма ====
===== Структура обращений в память и качественная оценка локальности =====
===== Количественная оценка локальности =====
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
Исследуем масштабируемость параллельной реализации строгого алгоритма C средних, соответствующей последовательной реализации, описанной в [[#Схема реализации последовательного алгоритма|п1.5]]. Для проведения расчётов использовались ресурсы вычислительного комплекса [http://hpc.cmc.msu.ru/bgp IBM Blue Gene/P] и [https://github.com/ykomarov94/AlgoWiki собственная реализация алгоритма].

<u>Основные характеристики реализации алгоритма и запуска программы:</u>
* программа написана на языке C с использованием [http://hpc.cmc.msu.ru/bgp/development/parallel библиотеки MPI];
* сборка программы произведена с использованием компилятора [http://hpc.cmc.msu.ru/bgp/development/compilers IBM XL C Advanced Edition for Blue Gene/P V9.0] версии 9.0.0.0005;
* при компиляции не использовались дополнительные опции;
* при тестировании использовались псевдослучайные наборы векторов заданной размерности;
* для каждой конфигурации проводилась серия экспериментов, после чего для получения конечных результатов производительности проводилось усреднение.

<u>Набор и границы значений изменяемых параметров запуска реализации алгоритма:</u>

* число процессоров <math>[16 : 256]</math> с шагом по степеням двойки;
* объём выборки <math>[10^3:10^7]</math>  с шагом по степеням десятки;
* размерность системы <math>n=10</math>;
* количество кластеров <math>C = 11</math>.

В качестве изменяемого параметра для построения графика был выбран объём выборки <math>N</math>, поскольку варьирование прочих параметров существенно хуже отражает возможности применения технологии параллельных вычислений. Это связано с тем, что при малом количестве векторов время, затраченное на пересылку данных, сравнимо со временем, затрачиваемым непосредственно на расчёты; это хорошо видно на графике: при объёмах выборки, составляющих <math>10^3</math> и <math>10^4</math> векторов, использование большего числа процессоров не приводит к приросту производительности (наоборот, производительность падает) ввиду высоких накладных расходов.

[[file:HardCMeansPerfomance.png|thumb|center|700px|Рис.3 Реализация строгого алгоритма C средних. График производительности (GFLOPS) в зависимости от числа процессоров и объёма выборки.]]

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===

1)  [https://github.com/artpar/Hard-Clustering-implementation Последовательная реализация на С++ ]

2) [http://projects.jasmeet.me/2013/03/HardCMeans.html Последовательная реализация на MatLab ]

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

== Литература ==
<references />