Участник:Екатерина/Алгоритм устойчивой кластеризации с использованием связей

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

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

Задача кластеризации заключается в следующем. Имеется множество [math]\displaystyle{ n }[/math] объектов [math]\displaystyle{ S=\{s_1,\ldots,s_n\} }[/math] и функция расстояния/сходства между объектами [math]\displaystyle{ \rho (x_i,x_j) }[/math]. Требуется разбить множество [math]\displaystyle{ S }[/math] на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из объектов, близких по метрике [math]\displaystyle{ \rho }[/math], а объекты разных кластеров существенно отличались. При этом каждому объекту [math]\displaystyle{ s\in S }[/math] приписывается метка(номер) кластера [math]\displaystyle{ y_i }[/math].

Алгоритм кластеризации - это функция [math]\displaystyle{ f:X\to Y }[/math], которая каждому объекту [math]\displaystyle{ s\in S }[/math] ставит в соответствие метку кластера [math]\displaystyle{ y\in Y }[/math].

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

• эвристические графовые алгоритмы;
• статистические алгоритмы;
• иерархические алгоритмы.

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

Одним из алгоритмов алгомеративной иерархической кластеризации является алгоритм устойчивой кластеризации с использованием связей (robust clustering using links, ROCK), предложенный Sudipto Guha (Stanford University), Rajeev Rastogi (Bell Laboratories) и Kyuseok Shim (Bell Laboratories) в 2000 году. Этот алгоритм отлично подходит для объектов с категориальными признаками, то есть признаками, принимающими свои значения из некоторого множества.

В 2005 году была предложена быстрая версия этого ROCK-алгоритма "Quick ROCK"(QROCK).

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

Входные данные: множество объектов [math]\displaystyle{ P = \{p_1, \ldots , p_n\} }[/math], целое число [math]\displaystyle{ k }[/math], пороговое значение [math]\displaystyle{ \theta }[/math], где [math]\displaystyle{ 0\leq\theta\leq 1 }[/math].

Необходимо разделить множество [math]\displaystyle{ P }[/math] на [math]\displaystyle{ k }[/math] непересекающиеся кластеры.

Определение 1. Два объекта [math]\displaystyle{ p_i }[/math] и [math]\displaystyle{ p_j }[/math] будут соседями, если они значительно близки друг к другу, то есть [math]\displaystyle{ p_i }[/math] и [math]\displaystyle{ p_j }[/math] являются соседями, если [math]\displaystyle{ sim(p_i, p_j)\geq \theta }[/math].

Мера схожести пары объектов [math]\displaystyle{ sim(p_i, p_j) }[/math] может быть основана на евклидовом расстоянии, коэффициенте Жаккарда или на какой-либо другой неметрической функции схожести. Предполагается, что мера [math]\displaystyle{ sim }[/math] принимает значения от 0 до 1, при этом чем ближе значение меры к 1, тем больше схожи точки. Чаще всего используется мера схожести, основанная на коэффициенте Жаккарда (Jaccard coefficient): [math]\displaystyle{ sim(p_i,p_j) = \frac{|p_i\cap p_j|}{|p_i\cup p_j|} }[/math], где [math]\displaystyle{ |p_i| }[/math] - количество атрибутов в [math]\displaystyle{ p_i }[/math]. Чем больше атрибутов в [math]\displaystyle{ p_i }[/math] и [math]\displaystyle{ p_j }[/math] похожи, тем больше значение [math]\displaystyle{ |p_i\cap p_j| }[/math]. Деление его на [math]\displaystyle{ |p_i\cup p_j| }[/math] гарантирует, что значение [math]\displaystyle{ \theta }[/math] окажется в промежутке от 0 до 1. Таким образом уравнение вычисляет похожесть двух объектов на основе схожести атрибутов, входящих в них.

Определение 2. Назовем функцией связи между объектами функцию [math]\displaystyle{ link(p_i, p_j) }[/math], определяемую как число общих соседей между объектами [math]\displaystyle{ p_i }[/math] и [math]\displaystyle{ p_j }[/math]. Из определения следует, что чем больше значение функции связи, тем вероятнее, что объекты [math]\displaystyle{ p_i }[/math] и [math]\displaystyle{ p_j }[/math] принадлежат одному кластеру.

Определение 3. Назовем функцией связи между кластерами функцию [math]\displaystyle{ link[C_i, C_j] = \sum_{\begin{smallmatrix}p_q\in C_i,\; p_r\in C_j\end{smallmatrix}}^{ } link(p_{q},p_r) }[/math], хранящую число связей между объектами кластеров [math]\displaystyle{ C_i }[/math] и [math]\displaystyle{ C_j }[/math].

