Участник:Astandrik/Алгоритм кластеризации, основанный на минимальном остовном дереве: различия между версиями
Перейти к навигации
Перейти к поиску
Astandrik (обсуждение | вклад) (Создание шаблона страницы) |
Zhum (обсуждение | вклад) |
||
| (не показано 105 промежуточных версий 3 участников) | |||
| Строка 1: | Строка 1: | ||
| + | {{Assignment|Zhum}} | ||
| + | |||
{{algorithm | {{algorithm | ||
| − | | name = | + | | name = Алгоритм кластеризации, основанный на минимальном остовном дереве |
| − | | serial_complexity = <math>O( | + | | serial_complexity = <math> O(N^2 log(N)) </math> |
| − | | pf_height = <math>O( | + | | pf_height = <math> O(log(N)) </math> |
| − | | pf_width = <math>O( | + | | pf_width = <math> O(N) </math> |
| − | | input_data = <math> | + | | input_data = <math> N </math> |
| − | | output_data = <math> | + | | output_data = <math> N </math> |
}} | }} | ||
| − | Основные авторы описания: [[Участник: | + | |
| + | Основные авторы описания: [[Участник:Astandrik|Стандрик Антон Сергеевич]] | ||
== Свойства и структура алгоритма == | == Свойства и структура алгоритма == | ||
| Строка 14: | Строка 17: | ||
=== Общее описание алгоритма === | === Общее описание алгоритма === | ||
| − | + | Данный алгоритм предназначен для решения задач кластеризации. Говоря в общем, кластеризация — многомерная статистическая процедура, выполняющая сбор данных, содержащих информацию о выборке объектов, и затем упорядочивающая объекты в сравнительно однородные группы <ref>Классификация и кластер. Под ред. Дж. Вэн Райзина. М.: Мир, 1980. 390 с.</ref>. | |
| + | |||
| + | Впервые идея использования минимального остовного дерева для кластеризации данных была преждожена C. T. Zah в 1971 году<ref name=zakh>C. T. Zahn, “graph-Theriotical methos for detecting and describing getalt clusters,” IEEE Trans. Computers, vol. 20, no. 1, pp. 68-86, 1971 </ref>. Этот метод широко используется исследователями при анализе изображений <ref>Y. Xu, V. Olman and E. C. Uberbacher, “A mentation using minimum spanning trees,” Image and Vission Computing, vol. 15, pp. 47-57,1997.</ref><ref>Zhi Min Wang, Yeng Chai Soh, Qing Song, Kang Sim, “Adaptive spatial information-theoretic clustering for image segmentation,” Pattern Recognition, vol. 42, no. 9, pp. 2029-2044, 2009.</ref>, а также при анализе биологических данных <ref>Y. Xu, V. Olman and D. Xu, “Clustering gene expression data using a graph-Theriotic approach: An application of minimum spanning trees,” Bioinformatics, vol. 18, no. 4, pp. 536-545, 2002.</ref>. | ||
| + | |||
| + | В общем описании работы алгоритма можно выделить следующие шаги: | ||
| − | + | # Разбиение входных векторов на пары и вычисление расстояний в них. Для этого небходим выбор метрики пространства, в котором осуществляется кластерный анализ. В качестве примеров метрик можно привести: | |
| + | #* Метрика Евклида: <math>d_{ij} = \sqrt{\sum^n_{k=1}(x_{ik} - x_{jk})^2}</math> | ||
| + | #* Метрика city-block (манхэттенская метрика): <math>d_{ij} = \sum^n_{k=1}|x_{ik} - x_{jk}| </math> | ||
| + | #* Метрика Брея-Картиса: <math>d_{ij} = \sum^n_{k=1}\frac{|x_{ik} - x_{jk}|}{|x_{ik} + x_{jk}| }</math> | ||
| + | # Построение минимального остовного дерева по заданному множеству. Существует три основных алгоритма [[Построение минимального остовного дерева (MST) | построения минимального остовного дерева]]: | ||
| + | #* [[Алгоритм Прима]] | ||
| + | #* [[Алгоритм Крускала]] | ||
| + | #* [[Алгоритм Борувки]] | ||
| + | # Разделение полученного дерева на несвязное множество поддеревьев в соответствии с некоторой стратегией. Можно выделить три популярные стратегии кластеризации: | ||
| + | #* SEMST<ref>Asano, M.K.T., Bhattacharya, B., Yao, F.: Clustering algorithms based on minimum and maximum spanning trees. In: In Proceedings of the fourth annual symposium on Computational Geometry. (1998) 252–257</ref> (standard Euclidean minimum spanning tree clustering algorithm) - предполагается, что нам заранее известно число кластеров. В таком случае, если число кластеров равно <math>k</math>, то из полученного в первом пункте остовного дерева удаляются <math>k-1</math> связей с наибольшим весом. | ||
| + | #* CEMST<ref>Ye, B., Chao, K.M.: Spanning Trees and Optimization Problems. Chapman and | ||
| + | Hall (2004)</ref>(The maximum cost minimum spanning tree) - Работает так же как и SEMST, однако вместо веса вершины используется её стоимость, которая вычисляется следующим образом: Пусть ребро <math>e = (x_{i}, x_{j})</math>, тогда <math>cost(e) = w(e) * degree(x_{i}) * degree(x_{j})</math>. Далее, зная число кластеров <math>k</math> и стоимости рёбер, удаляются <math>k-1</math> рёбер с максимальной стоимостью. | ||
| + | #* ZEMST<ref name=zakh /> (The Zahn’s minimum spanning tree) - от двух предыдущих эта стратегия отличается тем, что здесь число кластеров не известно заранее, а вычисляется в процессе кластеризации. Его суть в том, что удаляются те рёбра, веса которых значительно превосходят средний вес ребра в соседних поддеревьях. | ||
| − | + | ==== Достоинства ==== | |
| − | + | Основным достоинством алгоритма является его способность выделять кластеры произвольной формы. На следующем примере можно видеть, как алгоритм даже в случае SEMST алгоритма корректно разбивает вложенные окружности на три кластера: | |
| − | + | [[File:Selection_117_.png|frame|center|alt=Point-connected clusters| Слева направо: Рисунок 1: Начальные данные. Рисунок 2: построенное минимальное остовное дерево. Рисунок 3: разбиение на кластеры путём удаления самых длинных связей.]] | |
| − | + | ==== Недостатки ==== | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | Недостатком является трудность определения четких критериев, при которых осуществляется разбиение данных на разные кластеры. | |
| − | + | Кроме того, существуют наборы данных, на которых популярные алгоритмы (SEMST, CEMST, ZEMST) не дают правильного результата, и для которых требуется применение более сложных алгоритмов, к числу которых можно отнести, например, алгоритм HEMST <ref>Clustering with Minimum Spanning Trees Yan Zhou, Oleksandr Grygorash, Thomas F. Hain, Meador Inge Zach Jorgensen School of Computer and Information Sciences University of South Alabama, Mobile, AL 36688 USA Urban Insight Inc. 5657 Wilshire Blvd Ste 290, Los Angeles, CA 90036</ref> с временем работы <math>O(n^2 * log(n))</math>. | |
| − | |||
| − | + | В качестве иллюстрации можно привести два похожих визуально, но различных с точки зрения разбиения на кластеры с помощью описываемого алгоритма, примера. (В приведенных примерах первый рисунок обозначает начальный набор данных. Второй - построенное на основе данных минимальное остовное дерево. Третий рисунок - разбиение на кластеры.): | |
| − | + | [[File:Selection_116_.png|frame|center|alt=Point-connected clusters| Правильное разбиение на два кластера.]] | |
| − | = | + | [[File:Selection_115_.png|frame|center|alt=Point-connected clusters| Неверное разбиение на один кластер.]] |
| − | + | Отличие второго примера - в наличии плотной цепочки связей между двумя кластерами, которая и приводит к ошибке наивного алгоритма удаления самых длинных связей, и таким образом здесь потребуется более сложный алгоритм с более значительным временем работы для правильного нахождения кластеров. | |
| − | + | === Математическое описание алгоритма === | |
| − | + | Исходные данные: набор векторов размерности <math>n</math> : <math>(a_1, a_2, ..., a_n)</math>, где <math> a_j \in \mathbb {R} \quad \forall j \in [1,n] </math>. | |
| − | :<math> | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | </math> | ||
| − | + | Вычисляемые данные: Сгруппированные по кластерам данные. | |
| − | + | ==== Основные обозначения ==== | |
| − | = | + | * <math>E</math> - множество рёбер. |
| + | * <math>V</math> - множество вершин. | ||
| + | * <math>G = (V,E) </math> - неориентированный граф, заданный множеством вершин <math>V</math> и множеством рёбер <math>E</math> | ||
| + | * <math>MST(G)</math> (Minimum Spanning Tree) - минимальное остовное дерево графа <math></math>. Остовное дерево — ациклический связный подграф данного связного неориентированного графа, в который входят все его вершины. Минимальное остовное дерево в связанном взвешенном неориентированном графе — это остовное дерево этого графа, имеющее минимальный возможный вес, где под весом дерева понимается сумма весов входящих в него рёбер. | ||
| + | * <math>MSF(G)</math> (Minimum Spanning Forest) - набор минимальных остовных деревьев для всех компонент связности графа <math> G </math> | ||
| − | + | В описании алгоритмов этой статьи часто используется термин "построение MST", подразумевая, что на самом деле в результате построения может получиться как одна компонента связности (MST), так и множество компонент связности (MSF). | |
| − | :<math>\ | + | ==== Алгоритмы кластеризации ==== |
| + | '''SEMST'''(standard Euclidean minimum spanning tree) - простейший алгоритм, основанный на заранее известном числе кластеров <math> k </math>. Путём удаления <math> k - 1 </math> связей с наибольшим весом предполагается разбиение остовного дерева на нужное количество кластеров. Иными словами минимальное остовное дерево разбивается на <math> k </math> поддеревьев: <math> G = \{G_1,G_2,...,G_k\} </math> путём минимизации суммы всех весов: <math> W = \sum ^k_{l=1} \sum _{x_i,x_j \in V_l} w_{x_i,x_j} -> min </math>, где <math> V_l </math> - множество вершин в поддереве <math> G_l </math>. В конце такой процедуры каждое множество <math> V_l </math> становится кластером <math> C_l </math>. Поскольку (в случае использования алгоритма Крускала для построения остовного дерева) рёбра перед началом кластеризации уже отсортированы, для удаления наибольших <math> k </math> рёбер требуется <math> O(k) </math> операций, и для распределения получившихся вершин по кластерам требуется время <math> O(n) </math>.''' | ||
| − | в | + | '''CEMST'''(maximum cost Euclidean minimum spanning tree) - Работает в точности, как и SEMST, но здесь вместо весов рёбер используется т.н. "стоимость", которая высчитывается следующим образом: если ребро <math> (x_i, x_j) </math> соединяет вершины <math> x_i </math> и <math> x_j </math>, то его "стоимость" считается как: <math> Cost_{x_i, x_j} = w_{x_i,x_j} * Deg(x_i) * Deg(x_j) </math>, где <math> Deg(x_i) </math> - степень вершины <math> x_i </math> (т.е. количество рёбер графа G, инцидентных вершине <math> x_i </math>). В остальном алгоритм действует точно так же, как и SEMST, удаляя <math> k-1 </math> рёбер с наибольшей "стоимостью". Сложность <math> O(k) </math> . Алгоритм направлен на удаление рёбер, соединяющих высоко нагруженные (с точки зрения количества инцидентных рёбер) вершины. |
| − | :<math> | + | '''ZEMST''' (The Zahn’s maximum spanning tree (Zahn, 1971)) <ref> Attribute Clustering Based on Heuristic Tree Partition Jorge Cordero H. and Yifeng Zeng Department of Computer Science Aalborg University 9220 Aalborg, Denmark</ref> - отличается от двух предыдущих алгоритмов в первую очередь тем, что здесь количество кластеров заранее не известно, но определяется исходя из данных. Данный алгоритм не является жадным, в отличие от первых двух, так как в нём рассматривается не только само ребро и его вес, а также и поддеревья, которые соединены этим ребром. То есть для ребра <math> (x_i, x_j) </math> рассматриваются поддеревья <math> N_i = (V_{N_i}, E_{N_i}), N_j = (V_{N_j}, E_{N_j}) </math>, где <math> N_i </math> - поддерево глубины(высоты) <math> d </math>(величина d определяется эмпирически), состоящее из всех достижимых из вершины <math> x_j </math> вершин за <math> d </math> шагов через все достижимые из неё рёбра (исключая само ребро <math> (x_i, x_j) </math>). |
| + | Введём обозначения: | ||
| + | *<math>\overline{w_{N_i}} = \frac{\sum_{(x_r, x_s) \in E_{N_j}}w_{x_r,x_s}}{|E_{N_j}|} </math> - средний вес ребра в поддереве <math> N_i </math> | ||
| + | *<math>\sigma_{N_i} = \sqrt{\frac{\sum_{(x_r, x_s) \in E_{N_j}}(w_{x_r,x_s} - \overline{w_{N_i}})^2}{|E_{N_j}|}} </math> - стандартное отклонение веса ребра в поддереве <math> N_i </math> | ||
| + | Ребро <math> (x_i, x_j) </math> удаляется из списка в случае, если выполняется хотя бы одно из двух условий: | ||
| + | #<math> w_{x_i,x_j} > \overline{w_{N_i}} - \sigma_{N_i} </math> | ||
| + | #<math> w_{x_i,x_j} > \overline{w_{N_j}} - \sigma_{N_j} </math> | ||
| + | #<math> w_{x_i,x_j} > \frac{\sum_{(x_r, x_s) \in E_i \cup E_j}w_{x_r,x_s}}{|E_i \cup E_j|} </math> | ||
| − | в | + | Сложность этого алгоритма составляет в худшем случае <math> O(N^2) </math> (поскольку может оказаться нужным перебрать все рёбра), а в лучшем <math> O(Nd) </math>. |
| − | + | === Вычислительное ядро алгоритма === | |
| − | + | Процесс решения задачи кластеризации при помощи минимального остовного дерева может быть разделен на три независимых этапа: | |
| − | + | # Приведение задачи наборе векторов данных к задаче о поиске минимального остовного дерева в графе (то есть нахождение расстояний между всеми парами данных векторов с использованием выбранной метрики, и запись этого расстояния как веса ребра, соединяющего эти вершины). Последовательная сложность: <math> O(n^2) </math> | |
| + | # Поиск минимального остовного дерева. Последовательная сложность: <math> O(E*log(E)) </math> | ||
| + | # Разбиение остовного дерева на кластеры. Последовательная сложность зависит от реализации. (от <math> O(1) </math> в случае SEMST до <math> O(n^2) </math> в случае ZEMST). | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
| − | + | В реализации алгоритма используются следующие структуры данных: | |
| + | #'''Узел'''(Node)(Вершина) - Представляется набором координат в выбранном пространстве. | ||
| + | #'''Ребро'''(Edge) - Представляется: двумя индексами вершин, которые соединяет; своим весом; флагом, который указывает, входит ли данное ребро в минимальное остовное дерево. | ||
| − | : | + | В дальнейшем описании алгоритма используются следующие макрооперации: |
| + | #Вычисление расстояния между вершинами. | ||
| + | #Добавление ребра в граф - Осуществляется путём выставление специального флага у соответствующего ребра в true; | ||
| + | #Удаление ребра из графа - Операция, обратная к добавление: соответствующий флаг устанавливается в false. В начальный момент времени у всех рёбер флаг установлен в false, что означает отсутствие рёбер в графе. | ||
| − | в | + | Следующие три используемые операции относятся к реализации алгоритма Крускала при помощи леса непересекающихся множеств. Подробнее об этой структуре и связанных операциях можно почитать, например, здесь <ref> Томас Кормен и др. Алгоритмы: построение и анализ = INTRODUCTION TO ALGORITHMS. — 2-е изд. — М.: «Вильямс», 2006. — С. 1296. — ISBN 0-07-013151-1. (Глава 21. Раздел 21.3: "Леса непересекающихся множеств) </ref>, в данном пункте они будут только перечислены: |
| + | #Создание дерева с одним узлом - Make-Set | ||
| + | #Объединение компонент связности - Union | ||
| + | #Поиск компонент связности - Find-Set | ||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
| − | + | <source lang="cpp"> | |
| + | int N; // Размер входных данных | ||
| + | int number_of_clusters; // Число кластеров | ||
| + | int size = N* (N-1) / 2; // Количество связей в полном графе | ||
| + | Node nodes[N]; // Вершины(точки), заданные начальными условиями | ||
| + | Edge edges[size]; // Все расстояния между вершинами | ||
| + | vector<Edge> MST; // Рёбра минимального остовного дерева | ||
| + | read_data(&nodes, file_name); // Считывание начальных данных | ||
| − | 1 | + | //----------------Часть 1--------------// |
| + | for(int i = 0; i < N-1; i++) { | ||
| + | for(int j = i+1; j < N; j++) { | ||
| + | double distance = getDistance(node[i], node[j]); // Вычисление расстояния между вершинами | ||
| + | edges[i*N + j] = Edge(i,j,distance); // Запись значения расстояния в массив | ||
| + | } | ||
| + | } | ||
| − | |||
| − | + | //----------------Часть 2--------------// | |
| + | sort(edges, size); // Сортировка рёбер графа в порядке возрастания | ||
| + | for(int i = 0; i < size; i++) { // Обходим список рёбер по возрастанию | ||
| + | if(find_set(edges[i].first) != find_set(edges[i].second)) { // Если вершины принадлежат разным подграфам - связываем их | ||
| + | Union(edges[i].first, edges[i].second); // Если вершины принадлежат разным подграфам - связываем их | ||
| + | MST.push_back(edges[i]); // Запись использованного ребра в минимальное остовное дерево | ||
| + | edges[i]->is_in_MST = true; // Запись использованного ребра в минимальное остовное дерево | ||
| + | } | ||
| + | } | ||
| − | |||
| − | + | //----------------Часть 3--------------// | |
| − | + | for(int i = MST.size() - 1; i > MST.size() - 1 - number_of_clusters; i--) { | |
| − | + | remove_edge(MST, i); // Удаление самых длинных рёбер из MST - формирование кластеров | |
| + | } | ||
| + | write_data(MST); // Запись выходных данных | ||
| + | </source> | ||
| − | + | Как описано в пункте 1.3, весь описываемый алгоритм (начиная от момента получения данных и заканчивая моментом, когда данные отданы на вывод), может быть условно разделён на три части. В то время как первая часть может быть достаточно эффективно решена "в лоб", вторая и третьи части, как указано в пункте 1.1, могут иметь альтернативы. В данном фрагменте кода подробно описано использование [[Алгоритм Крускала | Алгоритма Крускала]] для построения MST, и алгоритма SEMST для разбиения на кластеры. | |
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
| − | + | Как и в остальных пунктах, здесь сложность алгоритма состоит их суммы сложностей трёх его составных частей. Если оценивать сложность асимптотически в зависимости от входных данных, то: | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | <math> T(n) = O(n^2) + O(\frac{n (n-1)}{2})log(\frac{n (n-1)}{2})) + O(1) </math> | |
| − | |||
| − | |||
| − | + | На первый взгляд второе слагаемое может показаться слишком большим. Однако, верны и следующие равенства: | |
| − | === | + | <math> O(\frac{n (n-1)}{2} log(\frac{n (n-1)}{2}) = O(n^2 log(n^2)) = O(2 n^2 log(n)) = O(n^2 log(n)) </math> |
| − | + | Отсюда видно, что для входных данных этого алгоритма асимптотические сложности алгоритмов Крускала, Прима и Борувки совпадают. | |
| − | + | Таким образом, последовательная сложность алгоритма: <math> O(n^2 log(n)) </math> | |
| − | + | === Информационный граф === | |
| − | |||
| − | + | [[File:Graph_algo10.png|frame|center|alt=Алгоритм| Макро-описание параллельного алгоритма кластеризации]] | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
| + | Данная задача хорошо распараллеливается благодаря тому, что все "тяжелые" вычислительные задачи, начиная с момента считывания входных данных, можно распределить по процессам практически равномерно. Это достигается за счет следующей реализации: | ||
| + | #Мастер-процесс считывает входные данные и отправляет их остальным процессам с помощью MPI_Broadcast. | ||
| + | #Все пары точек разбиваются так, что каждому процессу достается одинаковое количество пар. | ||
| + | #Каждый процесс-рабочий строит MST на своём наборе пар точек. | ||
| + | #Процессы объединяют свои деревья по алгоритму, изображенному в информационном графе. То есть объединение осуществляется последовательно между группами размера <math>2^k</math>, где <math> k \in 1,2,..log_2p</math>, где p - число процессов. | ||
| − | + | Таким образом '''высота Ярусно-Параллельной формы''' (при наличии <math> p = 2^i </math> процессоров) будет <math> 1 + i + 1 + 1 </math>. | |
| − | + | В то же время, как будет показано в пункте 1.10, оптимальное количетво процессоров может быть оценено по формуле <math> p(N) \approx \sqrt {\frac{N-1}{2}} </math>. Тогда '''ширину Ярусно-Параллельной формы''' можно оценить как <math> O(\sqrt{N}) </math>. Учитывая приближенную формулу для числа процессоров, оценим высоту ЯП при помощи соотношения <math> 2^i \approx \sqrt {\frac{N-1}{2}} </math>. Тогда высоту ЯП можно оценить как <math> O(log(N)) </math>. | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | Говоря о праллельной реализации вообще и, в частности, обращаясь к информационному графу, видно, что в основе параллельного алгоритма лежит массовый параллелизм, и единственным проблематичным местом может оказаться самый начальный этап считывания данных мастер-процессом с последующей рассылкой их всем процессам-рабочим. Проблема в том, что объем данных, пересылаемых каждому процессу, практически нельзя уменьшить из-за того, что в расчете MST на каждом наборе рёбер потенциально могут участвовать все вершины графа, и поэтому для их обработки требуется знание их пространственных координат. | |
| − | + | Самой "тяжелой" операцией здесь является самый первый расчёт MST, поскольку он выполняется чуть больше чем за квадратичное время. Все последующие объединения деревьев занимают асимптотически одинаковое количество процессорного времени и не зависят друг от друга (вплоть до последнего объединения в MST). | |
=== Входные и выходные данные алгоритма === | === Входные и выходные данные алгоритма === | ||
| − | + | ==== Входные данные ==== | |
| − | + | В качестве входных данных программе дается множество векторов <math> [x_1,x_2,...,x_N]</math>, где <math> x_j = (a_{1j}, a_{2j}, ..., a_{nj})</math> размерности <math>n</math>. Для удобства изучения свойств алгоритма и отображения результатов можно брать массив векторов размерности 2. В этом случае наглядность будет обеспечена естественной возможностью отобразить вектора в плоскости. | |
| − | |||
| − | |||
| − | + | Кроме того может оказаться полезным особым образом подготовить данные перед анализом (удаление "шумов", выбросов и масштабирование данных по тем или иным координатам). | |
| + | ==== Выходные данные ==== | ||
| + | Выходными данными программы будет множество <math> M = ((x_1,x_2,..,x_k),(x_{k+1},...),...) </math>, состоящее из сгруппированных по кластерам векторов из входных данных. | ||
| − | + | В конкретной программе это может быть реализовано путём записи <math> N </math> строк вида <math>(x_{i}, Cluster-Number)</math> в выходной файл. | |
| − | |||
| − | |||
=== Свойства алгоритма === | === Свойства алгоритма === | ||
| − | |||
| − | + | Одним из трудных мест параллельной реализации является правильный выбор количества процессов. Для удобства анализа рассматривается число процессоров <math> p = 2^i </math>, где <math> i \in \mathbb {N} </math>. Ясно, что такое допущение не повлияет на качество анализа, так как другое число процессоров лишь приведет к незначительным накладным расходам, сложность которых асимптотически стремится к нулю. В таком случае можно также ввести допущение о размере входных данных: <math> N = 2^j </math>, где <math> j \in \mathbb {N}, j > i </math>. В этом случае MST-часть дерева графа алгоритма является идеально сбалансированным бинарным деревом. Тогда сложность можно оценить следующим образом: | |
| + | #На первом этапе каждый процесс производит расчет на <math> 2^{j - i} * (2^j - 1)</math> рёбрах, получая локальное MST(или MSF) на не более чем <math> N </math> вершинах. | ||
| + | #Далее следует <math> i = log_2(2^i) </math> шагов объединения не более чем по N рёбер в каждом. | ||
| − | + | Таким образом, задача эффективного запуска алгоритма состоит в том, чтобы найти оптимальный баланс между количеством расчетов на каждом узле на первом шаге и общим количеством объединений поддеревьев. Пусть дано <math> p = 2^i </math> вычислительных узлов и <math> N = 2^j </math>. В таком случае каждый узел получает считает <math> 2^{j - i} * (2^j - 1)</math> рёбер на первом шаге. Высота части ЯП-формы, отвечающей за MST, равна <math> i </math>. Тогда всего будет <math> i </math> уровней объединений MST с числом рёбер не более чем <math> 2N </math> для каждого вычислительного узла. | |
| − | + | ==== Оценка скорости работы параллельного алгоритма ==== | |
| − | + | Грубая оценка времени вычисления может помочь изучить качественные характеристики алгоритма. Учитывая рассуждения выше, приближенная качественная формула может быть записана таким образом: | |
| − | + | *<math> T'(N,p) = \frac{N(N - 1)}{2 p} * logN + 2N*log(N) * log_2p</math>. | |
| − | |||
| − | |||
| − | </math> | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | Объяснение формулы: первое слагаемое, исходя из оценки для сложности алгоритма Прима: <math> O(ElogV) </math>, - время(число операций) расчета на каждом процессоре на первом шаге. После первого шага (поскольку построено минимальное остовное дерево) каждый процесс передает не более чем <math> N </math> рёбер дальше. Второе слагаемое: каждый процесс после первого шага получает не более чем <math> 2N </math> рёбер (от двух процессоров при объединении), и всего таких шагов объединения не более чем <math> log_2p </math>, при этом предполагается, что сложностью объединения <math> 2N </math> рёбер на каждом шаге можно пренебречь (или допустить его пренебрежительно малым), поскольку предполагается, что все процессы после первого шага получают свою часть уже отсортированных рёбер(то есть наибольшее ребро процессора <math> i </math> меньше наименьшего ребра процессора <math> i+1 </math>), и объединение будет просто конкатенацией блоков памяти. | |
| − | + | С помощью этой формулы можно построить графики зависимости времени вычислений от количества процессоров и входных данных.Поскольку рост вычислительного времени при увеличении объёма входных данных достаточно очевиден, этот объём будет зафиксирован, и будет исследоваться поведение приближенного априорного вычислительного времени алгоритма в зависимости от числа процессоров. Проведём два вычислительных эксперимента: | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | 1. <math> N = 10^4</math> | |
| + | [[File:График1_.JPG|frame|center|alt=Алгоритм| По оси абсцисс: количество вычислительных узлов. По оси ординат: время вычислений. При фиксированном <math> N </math>]] | ||
| − | + | 2. <math> N = 1.5*10^4 </math> | |
| + | [[File:График2_.JPG|frame|center|alt=Алгоритм| По оси абсцисс: количество вычислительных узлов. По оси ординат: время вычислений. При фиксированном <math> N </math>]] | ||
| − | + | Из этих графиков становится видно качественное поведение времени вычислений при увеличении числа процессоров. Можно заметить, что для каждого набора входных данных есть число процессоров, минимизирующее время вычислений. Но в то же время видно, что при слишком большом количестве процессоров алгоритм начинает работать хуже, чем при меньшем. Это вполне понятно: Слишком большое число процессоров порождает большое число обменов между ними, объединений деревьев, сортировок, и т.д. Таким образом, именно эта часть <math> 2N*log(N) * log_2p </math> оценки времени вычислений отвечает за ухудшение эффективности алгоритма при числе процессоров, которое больше оптимального. Кроме того, видно, что при увеличении объёма входных данных оптимальное количество процессоров также увеличивается: количество обменов и слияний деревьев компенсируется уменьшением вычислительной нагрузки на каждый отдельный процессор. | |
| − | + | Благодаря такой оценке становится возможным оценить примерное оптимальное число процессоров <math> p </math> в зависимости от объёма входных данных <math> N </math>. Продифференцировав приближенную формулу <math> T'(N,p) = \frac{N(N - 1)}{2 p} * logN + 2N*log(N) * log_2p </math> по <math> p </math> и приравняв производную нулю, получим следующее выражение для оптимального числа процессоров: | |
| − | < | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | </ | ||
| − | |||
| − | + | <math> p(N) \approx \frac{N-1}{4} </math> | |
| − | ==== | + | ==== Вычислительная мощность параллельного алгоритма ==== |
| − | + | Согласно определению, [[глоссарий#Вычислительная мощность|''вычислительная мощность'']] алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. В то же время, учитывая вышесказанное, видно, что число операций в данном алгоритме зависит не только от объема входных данных, но также от числа процессоров. Таким образом, имеет смысл оценить вычислительную мощность как в случае оптимального числа процессоров, так и в общем случае числа процессоров, равного <math> p </math>. | |
| − | + | В общем случае, исходя из предыдущего пункта, число операций будет равно | |
| − | + | <math> \frac{N(N - 1)}{2} * logN + \sum_{i = 0}^{log_2p}(2N*log(N) * \frac{p}{2^i}) = N logN (\frac{N-1}{2} + 2 * p * \sum_{i = 0}^{log_2p}( \frac{1}{2^i})) = </math> {переходя к пределу в сумме} <math> \approx N logN (\frac{N-1}{2} + 2 * p) </math> | |
| − | + | Таким образом, вычислительная мощность в общем случае: <math> \approx \frac{N logN (\frac{N-1}{2} + 2 * p)}{2N} = \frac{logN *(\frac{N-1}{2} + 2 * p)}{2} </math> | |
| − | + | В случае оптимального числа вычислительных узлов: <math> \approx \frac{N logN (\frac{N-1}{2} + \frac{N-1}{2})}{2N} = \frac{logN * (N-1)}{2}</math> | |
| − | + | == Программная реализация алгоритма == | |
| − | |||
| − | + | === Особенности реализации последовательного алгоритма === | |
| − | + | === Локальность данных и вычислений === | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | === | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
=== Возможные способы и особенности параллельной реализации алгоритма === | === Возможные способы и особенности параллельной реализации алгоритма === | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
=== Масштабируемость алгоритма и его реализации === | === Масштабируемость алгоритма и его реализации === | ||
| − | + | В настоящей статье приводится исследование самостоятельной реализации алгоритма: [http://pastebin.com/gH09k6id main.cpp]. Запуск производился на системе IBM Blue Gene/P. | |
| + | Для генерации данных использовался скрипт на питоне [http://pastebin.com/5YNxqqpA testGenerator]. Для отображения кластеров использовалась библиотека matplotlib и скрипт [http://pastebin.com/mTDZ6bqp graphDraw]. | ||
| − | + | IBM Blue Gene/P — массивно-параллельная вычислительная система, которая состоит из двух стоек, включающих 2048 четырехъядерных вычислительных узлов, с пиковой производительностью 27,9 терафлопс. | |
| − | |||
| − | + | Характеристики системы: | |
| + | * две стойки с вычислительными узлами и узлами ввода-вывода | ||
| + | * 1024 четырехъядерных вычислительных узла в каждой из стоек | ||
| + | * 16 узлов ввода-вывода в стойке (в текущей конфигурации активны 8, т.е. одна I/O-карта на 128 вычислительных узлов) | ||
| + | * программирование с использованием MPI, OpenMP/pthreads, POSIX I/O | ||
| + | * высокая энергоэффективность: ∼ 372 MFlops/W | ||
| − | + | Вычислительный узел включает в себя четырехъядерный процессор, 2 ГБ общей памяти и сетевые интерфейсы. | |
| − | |||
| − | + | Вычислительной узел: | |
| + | * четыре микропроцессорных ядра PowerPC 450 (4-way SMP) | ||
| + | * пиковая производительность: 4 ядра x 3,4 GFlop/sec на ядро = 13,6 GFlop/sec | ||
| + | * пропускная способность памяти: 13,6 GB/sec | ||
| + | * 2 ГБ общей памяти | ||
| + | * 2 x 4 МБ кэш-памяти 2-го уровня (в документации по BG/P носит название L3) | ||
| + | * легковесное ядро (compute node kernel, CNK), представляющее собой Linux-подобную операционную систему, поддерживающую значительное подмножество Linux-совместимых системных вызовов | ||
| + | * асинхронные операции межпроцессорных обменов (выполняются параллельно с вычислениями) | ||
| + | * операции ввода-вывода перенаправляются I/O-картам через сеть коллективных операций | ||
| − | |||
| − | |||
| − | + | Поскольку оптимальное число процессоров достаточно велико и линейно зависит от объема(размера) входных данных, задача хорошо масштабируется, и на графике ниже можно увидеть почти линейный рост производительности при фиксированном количестве входных данных и увеличении числа процессоров в том случае, когда число процессоров близко к оптимальному. В тоже время из графика видно, что когда входных данных слишком мало, а число процессоров возрастает, то происходит либо небольшой рост производительности, либо роста производительности нет. | |
| − | [[ | + | [[File:Performance3.JPG|frame|center|alt=Алгоритм| Производительность алгоритма]] |
| − | |||
| − | + | Говоря о вычислительном времени и проводя тесты на реальных данных, можно заметить концептуальную правильность сделанных качественных выводов о поведении времени вычислений для различного набора входных данных и числа процессоров. На графиках ниже видно, что время вычисления задачи уменьшается до определенного момента, когда количество процессоров, нужное для задачи, становится слишком большим, после чего начинает увеличиваться, когда число процессоров растет. | |
| − | |||
| − | |||
| − | |||
| − | [ | + | [[File:Timeline1.JPG|frame|center|alt=Алгоритм| Вычислительное время алгоритма для <math>N = 17000</math>]] |
| + | [[File:Снимок2.JPG|frame|center|alt=Алгоритм| Вычислительное время алгоритма для <math>N = 3000</math>]] | ||
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
=== Выводы для классов архитектур === | === Выводы для классов архитектур === | ||
| − | |||
| − | |||
| − | |||
| − | |||
=== Существующие реализации алгоритма === | === Существующие реализации алгоритма === | ||
| + | Как таковые, алгоритмы самого процесса кластеризации в популярных библиотеках (таких, как Scipy, Numpy или Boost) не реализованы (в отличие от других алгоритмов кластеризации). Однако, можно найти пользовательские любительские реализации (некоторые из которых приведены в списке ниже). Кроме того, существует масса реализаций основного коспонента данного алгоритма - вычисления MST - на многих языках программирования в разных библиотеках. Ссылки на некоторые из этих реализаций(перечислить все было бы невозможно) также приведены в списке ниже. | ||
| − | + | *[https://github.com/jakevdp/mst_clustering Minimum Spanning Tree Clustering in Python Author: Jake VanderPlas] | |
| − | + | *[https://github.com/amueller/information-theoretic-mst Information Theoretic Clustering using Minimum Spanning Trees] | |
| − | + | *[http://www.boost.org/doc/libs/1_51_0/libs/graph/doc/prim_minimum_spanning_tree.html C++ Boost library prim_minimum_spanning_tree] | |
| − | + | *[http://www.boost.org/doc/libs/1_35_0/libs/graph/doc/kruskal_min_spanning_tree.html C++ Boost library kruskal_minimum_spanning_tree] | |
| − | + | *[http://jgrapht.org/javadoc/org/jgrapht/alg/PrimMinimumSpanningTree.html Java library JGraphT PrimMinimumSpanningTree] | |
| − | + | *[http://jgrapht.org/javadoc/org/jgrapht/alg/KruskalMinimumSpanningTree.html Java library JGraphT KruskalMinimumSpanningTree] | |
| − | + | *[http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.sparse.csgraph.minimum_spanning_tree.html SciPy library MST] | |
== Литература == | == Литература == | ||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
Текущая версия на 14:43, 23 ноября 2016
 | Эта работа прошла предварительную проверку Дата последней правки страницы: 23.11.2016 Данная работа соответствует формальным критериям. Проверено Zhum. |
| Алгоритм кластеризации, основанный на минимальном остовном дереве | |
| Последовательный алгоритм | |
| Последовательная сложность | [math] O(N^2 log(N)) [/math] |
| Объём входных данных | [math] N [/math] |
| Объём выходных данных | [math] N [/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math] O(log(N)) [/math] |
| Ширина ярусно-параллельной формы | [math] O(N) [/math] |
Основные авторы описания: Стандрик Антон Сергеевич
Содержание
- 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 Общее описание алгоритма
Данный алгоритм предназначен для решения задач кластеризации. Говоря в общем, кластеризация — многомерная статистическая процедура, выполняющая сбор данных, содержащих информацию о выборке объектов, и затем упорядочивающая объекты в сравнительно однородные группы [1].
Впервые идея использования минимального остовного дерева для кластеризации данных была преждожена C. T. Zah в 1971 году[2]. Этот метод широко используется исследователями при анализе изображений [3][4], а также при анализе биологических данных [5].
В общем описании работы алгоритма можно выделить следующие шаги:
- Разбиение входных векторов на пары и вычисление расстояний в них. Для этого небходим выбор метрики пространства, в котором осуществляется кластерный анализ. В качестве примеров метрик можно привести:
- Метрика Евклида: [math]d_{ij} = \sqrt{\sum^n_{k=1}(x_{ik} - x_{jk})^2}[/math]
- Метрика city-block (манхэттенская метрика): [math]d_{ij} = \sum^n_{k=1}|x_{ik} - x_{jk}| [/math]
- Метрика Брея-Картиса: [math]d_{ij} = \sum^n_{k=1}\frac{|x_{ik} - x_{jk}|}{|x_{ik} + x_{jk}| }[/math]
- Построение минимального остовного дерева по заданному множеству. Существует три основных алгоритма построения минимального остовного дерева:
- Разделение полученного дерева на несвязное множество поддеревьев в соответствии с некоторой стратегией. Можно выделить три популярные стратегии кластеризации:
- SEMST[6] (standard Euclidean minimum spanning tree clustering algorithm) - предполагается, что нам заранее известно число кластеров. В таком случае, если число кластеров равно [math]k[/math], то из полученного в первом пункте остовного дерева удаляются [math]k-1[/math] связей с наибольшим весом.
- CEMST[7](The maximum cost minimum spanning tree) - Работает так же как и SEMST, однако вместо веса вершины используется её стоимость, которая вычисляется следующим образом: Пусть ребро [math]e = (x_{i}, x_{j})[/math], тогда [math]cost(e) = w(e) * degree(x_{i}) * degree(x_{j})[/math]. Далее, зная число кластеров [math]k[/math] и стоимости рёбер, удаляются [math]k-1[/math] рёбер с максимальной стоимостью.
- ZEMST[2] (The Zahn’s minimum spanning tree) - от двух предыдущих эта стратегия отличается тем, что здесь число кластеров не известно заранее, а вычисляется в процессе кластеризации. Его суть в том, что удаляются те рёбра, веса которых значительно превосходят средний вес ребра в соседних поддеревьях.
1.1.1 Достоинства
Основным достоинством алгоритма является его способность выделять кластеры произвольной формы. На следующем примере можно видеть, как алгоритм даже в случае SEMST алгоритма корректно разбивает вложенные окружности на три кластера:
Слева направо: Рисунок 1: Начальные данные. Рисунок 2: построенное минимальное остовное дерево. Рисунок 3: разбиение на кластеры путём удаления самых длинных связей.
1.1.2 Недостатки
Недостатком является трудность определения четких критериев, при которых осуществляется разбиение данных на разные кластеры.
Кроме того, существуют наборы данных, на которых популярные алгоритмы (SEMST, CEMST, ZEMST) не дают правильного результата, и для которых требуется применение более сложных алгоритмов, к числу которых можно отнести, например, алгоритм HEMST [8] с временем работы [math]O(n^2 * log(n))[/math].
В качестве иллюстрации можно привести два похожих визуально, но различных с точки зрения разбиения на кластеры с помощью описываемого алгоритма, примера. (В приведенных примерах первый рисунок обозначает начальный набор данных. Второй - построенное на основе данных минимальное остовное дерево. Третий рисунок - разбиение на кластеры.):
Правильное разбиение на два кластера.
Неверное разбиение на один кластер.
Отличие второго примера - в наличии плотной цепочки связей между двумя кластерами, которая и приводит к ошибке наивного алгоритма удаления самых длинных связей, и таким образом здесь потребуется более сложный алгоритм с более значительным временем работы для правильного нахождения кластеров.
1.2 Математическое описание алгоритма
Исходные данные: набор векторов размерности [math]n[/math] : [math](a_1, a_2, ..., a_n)[/math], где [math] a_j \in \mathbb {R} \quad \forall j \in [1,n] [/math].
Вычисляемые данные: Сгруппированные по кластерам данные.
1.2.1 Основные обозначения
- [math]E[/math] - множество рёбер.
- [math]V[/math] - множество вершин.
- [math]G = (V,E) [/math] - неориентированный граф, заданный множеством вершин [math]V[/math] и множеством рёбер [math]E[/math]
- [math]MST(G)[/math] (Minimum Spanning Tree) - минимальное остовное дерево графа [math][/math]. Остовное дерево — ациклический связный подграф данного связного неориентированного графа, в который входят все его вершины. Минимальное остовное дерево в связанном взвешенном неориентированном графе — это остовное дерево этого графа, имеющее минимальный возможный вес, где под весом дерева понимается сумма весов входящих в него рёбер.
- [math]MSF(G)[/math] (Minimum Spanning Forest) - набор минимальных остовных деревьев для всех компонент связности графа [math] G [/math]
В описании алгоритмов этой статьи часто используется термин "построение MST", подразумевая, что на самом деле в результате построения может получиться как одна компонента связности (MST), так и множество компонент связности (MSF).
1.2.2 Алгоритмы кластеризации
SEMST(standard Euclidean minimum spanning tree) - простейший алгоритм, основанный на заранее известном числе кластеров [math] k [/math]. Путём удаления [math] k - 1 [/math] связей с наибольшим весом предполагается разбиение остовного дерева на нужное количество кластеров. Иными словами минимальное остовное дерево разбивается на [math] k [/math] поддеревьев: [math] G = \{G_1,G_2,...,G_k\} [/math] путём минимизации суммы всех весов: [math] W = \sum ^k_{l=1} \sum _{x_i,x_j \in V_l} w_{x_i,x_j} -\gt min [/math], где [math] V_l [/math] - множество вершин в поддереве [math] G_l [/math]. В конце такой процедуры каждое множество [math] V_l [/math] становится кластером [math] C_l [/math]. Поскольку (в случае использования алгоритма Крускала для построения остовного дерева) рёбра перед началом кластеризации уже отсортированы, для удаления наибольших [math] k [/math] рёбер требуется [math] O(k) [/math] операций, и для распределения получившихся вершин по кластерам требуется время [math] O(n) [/math].
CEMST(maximum cost Euclidean minimum spanning tree) - Работает в точности, как и SEMST, но здесь вместо весов рёбер используется т.н. "стоимость", которая высчитывается следующим образом: если ребро [math] (x_i, x_j) [/math] соединяет вершины [math] x_i [/math] и [math] x_j [/math], то его "стоимость" считается как: [math] Cost_{x_i, x_j} = w_{x_i,x_j} * Deg(x_i) * Deg(x_j) [/math], где [math] Deg(x_i) [/math] - степень вершины [math] x_i [/math] (т.е. количество рёбер графа G, инцидентных вершине [math] x_i [/math]). В остальном алгоритм действует точно так же, как и SEMST, удаляя [math] k-1 [/math] рёбер с наибольшей "стоимостью". Сложность [math] O(k) [/math] . Алгоритм направлен на удаление рёбер, соединяющих высоко нагруженные (с точки зрения количества инцидентных рёбер) вершины.
ZEMST (The Zahn’s maximum spanning tree (Zahn, 1971)) [9] - отличается от двух предыдущих алгоритмов в первую очередь тем, что здесь количество кластеров заранее не известно, но определяется исходя из данных. Данный алгоритм не является жадным, в отличие от первых двух, так как в нём рассматривается не только само ребро и его вес, а также и поддеревья, которые соединены этим ребром. То есть для ребра [math] (x_i, x_j) [/math] рассматриваются поддеревья [math] N_i = (V_{N_i}, E_{N_i}), N_j = (V_{N_j}, E_{N_j}) [/math], где [math] N_i [/math] - поддерево глубины(высоты) [math] d [/math](величина d определяется эмпирически), состоящее из всех достижимых из вершины [math] x_j [/math] вершин за [math] d [/math] шагов через все достижимые из неё рёбра (исключая само ребро [math] (x_i, x_j) [/math]). Введём обозначения:
- [math]\overline{w_{N_i}} = \frac{\sum_{(x_r, x_s) \in E_{N_j}}w_{x_r,x_s}}{|E_{N_j}|} [/math] - средний вес ребра в поддереве [math] N_i [/math]
- [math]\sigma_{N_i} = \sqrt{\frac{\sum_{(x_r, x_s) \in E_{N_j}}(w_{x_r,x_s} - \overline{w_{N_i}})^2}{|E_{N_j}|}} [/math] - стандартное отклонение веса ребра в поддереве [math] N_i [/math]
Ребро [math] (x_i, x_j) [/math] удаляется из списка в случае, если выполняется хотя бы одно из двух условий:
- [math] w_{x_i,x_j} \gt \overline{w_{N_i}} - \sigma_{N_i} [/math]
- [math] w_{x_i,x_j} \gt \overline{w_{N_j}} - \sigma_{N_j} [/math]
- [math] w_{x_i,x_j} \gt \frac{\sum_{(x_r, x_s) \in E_i \cup E_j}w_{x_r,x_s}}{|E_i \cup E_j|} [/math]
Сложность этого алгоритма составляет в худшем случае [math] O(N^2) [/math] (поскольку может оказаться нужным перебрать все рёбра), а в лучшем [math] O(Nd) [/math].
1.3 Вычислительное ядро алгоритма
Процесс решения задачи кластеризации при помощи минимального остовного дерева может быть разделен на три независимых этапа:
- Приведение задачи наборе векторов данных к задаче о поиске минимального остовного дерева в графе (то есть нахождение расстояний между всеми парами данных векторов с использованием выбранной метрики, и запись этого расстояния как веса ребра, соединяющего эти вершины). Последовательная сложность: [math] O(n^2) [/math]
- Поиск минимального остовного дерева. Последовательная сложность: [math] O(E*log(E)) [/math]
- Разбиение остовного дерева на кластеры. Последовательная сложность зависит от реализации. (от [math] O(1) [/math] в случае SEMST до [math] O(n^2) [/math] в случае ZEMST).
1.4 Макроструктура алгоритма
В реализации алгоритма используются следующие структуры данных:
- Узел(Node)(Вершина) - Представляется набором координат в выбранном пространстве.
- Ребро(Edge) - Представляется: двумя индексами вершин, которые соединяет; своим весом; флагом, который указывает, входит ли данное ребро в минимальное остовное дерево.
В дальнейшем описании алгоритма используются следующие макрооперации:
- Вычисление расстояния между вершинами.
- Добавление ребра в граф - Осуществляется путём выставление специального флага у соответствующего ребра в true;
- Удаление ребра из графа - Операция, обратная к добавление: соответствующий флаг устанавливается в false. В начальный момент времени у всех рёбер флаг установлен в false, что означает отсутствие рёбер в графе.
Следующие три используемые операции относятся к реализации алгоритма Крускала при помощи леса непересекающихся множеств. Подробнее об этой структуре и связанных операциях можно почитать, например, здесь [10], в данном пункте они будут только перечислены:
- Создание дерева с одним узлом - Make-Set
- Объединение компонент связности - Union
- Поиск компонент связности - Find-Set
1.5 Схема реализации последовательного алгоритма
int N; // Размер входных данных
int number_of_clusters; // Число кластеров
int size = N* (N-1) / 2; // Количество связей в полном графе
Node nodes[N]; // Вершины(точки), заданные начальными условиями
Edge edges[size]; // Все расстояния между вершинами
vector<Edge> MST; // Рёбра минимального остовного дерева
read_data(&nodes, file_name); // Считывание начальных данных
//----------------Часть 1--------------//
for(int i = 0; i < N-1; i++) {
for(int j = i+1; j < N; j++) {
double distance = getDistance(node[i], node[j]); // Вычисление расстояния между вершинами
edges[i*N + j] = Edge(i,j,distance); // Запись значения расстояния в массив
}
}
//----------------Часть 2--------------//
sort(edges, size); // Сортировка рёбер графа в порядке возрастания
for(int i = 0; i < size; i++) { // Обходим список рёбер по возрастанию
if(find_set(edges[i].first) != find_set(edges[i].second)) { // Если вершины принадлежат разным подграфам - связываем их
Union(edges[i].first, edges[i].second); // Если вершины принадлежат разным подграфам - связываем их
MST.push_back(edges[i]); // Запись использованного ребра в минимальное остовное дерево
edges[i]->is_in_MST = true; // Запись использованного ребра в минимальное остовное дерево
}
}
//----------------Часть 3--------------//
for(int i = MST.size() - 1; i > MST.size() - 1 - number_of_clusters; i--) {
remove_edge(MST, i); // Удаление самых длинных рёбер из MST - формирование кластеров
}
write_data(MST); // Запись выходных данных
Как описано в пункте 1.3, весь описываемый алгоритм (начиная от момента получения данных и заканчивая моментом, когда данные отданы на вывод), может быть условно разделён на три части. В то время как первая часть может быть достаточно эффективно решена "в лоб", вторая и третьи части, как указано в пункте 1.1, могут иметь альтернативы. В данном фрагменте кода подробно описано использование Алгоритма Крускала для построения MST, и алгоритма SEMST для разбиения на кластеры.
1.6 Последовательная сложность алгоритма
Как и в остальных пунктах, здесь сложность алгоритма состоит их суммы сложностей трёх его составных частей. Если оценивать сложность асимптотически в зависимости от входных данных, то:
[math] T(n) = O(n^2) + O(\frac{n (n-1)}{2})log(\frac{n (n-1)}{2})) + O(1) [/math]
На первый взгляд второе слагаемое может показаться слишком большим. Однако, верны и следующие равенства:
[math] O(\frac{n (n-1)}{2} log(\frac{n (n-1)}{2}) = O(n^2 log(n^2)) = O(2 n^2 log(n)) = O(n^2 log(n)) [/math]
Отсюда видно, что для входных данных этого алгоритма асимптотические сложности алгоритмов Крускала, Прима и Борувки совпадают.
Таким образом, последовательная сложность алгоритма: [math] O(n^2 log(n)) [/math]
1.7 Информационный граф
Макро-описание параллельного алгоритма кластеризации
1.8 Ресурс параллелизма алгоритма
Данная задача хорошо распараллеливается благодаря тому, что все "тяжелые" вычислительные задачи, начиная с момента считывания входных данных, можно распределить по процессам практически равномерно. Это достигается за счет следующей реализации:
- Мастер-процесс считывает входные данные и отправляет их остальным процессам с помощью MPI_Broadcast.
- Все пары точек разбиваются так, что каждому процессу достается одинаковое количество пар.
- Каждый процесс-рабочий строит MST на своём наборе пар точек.
- Процессы объединяют свои деревья по алгоритму, изображенному в информационном графе. То есть объединение осуществляется последовательно между группами размера [math]2^k[/math], где [math] k \in 1,2,..log_2p[/math], где p - число процессов.
Таким образом высота Ярусно-Параллельной формы (при наличии [math] p = 2^i [/math] процессоров) будет [math] 1 + i + 1 + 1 [/math]. В то же время, как будет показано в пункте 1.10, оптимальное количетво процессоров может быть оценено по формуле [math] p(N) \approx \sqrt {\frac{N-1}{2}} [/math]. Тогда ширину Ярусно-Параллельной формы можно оценить как [math] O(\sqrt{N}) [/math]. Учитывая приближенную формулу для числа процессоров, оценим высоту ЯП при помощи соотношения [math] 2^i \approx \sqrt {\frac{N-1}{2}} [/math]. Тогда высоту ЯП можно оценить как [math] O(log(N)) [/math].
Говоря о праллельной реализации вообще и, в частности, обращаясь к информационному графу, видно, что в основе параллельного алгоритма лежит массовый параллелизм, и единственным проблематичным местом может оказаться самый начальный этап считывания данных мастер-процессом с последующей рассылкой их всем процессам-рабочим. Проблема в том, что объем данных, пересылаемых каждому процессу, практически нельзя уменьшить из-за того, что в расчете MST на каждом наборе рёбер потенциально могут участвовать все вершины графа, и поэтому для их обработки требуется знание их пространственных координат.
Самой "тяжелой" операцией здесь является самый первый расчёт MST, поскольку он выполняется чуть больше чем за квадратичное время. Все последующие объединения деревьев занимают асимптотически одинаковое количество процессорного времени и не зависят друг от друга (вплоть до последнего объединения в MST).
1.9 Входные и выходные данные алгоритма
1.9.1 Входные данные
В качестве входных данных программе дается множество векторов [math] [x_1,x_2,...,x_N][/math], где [math] x_j = (a_{1j}, a_{2j}, ..., a_{nj})[/math] размерности [math]n[/math]. Для удобства изучения свойств алгоритма и отображения результатов можно брать массив векторов размерности 2. В этом случае наглядность будет обеспечена естественной возможностью отобразить вектора в плоскости.
Кроме того может оказаться полезным особым образом подготовить данные перед анализом (удаление "шумов", выбросов и масштабирование данных по тем или иным координатам).
1.9.2 Выходные данные
Выходными данными программы будет множество [math] M = ((x_1,x_2,..,x_k),(x_{k+1},...),...) [/math], состоящее из сгруппированных по кластерам векторов из входных данных.
В конкретной программе это может быть реализовано путём записи [math] N [/math] строк вида [math](x_{i}, Cluster-Number)[/math] в выходной файл.
1.10 Свойства алгоритма
Одним из трудных мест параллельной реализации является правильный выбор количества процессов. Для удобства анализа рассматривается число процессоров [math] p = 2^i [/math], где [math] i \in \mathbb {N} [/math]. Ясно, что такое допущение не повлияет на качество анализа, так как другое число процессоров лишь приведет к незначительным накладным расходам, сложность которых асимптотически стремится к нулю. В таком случае можно также ввести допущение о размере входных данных: [math] N = 2^j [/math], где [math] j \in \mathbb {N}, j \gt i [/math]. В этом случае MST-часть дерева графа алгоритма является идеально сбалансированным бинарным деревом. Тогда сложность можно оценить следующим образом:
- На первом этапе каждый процесс производит расчет на [math] 2^{j - i} * (2^j - 1)[/math] рёбрах, получая локальное MST(или MSF) на не более чем [math] N [/math] вершинах.
- Далее следует [math] i = log_2(2^i) [/math] шагов объединения не более чем по N рёбер в каждом.
Таким образом, задача эффективного запуска алгоритма состоит в том, чтобы найти оптимальный баланс между количеством расчетов на каждом узле на первом шаге и общим количеством объединений поддеревьев. Пусть дано [math] p = 2^i [/math] вычислительных узлов и [math] N = 2^j [/math]. В таком случае каждый узел получает считает [math] 2^{j - i} * (2^j - 1)[/math] рёбер на первом шаге. Высота части ЯП-формы, отвечающей за MST, равна [math] i [/math]. Тогда всего будет [math] i [/math] уровней объединений MST с числом рёбер не более чем [math] 2N [/math] для каждого вычислительного узла.
1.10.1 Оценка скорости работы параллельного алгоритма
Грубая оценка времени вычисления может помочь изучить качественные характеристики алгоритма. Учитывая рассуждения выше, приближенная качественная формула может быть записана таким образом:
- [math] T'(N,p) = \frac{N(N - 1)}{2 p} * logN + 2N*log(N) * log_2p[/math].
Объяснение формулы: первое слагаемое, исходя из оценки для сложности алгоритма Прима: [math] O(ElogV) [/math], - время(число операций) расчета на каждом процессоре на первом шаге. После первого шага (поскольку построено минимальное остовное дерево) каждый процесс передает не более чем [math] N [/math] рёбер дальше. Второе слагаемое: каждый процесс после первого шага получает не более чем [math] 2N [/math] рёбер (от двух процессоров при объединении), и всего таких шагов объединения не более чем [math] log_2p [/math], при этом предполагается, что сложностью объединения [math] 2N [/math] рёбер на каждом шаге можно пренебречь (или допустить его пренебрежительно малым), поскольку предполагается, что все процессы после первого шага получают свою часть уже отсортированных рёбер(то есть наибольшее ребро процессора [math] i [/math] меньше наименьшего ребра процессора [math] i+1 [/math]), и объединение будет просто конкатенацией блоков памяти.
С помощью этой формулы можно построить графики зависимости времени вычислений от количества процессоров и входных данных.Поскольку рост вычислительного времени при увеличении объёма входных данных достаточно очевиден, этот объём будет зафиксирован, и будет исследоваться поведение приближенного априорного вычислительного времени алгоритма в зависимости от числа процессоров. Проведём два вычислительных эксперимента:
1. [math] N = 10^4[/math]
По оси абсцисс: количество вычислительных узлов. По оси ординат: время вычислений. При фиксированном [math] N [/math]
2. [math] N = 1.5*10^4 [/math]
По оси абсцисс: количество вычислительных узлов. По оси ординат: время вычислений. При фиксированном [math] N [/math]
Из этих графиков становится видно качественное поведение времени вычислений при увеличении числа процессоров. Можно заметить, что для каждого набора входных данных есть число процессоров, минимизирующее время вычислений. Но в то же время видно, что при слишком большом количестве процессоров алгоритм начинает работать хуже, чем при меньшем. Это вполне понятно: Слишком большое число процессоров порождает большое число обменов между ними, объединений деревьев, сортировок, и т.д. Таким образом, именно эта часть [math] 2N*log(N) * log_2p [/math] оценки времени вычислений отвечает за ухудшение эффективности алгоритма при числе процессоров, которое больше оптимального. Кроме того, видно, что при увеличении объёма входных данных оптимальное количество процессоров также увеличивается: количество обменов и слияний деревьев компенсируется уменьшением вычислительной нагрузки на каждый отдельный процессор.
Благодаря такой оценке становится возможным оценить примерное оптимальное число процессоров [math] p [/math] в зависимости от объёма входных данных [math] N [/math]. Продифференцировав приближенную формулу [math] T'(N,p) = \frac{N(N - 1)}{2 p} * logN + 2N*log(N) * log_2p [/math] по [math] p [/math] и приравняв производную нулю, получим следующее выражение для оптимального числа процессоров:
[math] p(N) \approx \frac{N-1}{4} [/math]
1.10.2 Вычислительная мощность параллельного алгоритма
Согласно определению, вычислительная мощность алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. В то же время, учитывая вышесказанное, видно, что число операций в данном алгоритме зависит не только от объема входных данных, но также от числа процессоров. Таким образом, имеет смысл оценить вычислительную мощность как в случае оптимального числа процессоров, так и в общем случае числа процессоров, равного [math] p [/math].
В общем случае, исходя из предыдущего пункта, число операций будет равно
[math] \frac{N(N - 1)}{2} * logN + \sum_{i = 0}^{log_2p}(2N*log(N) * \frac{p}{2^i}) = N logN (\frac{N-1}{2} + 2 * p * \sum_{i = 0}^{log_2p}( \frac{1}{2^i})) = [/math] {переходя к пределу в сумме} [math] \approx N logN (\frac{N-1}{2} + 2 * p) [/math]
Таким образом, вычислительная мощность в общем случае: [math] \approx \frac{N logN (\frac{N-1}{2} + 2 * p)}{2N} = \frac{logN *(\frac{N-1}{2} + 2 * p)}{2} [/math]
В случае оптимального числа вычислительных узлов: [math] \approx \frac{N logN (\frac{N-1}{2} + \frac{N-1}{2})}{2N} = \frac{logN * (N-1)}{2}[/math]
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
В настоящей статье приводится исследование самостоятельной реализации алгоритма: main.cpp. Запуск производился на системе IBM Blue Gene/P. Для генерации данных использовался скрипт на питоне testGenerator. Для отображения кластеров использовалась библиотека matplotlib и скрипт graphDraw.
IBM Blue Gene/P — массивно-параллельная вычислительная система, которая состоит из двух стоек, включающих 2048 четырехъядерных вычислительных узлов, с пиковой производительностью 27,9 терафлопс.
Характеристики системы:
- две стойки с вычислительными узлами и узлами ввода-вывода
- 1024 четырехъядерных вычислительных узла в каждой из стоек
- 16 узлов ввода-вывода в стойке (в текущей конфигурации активны 8, т.е. одна I/O-карта на 128 вычислительных узлов)
- программирование с использованием MPI, OpenMP/pthreads, POSIX I/O
- высокая энергоэффективность: ∼ 372 MFlops/W
Вычислительный узел включает в себя четырехъядерный процессор, 2 ГБ общей памяти и сетевые интерфейсы.
Вычислительной узел:
- четыре микропроцессорных ядра PowerPC 450 (4-way SMP)
- пиковая производительность: 4 ядра x 3,4 GFlop/sec на ядро = 13,6 GFlop/sec
- пропускная способность памяти: 13,6 GB/sec
- 2 ГБ общей памяти
- 2 x 4 МБ кэш-памяти 2-го уровня (в документации по BG/P носит название L3)
- легковесное ядро (compute node kernel, CNK), представляющее собой Linux-подобную операционную систему, поддерживающую значительное подмножество Linux-совместимых системных вызовов
- асинхронные операции межпроцессорных обменов (выполняются параллельно с вычислениями)
- операции ввода-вывода перенаправляются I/O-картам через сеть коллективных операций
Поскольку оптимальное число процессоров достаточно велико и линейно зависит от объема(размера) входных данных, задача хорошо масштабируется, и на графике ниже можно увидеть почти линейный рост производительности при фиксированном количестве входных данных и увеличении числа процессоров в том случае, когда число процессоров близко к оптимальному. В тоже время из графика видно, что когда входных данных слишком мало, а число процессоров возрастает, то происходит либо небольшой рост производительности, либо роста производительности нет.
Производительность алгоритма
Говоря о вычислительном времени и проводя тесты на реальных данных, можно заметить концептуальную правильность сделанных качественных выводов о поведении времени вычислений для различного набора входных данных и числа процессоров. На графиках ниже видно, что время вычисления задачи уменьшается до определенного момента, когда количество процессоров, нужное для задачи, становится слишком большим, после чего начинает увеличиваться, когда число процессоров растет.
Вычислительное время алгоритма для [math]N = 17000[/math]
Вычислительное время алгоритма для [math]N = 3000[/math]
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Как таковые, алгоритмы самого процесса кластеризации в популярных библиотеках (таких, как Scipy, Numpy или Boost) не реализованы (в отличие от других алгоритмов кластеризации). Однако, можно найти пользовательские любительские реализации (некоторые из которых приведены в списке ниже). Кроме того, существует масса реализаций основного коспонента данного алгоритма - вычисления MST - на многих языках программирования в разных библиотеках. Ссылки на некоторые из этих реализаций(перечислить все было бы невозможно) также приведены в списке ниже.
- Minimum Spanning Tree Clustering in Python Author: Jake VanderPlas
- Information Theoretic Clustering using Minimum Spanning Trees
- C++ Boost library prim_minimum_spanning_tree
- C++ Boost library kruskal_minimum_spanning_tree
- Java library JGraphT PrimMinimumSpanningTree
- Java library JGraphT KruskalMinimumSpanningTree
- SciPy library MST
3 Литература
- ↑ Классификация и кластер. Под ред. Дж. Вэн Райзина. М.: Мир, 1980. 390 с.
- ↑ 2,0 2,1 C. T. Zahn, “graph-Theriotical methos for detecting and describing getalt clusters,” IEEE Trans. Computers, vol. 20, no. 1, pp. 68-86, 1971
- ↑ Y. Xu, V. Olman and E. C. Uberbacher, “A mentation using minimum spanning trees,” Image and Vission Computing, vol. 15, pp. 47-57,1997.
- ↑ Zhi Min Wang, Yeng Chai Soh, Qing Song, Kang Sim, “Adaptive spatial information-theoretic clustering for image segmentation,” Pattern Recognition, vol. 42, no. 9, pp. 2029-2044, 2009.
- ↑ Y. Xu, V. Olman and D. Xu, “Clustering gene expression data using a graph-Theriotic approach: An application of minimum spanning trees,” Bioinformatics, vol. 18, no. 4, pp. 536-545, 2002.
- ↑ Asano, M.K.T., Bhattacharya, B., Yao, F.: Clustering algorithms based on minimum and maximum spanning trees. In: In Proceedings of the fourth annual symposium on Computational Geometry. (1998) 252–257
- ↑ Ye, B., Chao, K.M.: Spanning Trees and Optimization Problems. Chapman and Hall (2004)
- ↑ Clustering with Minimum Spanning Trees Yan Zhou, Oleksandr Grygorash, Thomas F. Hain, Meador Inge Zach Jorgensen School of Computer and Information Sciences University of South Alabama, Mobile, AL 36688 USA Urban Insight Inc. 5657 Wilshire Blvd Ste 290, Los Angeles, CA 90036
- ↑ Attribute Clustering Based on Heuristic Tree Partition Jorge Cordero H. and Yifeng Zeng Department of Computer Science Aalborg University 9220 Aalborg, Denmark
- ↑ Томас Кормен и др. Алгоритмы: построение и анализ = INTRODUCTION TO ALGORITHMS. — 2-е изд. — М.: «Вильямс», 2006. — С. 1296. — ISBN 0-07-013151-1. (Глава 21. Раздел 21.3: "Леса непересекающихся множеств)
