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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• [math]M = \{ x_{m_1}, x_{m_2}, ... , x_{m_k} \}[/math] - множество медоидов;
• [math]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[/math] - искомые кластеры.

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

• фаза BUILD (состоит из [math]k[/math] шагов последовательного выбора [math]k[/math] медоидов):
  • [math]O_0 = X, M_0 = \varnothing ;[/math]
  • [math]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}\} ; [/math]
  • [math]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}\} ; [/math]
  • [math]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}\} ; [/math]
  • ...
  • [math]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}\} ; [/math]

• фаза SWAP (итерационный процесс):
  • начальное приближение: [math]O_0 = O_k^{BUILD}, M_0 = M_k^{BUILD};[/math]
  • [math]i[/math]-й шаг итерации:
    • вычисление целевой функции при оптимальной замене (замене, при которой достигается лучшая кластеризация из возможных на текущей итерации):
      [math]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) ); [/math]
    • определение пары (медоид, не-медоид), на которой достигается оптимальная замена:
      [math] (x_{m_s}, x_{o_h}) = \arg F^i; [/math]
    • проверка критерия останова;
    • изменение множеств в соответствии с выбранной парой [math] (x_{m_s}, x_{o_h}):[/math]
      [math]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}\} ;[/math]
  • критерий останова:
    • [math]T^i = F^{i-1} - F^{i} \leq 0 .[/math]

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

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

[math]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)}.[/math]

Где:

• [math]\overline{MIN}_{1, x_{o_h} \in O_0} = \{ \rho (x_{o_h}, x_j ) \| 1 \leq j \leq n \};[/math]
• [math]\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.[/math]

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

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

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

[math]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) ) =[/math]
[math]= \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) )[/math]

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

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

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

[math]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} [/math], где [math]p, q \in R^w.[/math]

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 };

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

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

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

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

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

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

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

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

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

[math]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[/math]
[math] \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 [/math]
[math] \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} [/math]

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

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

[math]S = ( ( (1 mins)*(n-1) + (n-1) summs )*(n-k) + n*((k-1-1) mins) )*k \Rightarrow[/math]
[math] \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} [/math]

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

Полученные оценки позволяют говорить о том, что 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.