Определение 4. Целевая функция качества, использующаяся при выборе кластеров для объединения:

[math]\displaystyle{ g(C_i, C_j) = \frac{link[C_i,C_j]}{(n_i + n_j)^{1 + 2f(\theta )} - n_i^{1+2f(\theta )} - n_j^{1+2f(\theta )}} }[/math], где [math]\displaystyle{ n_i }[/math] - число объектов в кластере [math]\displaystyle{ C_i }[/math], функция [math]\displaystyle{ f(\theta ) }[/math] выбирается пользователем, обычно [math]\displaystyle{ f(\theta ) = \frac{\theta - 1}{\theta + 1} }[/math] , [math]\displaystyle{ n_i^{1+2f(\theta )} }[/math] - ожидаемое число связей между парами объектов кластера [math]\displaystyle{ C_i }[/math].

Определение 5. Симметричная матрица [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n }[/math] называется матрицей связей, где элемент матрицы [math]\displaystyle{ a_{ij} = link(p_i,p_j), i = \overline{1,n}, j = \overline{1,n} }[/math]

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

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

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

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

Как записано и в описании ядра алгоритма, основную часть метода составляют [math]\displaystyle{ O(n^2\log n) }[/math] обновлений содержимого локальной и глобальной куч.

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

Алгоритм иерархической кластеризации ROCK представлен ниже в виде двух процедур: собственно кластеризации и вычисления связей между точками.

Кластеризация

На вход алгоритм принимает некий набор [math]\displaystyle{ n }[/math] точек [math]\displaystyle{ S }[/math]. Кроме того, заданное значение [math]\displaystyle{ k }[/math] кластеров желаемого разбиения. Изначально каждая точка является отдельным кластером. На шаге 1 вычисляется количество связей между всеми парами точек. На шагах 2 и 3 для каждого [math]\displaystyle{ s }[/math] строится локальная куча [math]\displaystyle{ q[s] }[/math]. Локальная куча [math]\displaystyle{ q[s] }[/math] содержит те кластеры [math]\displaystyle{ t }[/math], для которых [math]\displaystyle{ link[s,t] }[/math] отлична от нуля, причем кластеры в куче упорядочены в порядке убывания, на основе целевой функции качества [math]\displaystyle{ g(s,t) }[/math]. На шаге 4 вычисляется глобальная куча [math]\displaystyle{ Q }[/math], содержащая все текущие кластеры, упорядоченные в порядке убывания на основе [math]\displaystyle{ g(s,max(q[s])) }[/math]. Таким образом на каждом шаге лучший кластер в [math]\displaystyle{ Q }[/math] и соответствующий ему максимальный кластер из [math]\displaystyle{ q }[/math] наилучшая пара кластеров для объединения их в один кластер.

На шаге 5 начинает выполняться цикл [math]\displaystyle{ while }[/math], до тех пор, пока в глобальной куче не останется только [math]\displaystyle{ k }[/math] кластеров. На каждой итерации [math]\displaystyle{ while }[/math]-цикла на шагах 6-9 из [math]\displaystyle{ Q }[/math] извлекается лучший кластер [math]\displaystyle{ u }[/math] и соответствующий ему максимальный кластер [math]\displaystyle{ v }[/math] из локальной кучи и удаляется из [math]\displaystyle{ Q }[/math], так как он будет объединен с [math]\displaystyle{ u }[/math], и больше в ней не нужен. Далее кластеры [math]\displaystyle{ u }[/math] и [math]\displaystyle{ v }[/math] объединяются в новый кластер [math]\displaystyle{ w }[/math].

Далее на шагах 10-15 в цикле [math]\displaystyle{ for }[/math] для каждого кластера, который содержит элементы [math]\displaystyle{ u }[/math] и [math]\displaystyle{ v }[/math] в своей локальной куче, они оттуда удаляются и заменяются на новый кластер [math]\displaystyle{ w }[/math]. Для [math]\displaystyle{ w }[/math] создается новая локальная куча.

На 17 шаге новый кластер добывляется в локальную кучу, а на 18 шаге уничтожаются локальные кучи для [math]\displaystyle{ u }[/math] и [math]\displaystyle{ v }[/math].

