Участник:SenderovichNikita/Алгоритм кластеризации, основанный на построении каркаса

MST-Кластеризация
Последовательный алгоритм
Последовательная сложность	[math]O(N^2\log N)[/math]
Объём входных данных	[math]O(N)[/math]
Объём выходных данных	[math]O(N)[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(\log N)[/math]
Ширина ярусно-параллельной формы	[math]O(N^2)[/math]

Данный документ содержит описание алгоритма кластеризации, основанного на построении минимального остовного дерева (каркаса, остова, кратчайшего незамкнутого пути, минимального покрывающего дерева, англ. Minimum Spanning Tree (MST)) - эвристического графового алгоритма построения разбиения множества объектов на такие непересекающиеся подмножества объектов, что объекты, принадлежащие одному подмножеству, схожи по заданной метрике в большей степени, чем объекты, принадлежащие разным группам.

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

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

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

Рассматривается выборка [math]X = \{x_i\}, i = 1, 2, \ldots, N[/math], состоящая из [math]N[/math] объектов с заданной симметричной неотрицательной функцией расстояния [math]\rho[/math]. Данную выборку можно рассматривать как связный неориентированный взвешенный граф [math]G[/math], множество вершин которого совпадает со множеством объектов выборки, веса рёбер соответствуют значениям функции расстояния между соответствующими вершинами.

Также должно быть задано искомое число кластеров [math]K[/math].

Рис.1. Ключевая идея алгоритма кластеризации на основе построения остовного дерева, разбиение на 3 кластера

Идея рассматриваемого алгоритма кластеризации очень проста:

• с помощью какого-либо алгоритма на определённом выше графе строится минимальное остовное дерево
• из полученного дерева выбрасывается [math]K - 1[/math] ребро с наибольшими весами, что приводит к появлению в остовном дереве [math]K[/math] компонент связности, каждая из которых объявляется кластером

Иллюстрация работы алгоритма приведена на рис. 1. В данном случае рассматриваются 13 объектов на евклидовой плоскости, построено остовное дерево, после чего из него выброшены 2 наиболее длинных ребра и получены три кластера.

Существуют различные модификации описанного алгоритма, связанные с другим алгоритмом конструирования исходного графа или условием выбрасывания рёбер из построенного остовного дерева. Например, число [math]K[/math] может быть полезно определить не изначально, а подобрать, рассматривая отсортированную по убыванию последовательность весов рёбер в остовном дереве: если имеется резкое уменьшение весов, то можно удалить все рёбра до него. Формализацией такого подхода является алгоритм ZEMST. Однако в дальнейшем ограничимся анализом определённого выше простейшего варианта.

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

Пусть задан связный неориентированный граф на множестве объектов [math]G = (X, E)[/math] с неотрицательными весами рёбер [math]\rho(e) = \rho(x_i, x_j) \geqslant 0[/math].

Наиболее частый используемый на практике сценарий следующий:

• каждый объект описывается [math]d[/math]-мерным вектором в пространстве
• используется функция расстояния [math]\ell_p[/math]

Например, случай, когда в качестве метрики выбирается евклидова метрика [math]\ell_2[/math] в литературе называется EMST - Euclidean Minimum Spanning Tree.

Хотя в некоторых задачах возможна ситуация, когда данный граф не является полным (по каким-то причинам определены не все расстояния между вершинами), но в описанных типовых случаях граф является полным, содержит [math]\frac{N(N - 1)}{2} = O(N^2)[/math] рёбер.

1.2.1 Обзор алгоритмов поиска минимального остова

Для поиска минимального остовного дерева на практике используется следующие три похожих друг на друга алгоритма и их комбинации:

• Алгоритм Прима

• Алгоритм Крускала

• Алгоритм Борувки

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

В алгоритме Прима рассматривается текущее остовное дерево [math]T[/math]. Изначально [math]T[/math] состоит из одной вершины, выбранной произвольно. Далее на каждом шаге отыскивается минимальное ребро, ровно один из концов которого принадлежит [math]T[/math], оно добавляется в [math]T[/math] вместе со вторым своим концом. Алгоритм завершается через [math]N - 1[/math] шаг - ровно столько рёбер в дереве на [math]N[/math] вершинах. Последовательная реализация, при которой минимальное ребро на каждом шаге ищется полным перебором рёбер, исходящих из [math]T[/math], имеет сложность [math]O(N^3)[/math] на полном графе. Алгоритм можно ускорить, если поддерживать двоичное дерево поиска по рёбрам, ровно один из концов которых принадлежит [math]T[/math]. Тогда при добавлении очередной вершины в [math]T[/math] необходимо добавить в двоичное дерево не более чем [math]N - 1[/math] новое ребро за [math]O(N \log N)[/math], а выбор минимума осуществляется за [math]O(\log N)[/math]. Поскольку шаг добавления вершины повторяется [math]N[/math] раз, итоговая сложность [math]O(N^2 \log N)[/math].

В алгоритме Крускала остов [math]T[/math] изначально пуст, и в остов на каждом шаге добавляется вместе со своими концами минимальное ребро, имеющееся в графе, не нарушающее свойства ацикличности уже выбранного на предыдущих шагах леса [math]T[/math]. Через [math]N - 1[/math] шаг лес [math]T[/math] гарантированно станет связным деревом. Если искать минимум полным перебором по [math]O(N^2)[/math] рёбрам полного графа, сложность получится [math]O(N^3)[/math]. Если же предварительно отсортировать все рёбра по весу за [math]O(N^2 \log N^2) = O(N^2 \log N)[/math], то дальше все рёбра искомого остова можно будет добавить в процессе одного прохода по этому массиву рёбер в порядке возрастания весов - для каждого ребра проверяем, не нарушает ли оно свойство ацикличности [math]T[/math] и добавляем его в [math]T[/math], если не нарушает. Для того чтобы быстро проверять свойство принадлежности обоих концов ребра одной из компонент [math]T[/math], используется эффективная структура данных - система непересекающихся множеств. Сложность одной такой проверки [math]O(\alpha(N))[/math], где [math]\alpha(N)[/math] - обратная функция Аккермана, значение которой не превосходит 4 для [math]N \leqslant 10^{500}[/math]. Итоговая сложность такого алгоритма равна сложности сортировки и также составляет [math]O(N^2 \log N)[/math].

В алгориме Борувки, как и в алгоритме Крускала, поддерживается лес [math]T[/math]. Изначально лес состоит из [math]N[/math] изолированных вершин графа. На каждом шаге путём перебора всех рёбер графа для каждого дерева в [math]T[/math] ищется минимальное ребро, ровно один конец которого принадлежит этому дереву. После этого все найденные рёбра добавляются в [math]T[/math]. Поскольку после каждой итерации число деревьев в лесе [math]T[/math], исходно равное [math]N[/math], уменьшается не менее чем в 2 раза, таких шагов будет произведено не более [math]\log N[/math], поэтому итоговая сложность алгоритма составляет [math]O(N^2\log N)[/math].

1.2.2 Анализ ресурса параллелизма алгоритмов поиска минимального остова

Итак, каждый из данных алгоритмов имеет сложность [math]O(N^2\log N)[/math], однако с точки зрения параллельной реализации они различаются. Эффективные реализации алгоритмов Прима и Крускала используют структуры данных, рассчитанные на последовательный характер алгоритма. Так, двоичное дерево поиска, содержащее все рёбра, исходящие из построенного на данный момент остова [math]T[/math], трудно приспособить к распределённым архитектурам. А при работе с полным графом с большим числом вершин хранить его становится практически невозможно. В алгоритме Крускала можно распараллелить этап сортировки рёбер, но на последнем шаге алгоритма необходимо последовательно просмотреть отсортированный массив [math]O(N^2)[/math] рёбер, поддерживая СНМ, что сводит на нет параллелизм. Таким образом, для этих двух алгоритмов необходима существенная переработка для получения эффективного параллельного решения .

На этом фоне выгодно выделяется исторически первый алгоритм Борувки. Он был разработан чешским учёным Отакаром Борувкой в 1926 году для решения задачи построения оптимальной электрической сети в Моравии. Он эффективнее всего использует тат факт, что остовные деревья сперва можно найти для отдельных частей графа, после чего объединить полученные деревья. Далее будем рассматривать именно этот метод для построения минимального остовного дерева в алгоритме кластеризации.

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

Вычислительное ядро всех алгоритмов поиска каркаса - нахождение минимального из рёбер, не приводящего к возникновению ацикличности графа. Именно для оптимизации этой операции используются эффективные структуры данных и алгоритмы. В дальнейшем будет показано, что для выбранного алгоритма Борувки на этом этапе возможна эффективная параллельная реализация. В случае использования этого алгоритма для поиска минимального остовного дерева вычисление ядра состоит в нахождении для каждого дерева из уже имеющегося остовного леса минимального исходящего из него ребра, т.е. минимального ребра исходящего из вершин данного дерева.

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

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

• вычисление расстояния [math]\rho(x_i, x_j)[/math] между объектами (вершинами полного графа)
• функции для работы с графом - инициализация добавление, удаление ребра
• использование упомянутого выше алгоритма СНМ (система непересекающихся множеств) для поддержки компонент связности: можно объединять компоненты при добавлении ребра и проверять принадлежность двух вершин одной компоненте за [math]O(1)[/math] (функции компонента, обновить_компоненты)
• процедура сравнения рёбер: рёбра сравниваются по весу, а при равенстве весов по идентификатору ребра. Это важная составная часть алгоритма Борувки, поскольку иначе на этапе добавления рёбер могут возникнуть циклы. Простейшим примером является граф из трёх вершин, все попарные расстояния между которыми равны 1: если по ходу работы для первой окажется ближайшей вторая, для второй третья, а для третьей первая, то в граф будет добавлен цикл.
• сортировка массива рёбер по их весу

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

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

Алгоритм Борувки поиска минимального остовного дерева в полном графе

 1 Вход: Выборка объектов X размера N
 2 Выход: рёбра минимального остовного дерева T
 3 Алгоритм:
 4   T := пустое множество рёбер
 5   инициализировать компоненты числами от 1 до N
 6   есть_разные_компоненты = 1
 7   Пока (есть_разные_компоненты == 1): 
 8     есть_разные_компоненты = 0
 9     инициализировать массив D кратчайших рёбер, исходящих из компонент, бесконечно длинными рёбрами
10     цикл по всем парам вершин (x, y) из X:
11       если компонента(x) не равно компонента(y):
12         есть разные компоненты = 1
13         d = расстояние(x, y)
14         если (ребро (x, y) < ребро D(компонента(x))):
15             D(компонента(x)) = (x, y)
16         если (ребро (x, y) < ребро D(компонента(y))):
17             D(компонента(y)) = (x, y)
18     цикл по рёбрам E из D:
19        обновить_компоненты(E)
20        добавить E в множество T
21   вернуть Т

Итоговый алгоритм кластеризации:

 1 Вход: Выборка объектов X, число кластеров K
 2 Выход: метки кластеров для каждого объекта из X
 3 Алгоритм:
 4    T := рёбра минимального остовного дерева для набора объектов X
 5    E := отсортировать рёбра T по убыванию
 6    цикл по первым K рёбрам в E:
 7       удалить ребро из T
 8    цикл по вершинам v из T:
 9      компоненты[v] = компонента(v)
10    вернуть компоненты

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

В разделе "Обзор алгоритмов поиска минимального остова" уже было показано, что алгоритм Борувки поиска остова имеет сложность [math]O(N^2 \log N)[/math]: поскольку на каждом шаге число компонент уменьшается не менее чем в 2 раза, цикл while выполнится не более [math]\log N[/math] раз. Внутри цикла находятся:

• инициализация массива D кратчайших рёбер, исходящих из компонент связности, бесконечностями - [math]O(N)[/math]
• перебор по всем рёбрам полного графа и обновление содержимого ячеек массива D - [math]O(N^2)[/math]
• добавление найденных рёбер из массива D в граф [math]T[/math] - [math]O(N)[/math]
• обновление компонент связности в графе [math]T[/math] - [math]O(N)[/math]

После поиска остова рёбра остовного дерева сортируются за [math]O(N \log N)[/math] операций, за [math]O(N)[/math] производится выбрасывание рёбер, за [math]O(N)[/math] определяются финальные метки кластеров.

Итоговая сложность [math]O(N^2 \log N)[/math] операций. В ходе работы основными являются операции сравнения и обращения в память. Обратим внимание, что "узким местом" алгоритма является перебор по всем рёбрам для нахождения минимальных рёбер, исходящих из каждой компоненты связности остовного леса [math]T[/math].

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

Приведём информационный граф, оставив за скобками детали реализации последних шагов алгоритма - сортировки и выбрасывания рёбер:

Рис.2. Информационная структура алгоритма кластеризации, основанного на поиске минимального остовного дерева

Кратко опишем шаги схемы. Для вычисления массива D необходимо обработать некоторые из перебираемых рёбер (конкретные связи зависят от структуры текущего остовного дерева [math]T[/math]), поэтому на этом этапе присутствует нерегулярность по данным; кроме того необходима инициализация массива D. Операции обновления структуры остовного графа независимы, могут быть распараллелены и склеены в шаг 3. Алгоритм сортировки теоретически имеет высоту ярусно-параллельной формы [math]O(\log N)[/math]. Независимо могут быть выполнены и последние этапы алгоритма, связанные с выбрасыванием рёбер и получением финального ответа, хотя на практике это не так уж важно.

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

Как показал анализ сложности узким местом является нахождение минимальных рёбер, исходящих из каждой компоненты связности остовного леса [math]T[/math]. Информационный граф говорит о том, что этот этап можно эффективно распаралелить, так как выбор минимального исходящего ребра для каждой компоненты связности остовного дерева и даже для каждой вершины может исполнятся независимо. Таким образом [math]O(N^2)[/math] рёбер - ширина ярусно-параллельной формы алгоритма. Процедуры инициализации массива D, добавления найденных рёбер в граф могут быть произведены независимо по элементам. Параллельно может быть проведена и сортировка, и выбрасывание рёбер в конце, хотя это не так принципиально в реальности. Но это позволяет говорить, что высота ярусно-параллельной формы имеет порядок [math]O(\log N)[/math].

Операции на ключевом этапе алгоритма достаточно сбалансированы - в полном графе для каждой из вершин необходимо просмотреть всех её соседей, поэтому в случае равномерного по вершинам графа распределения работ, эффективность параллельной реализации будет высокой. В целом, имеет место массовый параллелизм на итерациях вычисления ближайших вершин для текущего набора компонент связности [math]T[/math].

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

Входными данными алгоритма является описание выборки из [math]N[/math] объектов. В типовом случае каждый из них описывается [math]d[/math] координатами, однако может описываться каким-то иным способом, позволяющим находить расстояния между объектами. Также на вход передаётся число искомых кластеров [math]K[/math].

Выходом является массив из [math]N[/math] целых чисел от 1 до [math]K[/math] - номера кластеров, найденных алгоритмом кластеризации.

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

Известно, что минимальное остовное дерево определено однозначно, если в графе все рёбра различны. Очень важно ввести операцию сравнения рёбер и добиться, чтобы все рёбра были различны. Пример, когда алгоритм Борувки зацикливается из-за этого, приведён выше.

Распараллеливание описанного алгоритма может принести значительный выигрыш - без него алгоритм фактически невозможно запустить на данных с числом объектов более 100 тысяч. Теоретические возможности распараллеливания на самом сложном этапе - вычислении минимального остовного дерева - были подробно обсуждены выше. Однако на практике может возникнуть ряд трудностей, связанных с эффективным обменом информацией между процессами. Так, всем параллельно работающим процессам необходимо быстро отвечать на вопрос о расположении вершин графа в одной компоненте связности - решения этой проблемы можно добиться, если поддерживать в каждом вычислительном узле данные сведения. Кроме того, важно помнить, что в условиях полного графа каждому узлу необходима и информация обо всех кластеризуемых объектах.

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

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

Как правило, в библиотеках реализуют не собственно алгоритм кластеризации с использованием остовного дерева, а только алгоритм остовного дерева, предоставляя пользователю возможность самостоятельно реализовать функции построения графа и обработки результатов работы алгоритма.

Параллельная реализация алгоритма MST:

• C++, MPI: Parallel Boost Graph Library; функции dense_boruvka_minimum_spanning_tree, boruvka_then_merge, boruvka_mixed_merge сочетают алгоритм Борувки и алгоритм Крускала.

Последовательные реализации доступны во многих библиотеках:

• NetworkX MST - реализация в библиотеке для работы с графами на Python
• Алгоритм Прима, Алгоритм Крускала - C++ реализация в boost
• Алгоритм Прима, Алгоритм Крускала - Java реализация в JGraphT