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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ 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]\displaystyle{ \overline{MIN}_{1, x_{o_h} \in O_0} = \{ \rho (x_{o_h}, x_j ) \| 1 \leq j \leq n \}; }[/math]
• [math]\displaystyle{ \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]\displaystyle{ i }[/math]-м шаге целевую функцию нужно вычислить лишь [math]\displaystyle{ |O_{i-1}| = n-i+1 }[/math] раз.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ 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]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ O(kn^2) }[/math] как по количеству найденных минимумов, так и по количеству операций сложения.

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

[math]\displaystyle{ S = ( ( (1 mins)*(n-1) + (n-1) summs )*(n-k) + n*((k-1-1) mins) )*k \Rightarrow }[/math]
[math]\displaystyle{ \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]\displaystyle{ O(kn^2) }[/math] как по количеству найденных минимумов, так и по количеству операция сложения.

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

Литература

1. Нейский И. М. Классификация и сравнение методов кластеризации //ББК 32.813 И 76 Составитель: ЮН Филиппович. – 2006. – С. 130.