[math]\displaystyle{ \begin{align} &\mathbf{procedure}\ cluster (S,k) \\ &\mathbf{begin} \\ &1.\ \ link := compute\_links(S) \\ &2.\ \ \mathbf{for\ each}\ s\in S\ \mathbf{do} \\ &3.\ \ \ \ q[s] := build\_local\_heap(link,s) \\ &4.\ \ Q := build\_global\_heap(S,q) \\ &5.\ \ \mathbf{while}\ size(Q)\ \gt k\ \mathbf{do}\ \{ \\ &6.\ \ \ \ u := extract\_max(Q) \\ &7.\ \ \ \ v := max(q[u]) \\ &8.\ \ \ \ delete(Q,v) \\ &9.\ \ \ \ w := merge(u,v) \\ &10.\ \ \ \mathbf{for\ each}\ x\ \in \ q[u]\ \cup \ q[v]\\ &11.\ \ \ \ \ link[x,w] := link[x,u]\ +\ link[x,v] \\ &12.\ \ \ \ \ delete(q[x],u);\ delete(q[x],v) \\ &13.\ \ \ \ \ insert(q[x],w,g(x,w));\ insert(q[w],x,g(x,w) \\ &14.\ \ \ \ \ update(Q,x,q[x]) \\ &15.\ \ \ \} \\ &16.\ \ \ insert(Q,w,q[w]) \\ &17.\ \ \ deallocate(q[u]);\ deallocate(q[v]) \\ &18.\ \ \} \\ &\mathbf{end} \end{align} }[/math]

Вычисление связей

Для вычисления связи между всеми парами точек в [math]\displaystyle{ S }[/math] рассматривается матрица смежности [math]\displaystyle{ A }[/math] размера [math]\displaystyle{ n*n }[/math], где [math]\displaystyle{ A[i,j] }[/math] равна 1 или 0, если [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] являются или не являются соседями соответственно. В результате умножения [math]\displaystyle{ A*A }[/math] в позиции [math]\displaystyle{ A[i,j] }[/math] получим количество связей между [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math]. Для умножения матриц использовался алгоритм Штрассена.

[math]\displaystyle{ \begin{align} &\mathbf{procedure}\ compute\_links(S) \\ &\mathbf{begin} \\ &1.\ \ Compute\ nbrlist[i]\ for\ every\ point\ i\ in\ S \\ &2.\ \ Set\ link[i,j]\ to\ be\ zero\ for\ all\ i,j \\ &3.\ \ \mathbf{for}\ i := 1\ \mathbf{to}\ n\ \mathbf{do}\ \{ \\ &4.\ \ \ \ N := nbrlist[i] \\ &5.\ \ \ \ \mathbf{for}\ j := 1\ \mathbf{to}\ |N| - 1\ \mathbf{do} \\ &6.\ \ \ \ \ \ \mathbf{for}\ l := j+1\ \mathbf{to}\ |N|\ \mathbf{do} \\ &7.\ \ \ \ \ \ \ \ link[N[j],N[l]] := link[N[j],N[l]]+1 \\ &8.\ \ \} \\ &\mathbf{end} \end{align} }[/math]

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

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

• Инициализация глобальной кучи выполняется за время [math]\displaystyle{ O(n) }[/math].
• Инициализация локальной кучи также требует [math]\displaystyle{ O(n) }[/math].
• Связи между парой точек можно вычислить за [math]\displaystyle{ O(n^2 m_\alpha) }[/math] шагов, где [math]\displaystyle{ m_\alpha }[/math] - среднее число связей.
• В процедуре кластеризации основные вычисления происходят в цикле [math]\displaystyle{ while }[/math]. Сам цикл выполняется [math]\displaystyle{ O(n) }[/math] раз. Внутренний [math]\displaystyle{ for }[/math]-цикл доминирует над сложностью в [math]\displaystyle{ while }[/math]-цикле. Его сложность [math]\displaystyle{ O(n\log n) }[/math], так как размер каждой очереди в худшем случае [math]\displaystyle{ n }[/math] и новый объединенный кластер [math]\displaystyle{ w }[/math] может нуждаться в [math]\displaystyle{ O(n) }[/math] локальных очередей. Следовательно весь [math]\displaystyle{ while }[/math]-цикл выполняется за [math]\displaystyle{ O(n^2\log n) }[/math].

Время работы ROCK-алгоритма равна [math]\displaystyle{ O(n^2+nm_m m_\alpha + n^2\log n) }[/math], где [math]\displaystyle{ m_m }[/math] - максимально возможное количество соседей, [math]\displaystyle{ m_\alpha }[/math] - среднее число соседей, [math]\displaystyle{ n }[/math] - количество объектов.

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

