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

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

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

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

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

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

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

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

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

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

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

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

Необходимо разделить множество $P$ на $k$ непересекающиеся кластеры.

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

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

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

Определение 3. Назовем функцией связи между кластерами функцию $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)$, хранящую число связей между объектами кластеров $C_i$ и $C_j$.

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

$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 )}}$, где $n_i$ - число объектов в кластере $C_i$, функция $f(\theta )$ выбирается пользователем, обычно $f(\theta ) = \frac{\theta - 1}{\theta + 1}$ , $n_i^{1+2f(\theta )}$ - ожидаемое число связей между парами объектов кластера $C_i$.

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

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

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

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

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

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

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

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

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

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

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

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

На 17 шаге новый кластер добывляется в локальную кучу, а на 18 шаге уничтожаются локальные кучи для $u$ и $v$.

$\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}$

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

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

$\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}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Выходные данные представлены в виде массива из $n$ элементов, где каждому элементу сопоставлен номер кластера, в который он попал.

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

Отношение последовательной и параллельной сложностей: $\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})}$.

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

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

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

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

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