Участник:Avasilenko/Partitioning Around Medoids (PAM)

Авторы: А.Э.Василенко, А.В.Тузикова

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

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

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

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

PAM (Partitioning Around Medoids) - плоский алгоритм кластеризации, позволяющй выделить $k$ кластеров (множеств), удовлетворяющих следующим свойствам:

• каждая группа содержит хотя бы один объект;
• каждый объект принадлежит ровно одной группе.

Модель , лежащая в основе алгоритма PAM, основана на выборе $k$ объектов, являющихся характерными точками соответствующего кластера. Таким образом, кластерам ставятся в соответствие принадлежащие им объекты, называемые медоидами, на основании которых распределяются остальные объекты по принципу наибольшего сходства. Cама модель называется моделью k-medoids (не путать с моделью k-median).

PAM состоит из двух фаз: BUILD и SWAP:

• на фазе BUILD выполняется первичная кластеризация, в ходе которой последовательно выбираются $k$ объектов в качестве медоидов;
• фаза SWAP представляет собой итерационный процесс, в ходе которого производятся попытки улучшить множество медоидов. Алгоритм выполняет поиск пары объектов (медоид, не-медоид), минимизирующих целевую функцию при замене, после чего обновляет множество медоидов.

На каждой итерации алгоритма выбирается пара (медоид, не-медоид) такая, что замена медоида $x_{m_i}$ на не-медоид $x_{o_h}$ дает лучшую кластеризацию из возможных. Оценка кластеризации выполняется с помощью целевой функции, вычисляемой как сумма расстояний от каждого объекта до ближайшего медоида:

$E = \sum_{j = 1}^n \min_{1 \leq i \leq k} \rho ( x_{m_i},x_{o_h} )$

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

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

Исходные данные:

• множество кластеризуемых объектов $X = \{ x_1, x_2, ... , x_n \}$, где $x_i$ - кортеж, состоящий из $w$ вещественных чисел;
• симметрическая матрица $D$, элементы матрицы $d_{ij}=d(i,j)$ - вычисленные расстояния между объектами $x_i, x_j$;
• количество кластеров $k \leq n.$

Обозначения:

• $M = \{ x_{m_1}, x_{m_2}, ... , x_{m_k} \}$ - множество медоидов, $M \subseteq X ;$
• $O = \{ x_{o_1}, x_{o_2}, ... , x_{o_r} \}$ - множество не-медоидов, $O \subseteq X,$ где $r=n-k;$
• $d_{ij}=\rho (x_i, x_j )$, где $\rho : X \times X \rightarrow R$ — метрика расстояния;
• $F_{MO} = \sum_{j = 1}^{n} \min_{ x_{m_l} \in M} \rho ( x_{m_l},x_{o_j} )$ - значение целевой функции для заданных множеств $O$ и $M;$
• $T^i$ - изменение целевой функции на $i$-м шаге фазы SWAP.

Выходные данные:

• $M = \{ x_{m_1}, x_{m_2}, ... , x_{m_k} \}$ - множество медоидов;
• $K_{m_i} = \{ x_{o_h} \in O \| x_{m_i} = \arg \min_{x_{m_s} \in M} \rho (x_{o_h}, x_{m_s}) \}, 1 \leq i \leq k$ - искомые кластеры.

Вычислительные формулы метода:

• фаза BUILD (состоит из $k$ шагов последовательного выбора $k$ медоидов):
  • $O_0 = X, M_0 = \varnothing ;$
  • $x_{m_1} = \arg \min_{x_{o_h} \in O_0} \sum_{j = 1}^n \rho \left ( x_{o_h},x_j \right ), O_1 = O_0 \backslash \{x_{m_1}\}, M_1 = M_1 \cup \{x_{m_1}\} ;$
  • $x_{m_2} = \arg \min_{x_{o_h} \in O_1} \sum_{j = 1}^{n} \min ( \rho ( x_{m_1},x_j ), \rho ( x_{o_h},x_j ) ), O_2 = O_1 \backslash \{x_{m_2}\}, M_2 = M_1 \cup \{x_{m_2}\} ;$
  • $x_{m_3} = \arg \min_{x_{o_h} \in O_2} \sum_{j = 1}^{n} \min_{1 \leq l \leq 2} ( \rho ( x_{m_l},x_j ) , \rho ( x_{o_h},x_j ) ) , O_3 = O_2 \backslash \{x_{m_3}\}, M_3 = M_3 \cup \{x_{m_3}\} ;$
  • ...
  • $x_{m_k} = \arg \min_{x_{o_h} \in O_{k-1}} \sum_{j = 1}^{n} \min_{ x_{m_l} \in M_{k-1}} ( \rho ( x_{m_l},x_j ), \rho ( x_{o_h},x_j ) ) , O_k = O_{k-1} \backslash \{x_{m_k}\}, M_k = M_{k-1} \cup \{x_{m_k}\} ;$