Опишем вершины графа алгоритма:

• [math]\displaystyle{ s_i }[/math] - входные данные.
• [math]\displaystyle{ cml }[/math] - процедура, вычисляющая матрицу связей.
• [math]\displaystyle{ blh }[/math] - функция, строящая начальную локальную очередь для объекта [math]\displaystyle{ s_i }[/math].
• [math]\displaystyle{ bgh }[/math] - строит начальную глобальную очередь входных объектов.
• [math]\displaystyle{ ex\_max }[/math] - достает первый кластер из глобальной кучи.
• [math]\displaystyle{ max }[/math] - достает первый кластер из локальной кучи некоторого кластера [math]\displaystyle{ s_i }[/math].
• [math]\displaystyle{ delete }[/math] - удаляет кластер из кучи.
• [math]\displaystyle{ merge }[/math] - производит слияние двух кластеров.
• [math]\displaystyle{ link\_sum }[/math] - сумма двух чисел.
• [math]\displaystyle{ delete }[/math] - удаляет кластер из кучи.
• [math]\displaystyle{ insert }[/math] - добавляет в кучу новый элемент учитывая целевую функцию качества.
• [math]\displaystyle{ update }[/math] - обновление кучи учитывая целевую функцию качества.
• [math]\displaystyle{ deallocate }[/math] - удаление кучи.

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

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

В частности можно выполнять параллельно построение начальной локальной кучи (функция [math]\displaystyle{ build\_local\_heap }[/math]). Т.е. для входа размерности [math]\displaystyle{ n }[/math] и [math]\displaystyle{ p }[/math] процессов, каждый посчитает [math]\displaystyle{ \frac{n}{p} }[/math] локальных куч.

Еще можно выполнять параллельно цикл, связанный обновлением матрицы связей, локальной и глобальной кучи после проведения операции слияния. Но после каждой итерации необходимо делать барьерную синхронизацию для операции [math]\displaystyle{ update }[/math], которая обновляет глобальную кучу.

Как видно из пункту 1.6 сложность последовательного алгоритма была равна [math]\displaystyle{ O(n^2+nm_m m_\alpha + n^2\log n) }[/math]. С учетом распараллеливания описанных выше участков кода на [math]\displaystyle{ p }[/math] процессов, теоретически можно получить время выполнения [math]\displaystyle{ O(\frac{n^2}{p}+nm_m m_\alpha + \frac{n^2\log n}{p}) }[/math].

Также стоит отметить, что можно распараллелить процедуру [math]\displaystyle{ compute\_links }[/math], которая умножает матрицы алгоритмом Штрассена. Данный алгоритм имеет смысл выполнять параллельно только если порядок матрицы не менее 1287.

Сама функция [math]\displaystyle{ build\_local\_heap }[/math] по сути представляет задачу сортировки, которую также можно выполнить параллельно согласно http://www.intuit.ru/studies/courses/1156/190/lecture/4958 или М. В. Якобовский "Введение в параллельные методы решения задач: Учебное пособие".

Распараллеливание функций [math]\displaystyle{ compute\_links }[/math] и [math]\displaystyle{ build\_local\_heap }[/math] также может дать некоторое ускорение выполнению алгоритма, но обе функции выполняются ровно один раз.

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

Входные данные: на вход подается [math]\displaystyle{ n }[/math] объектов и [math]\displaystyle{ k }[/math] - количество кластеров, на которое необходимо разделить объекты и некоторое пороговое значение [math]\displaystyle{ \theta }[/math], принимающее значение от 0 до 1.

Задаются входные данные в виде матрицы размера [math]\displaystyle{ n*r }[/math], состоящей из 0 и 1, где [math]\displaystyle{ r }[/math] - количество признаков у объектов.

Выходные данные: Группа кластеризованных данных.

Выходные данные представлены в виде массива из [math]\displaystyle{ n }[/math] элементов, где каждому элементу сопоставлен номер кластера, в который он попал.

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

Отношение последовательной и параллельной сложностей: [math]\displaystyle{ \frac{O(n^2+nm_m m_\alpha + n^2\log n)}{O(\frac{n^2}{p}+nm_m m_\alpha + \frac{n^2\log n}{p})} }[/math].

Алгоритм является детерменированным.

Вычислительная сложность: [math]\displaystyle{ \frac{O(n^2+nm_m m_\alpha + n^2\log n)}{2n+2} }[/math].

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

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

Алгоритм устойчивой кластеризации с использованием связей реализован в библиотеке CBA в июне 2016 года.