• фаза SWAP (итерационный процесс):
  • начальное приближение: $O_0 = O_k^{BUILD}, M_0 = M_k^{BUILD};$
  • $i$-й шаг итерации:
    • вычисление целевой функции при оптимальной замене (замене, при которой достигается лучшая кластеризация из возможных на текущей итерации):
      $F^i = \min_{x_{m_s} \in M_{i-1} , x_{o_h} \in O_{i-1}} \sum_{j = 1}^{n} \min_{ x_{m_l} \in M_{i-1} \backslash \{x_{m_s}\} } ( \rho ( x_{m_l},x_j ), \rho ( x_{o_h},x_j) );$
    • определение пары (медоид, не-медоид), на которой достигается оптимальная замена:
      $(x_{m_s}, x_{o_h}) = \arg F^i;$
    • проверка критерия останова;
    • изменение множеств в соответствии с выбранной парой $(x_{m_s}, x_{o_h}):$
      $M_i = M_{i-1} \backslash \{x_{m_s}\} \cup \{x_{o_h}\} , O_i = O_{i-1} \cup \{x_{m_s}\} \backslash \{x_{o_h}\} ;$
  • критерий останова:
    • $T^i = F^{i-1} - F^{i} \leq 0 .$

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

Вычислительное ядро алгоритма PAM на этапе BUILD состоит из множественных вычислений значений целевой функции при различных вариантах выбора дополнительного медоида (шаг алгоритма $1 \leq i \leq k$):

$F^i = \min_{x_{o_h} \in O_{i-1}} \sum_{j = 1}^{n} \min_{ x_{m_l} \in M_{i-1}} ( \rho ( x_{m_l},x_j ), \rho ( x_{o_h},x_j ) ) = \min_{x_{o_h} \in O_{i-1}} \sum_{j = 1}^{n} \overline{MIN}_{i, x_{o_h}}^{(j)}.$

Где:

• $\overline{MIN}_{1, x_{o_h} \in O_0} = \{ \rho (x_{o_h}, x_j ) \| 1 \leq j \leq n \};$
• $\overline{MIN}_{i, x_{o_h} \in O_{i-1}} = \{ \min ( \overline{MIN}_{i-1, x_{m_{i-1}}}^{(j)}, \rho (x_{o_h}, x_j) ) \| 1 \leq j \leq n \}, \forall i: 1 \lt i \leq k.$

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

С каждой итерацией количество вариантов выбора дополнительных медоидов уменьшается на единицу, таким образом, на $i$-м шаге целевую функцию нужно вычислить лишь $|O_{i-1}| = n-i+1$ раз.

Вычислительное ядро алгоритма PAM на этапе SWAP состоит из множественных вычислений значений целевой функции при различных вариантах замены некоторого существующего медоида на элемент из числа не-медоидов ($i$-я итерация фазы SWAP):

$F^i = \min_{x_{m_s} \in M_{i-1} , x_{o_h} \in O_{i-1}} \sum_{j = 1}^{n} \min_{ x_{m_l} \in M_{i-1} \backslash \{x_{m_s}\} } ( \rho ( x_{m_l},x_j ), \rho ( x_{o_h},x_j) ) =$
$= \min_{x_{m_s} \in M_{i-1} } \min_{x_{o_h} \in O_{i-1}} \sum_{j = 1}^{n} \min ( \min_{ x_{m_l} \in M_{i-1} \backslash \{x_{m_s}\} } ( \rho ( x_{m_l},x_j )) , \rho ( x_{o_h},x_j) )$

Таким образом, при поиске пары (медоид, не-медоид), наилучшим образом минимизирующей целевую функцию, происходит поиск минимума по всем возможным элементам $x_{o_h}$, для которых вычисляется $\sum_{j = 1}^{n} \min (a, b)$, где $a$ есть $\min_{ x_{m_l} \in M_{i-1} \backslash \{x_{m_s}\} } ( \rho ( x_{m_l},x_j )),$ и является инвариантом по отношению к $\min_{x_{o_h} \in O_{i-1}} \sum_{j = 1}^{n} \min ( \min_{ x_{m_l} \in M_{i-1} \backslash \{x_{m_s}\} } ( \rho ( x_{m_l},x_j )) , \rho ( x_{o_h},x_j) ),$ что позволяет сократить количество вычисляемых минимумов.

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

Нахождение расстояния между двумя объектами $\rho \left ( p, q \right )$ является макрооперацией в том случае, если на вход алгоритму подавалась не матрица расстояний, а набор точек с их координатами в пространстве $R^w$. В этом случае может быть использована, например, евклидова метрика:

$d(p,q)=\sqrt{(p_1-q_1)^2+(p_2-q_2)^2+\dots+(p_w-q_w)^2} = \sqrt{\sum_{j=1}^w (p_j-q_j)^2}$, где $p, q \in R^w.$

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

Псевдокод алгоритма:

Вход  : Множество объектов [math]X[/math] или матрица расстояний [math]D[/math], количество кластеров [math]k[/math]
Выход : Множество медоидов [math]M[/math]
0  Если на вход было подано множество объектов, необходимо предварительно вычислить матрицу расстояний;
1  Инициализация множеств [math]M, O[/math]; // фаза BUILD
2  while (1) { // фаза SWAP
3  	Вычислить [math]F^i[/math] и соответствующие ей [math](x_{m_s}, x_{o_h})[/math];// вычисление проводить с учётом приведённой ранее оптимизации (см. 1.3 Вычислительное ядро алгоритма);
4  	Вычислить [math]T^i[/math];
5  	if ([math]T^i \gt  0[/math]) {
6      	Обновить множества [math]M, O[/math];
7  	} else {
8      	break;
9  	}
10 };

Оптимизации:

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

• Минимизация кеш-промахов. Матрица расстояний, к которой постоянно выполняется доступ на чтение для вычисления минимумов и их суммирования, должна храниться в памяти как матрица размера $n*n$, несмотря на то, что она симметрическая и достаточно хранить лишь половину матрицы (элементы над или под главной диагональю). Данная оптимизация потребует использования в 2 раза большего количества памяти для хранения расстояний (а это является главным расходом памяти программы, так как все остальные используемые данные являются временными (локальными для каждой операции) и представляют из себя несколько векторов длинной $n$ элементов). Однако так как в данном алгоритме доступ к матрице осуществляется построчно (в вычислительном ядре алгоритма поэлементно ищутся минимумы 2-х векторов, после чего все они складываются), то в этом случае из матрицы можно будет читать данные по строкам подряд, что снижает количество кеш-промахов. Использование только половины матрицы приводит к частым промахам в кеше в связи с необходимостью часть векторов загружать по строкам, а часть - по столбцам.

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

Для вычисления матрицы расстояний требуется:

• $\frac{n*(n-1)}{2}$ операций вычисления расстояний;
• сложность каждой операции вычисления расстояния (для случая использования Евклидовой метрики) равна $2*w-1$ операций сложения и вычитания, $w$ операций возведения в квадрат и $1$ операции извлечения квадратного корня.

На основании формулы из описания ядра на этапе BUILD определим сложность выполнения всего этапа. Обозначим:

• $summs$ - операция сложения двух вещественных чисел;
• $mins$ - операция нахождения минимума для двух вещественных чисел.

Тогда, количество операций равно:

$S = ( ( (n-1) summs ) * (n-1+1) + (n-1+1-1)mins) + \sum_{i = 2}^{k} ( ( (1 mins)*(n-1) + (n-1) summs ) * (n-i+1) + (n-i+1-1)mins) \Rightarrow$
$\Rightarrow S = \begin{cases} S(mins) = ( ( 0 ) * (n-1+1) + (n-1+1-1)mins) + \sum_{i = 2}^{k} ( ( (1 mins)*(n-1)) * (n-i+1) + (n-i+1-1)mins ) & \quad\\ S(summs) = \sum_{i = 1}^{k} ( ( (n-1) summs ) * (n-i+1) ) & \quad\\ \end{cases} \Rightarrow$
$\Rightarrow S = \begin{cases} S(mins) = n^2(k-1) - \frac{1}{2}(nk^2 - nk) + n - k & \quad\\ S(summs) = n^2k - \frac{1}{2}(nk^2 + nk - k^2 + k ) & \quad\\ \end{cases}$

Таким образом, сложность фазы BUILD есть $O(kn^2)$ как по количеству найденных минимумов, так и по количеству операций сложения.

На основании формулы из описания ядра на этапе SWAP определим сложность вычисления целевой функции $F^i$ при оптимизации процесса нахождения минимума (см. раздел Вычислительное ядро алгоритма). Количество операций равно:

$S = ( ( (1 mins)*(n-1) + (n-1) summs )*(n-k) + n*((k-1-1) mins) )*k \Rightarrow$
$\Rightarrow S = \begin{cases} S(mins) = ( ((1 mins)*(n-1) )*(n-k) + n*((k-1-1) mins) )*k = kn^2 - 3nk + k^2 & \quad\\ S(summs) = ( ( (n-1) summs )*(n-k) )*k = (n-1)(n-k)k & \quad\\ \end{cases}$

Таким образом, сложность одной итерации фазы SWAP есть $O(kn^2)$ как по количеству найденных минимумов, так и по количеству операция сложения.

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

3 Литература

1. Нейский И. М. Классификация и сравнение методов кластеризации //ББК 32.813 И 76 Составитель: ЮН Филиппович. – 2006. – С. 130.
2. Rousseeuw P. J., Kaufman L. Finding Groups in Data. – Wiley Online Library, 1990.
3. Речкалов Т.В. Параллельный алгоритм кластеризации для многоядерного сопроцессора Intel Xeon Phi // Суперкомпьютерные дни в России: Труды vеждународной конференции. – М.: Изд-во МГУ, 2015, c. 532.
4. Kaufman, L. and Rousseeuw, P.J., Clustering by means of Medoids. In: Y. Dodge and North-Holland, editor. Statistical Data Analysis Based on the L1-Norm and Related Methods. Springer US; 1987, p. 405–416.