Участник:Avasilenko/Partitioning Around Medoids (PAM): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
Строка 1: Строка 1:
−
{{algorithm
+
Авторы: [[Участник:Avasilenko|А.Э.Василенко]], [[Участник:Anastasia|А.В.Тузикова]]
−
| name              = Алгоритм кластризации, основанный на разделении вокруг медоидов (Partitioning Around Medoids, PAM)
 
−
| serial_complexity = <math>??? N * O(n^3)</math>
 
−
| pf_height        = <math>??? Unknown</math>
 
−
| pf_width          = <math>??? Unknown</math>
 
−
| input_data        = <math>??? n^2</math>
 
−
| output_data      = <math>??? n</math>
 
−
}}
 
  
−
Основные авторы описания: [[Участник:Avasilenko|А.Э.Василенко]], [[Участник:Anastasia|А.В.Тузикова]]
+
==Свойства и структура алгоритма==
  
−
== Свойства и структура алгоритма ==
+
=== Общее описание алгоритма ===
 +
 +
Объединение схожих объектов в группы - один из важнейших видов человеческой деятельности<ref>Rousseeuw P. J., Kaufman L. Finding Groups in Data. – Wiley Online Library, 1990.</ref>. Более того - это важная часть процесса познания мира, в котором мы живем. Мы учим детей различать кошек и собак, воду и песок, мужчин и женщин путем постоянного развития и улучшения подсознательных схем классификации. Примеры есть и во взрослой жизни – в химии необходимо классифицировать соединения, в археологии – объединять находки по периодам и т.д.
 +
 
 +
 
 +
'''Задача кластеризации'''  - разделение множества объектов на заданное число множеств с минимизацией некоторого функционала - является популярной задачей, относится к разделу машинного обучения и находит применение в таких областях, как анализ текста, биоинформатика, интеллектуальные транспортные системы и др. <ref> Речкалов Т.В. Параллельный алгоритм кластеризации для многоядерного сопроцессора Intel Xeon Phi // Суперкомпьютерные дни в России: Труды vеждународной конференции. – М.: Изд-во МГУ, 2015, c. 532.</ref>
 +
 
 +
'''PAM''' ('''Partitioning Around Medoids''') - плоский алгоритм кластеризации, позволяющй выделить <math>k</math> кластеров (множеств), удовлетворяющих следующим свойствам:
 +
 +
*  каждая группа содержит хотя бы один объект;
 +
*  каждый объект принадлежит ровно одной группе.
 +
 
 +
 
 +
Модель <ref>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.</ref>, лежащая в основе алгоритма 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>
  
−
=== Общее описание алгоритма ===
+
Процедура изменения множества медоидов повторяется, пока есть возможность улучшения значения целевой функции.
 +
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
 +
 +
''Исходные данные'':
 +
 +
*  множество кластеризуемых объектов <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>-й шаг итерации:
 +
*** вычисление целевой функции при оптимальной замене (замене, при которой достигается лучшая кластеризация из возможных на текущей итерации): <br><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>
 +
*** определение пары (медоид, не-медоид), на которой достигается оптимальная замена: <br><math> (x_{m_s}, x_{o_h}) = \arg F^i; </math>
 +
*** проверка критерия останова;
 +
*** изменение множеств в соответствии с выбранной парой <math> (x_{m_s}, x_{o_h}):</math> <br><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>
 +
  
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
 +
 +
'''Вычислительное ядро алгоритма 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 < 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><br><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> что позволяет сократить количество вычисляемых минимумов.
 +
  
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
 +
 +
Нахождение расстояния между двумя объектами <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>
 +
  
 
=== Схема реализации последовательного алгоритма ===
 
=== Схема реализации последовательного алгоритма ===
 +
 +
''Псевдокод алгоритма'':
 +
 +
'''Вход'''  : Множество объектов <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 > 0</math>) {
 +
6      Обновить множества <math>M, O</math>;
 +
7  } '''else''' {
 +
8      '''break''';
 +
9  }
 +
10 };
 +
 +
''Оптимизации'':
 +
 +
* '''Предвычисление матрицы расстояний'''. Необходимо предварительно вычислить матрицу расстояний и сохранить её в памяти, так как в ходе вычисления алгоритма многократно используются всевозможные расстояния между элементами, а вычисление каждого из них является ресурсоёмкой операцией из нескольких умножений и одного извлечения квадратного корня.
 +
 +
* '''Минимизация кеш-промахов'''. Матрица расстояний, к которой постоянно выполняется доступ на чтение для вычисления минимумов и их суммирования, должна храниться в памяти как матрица размера <math>n*n</math>, несмотря на то, что она симметрическая и достаточно хранить лишь половину матрицы (элементы над или под главной диагональю). Данная оптимизация потребует использования в 2 раза большего количества памяти для хранения расстояний (а это является главным расходом памяти программы, так как все остальные используемые данные являются временными (локальными для каждой операции) и представляют из себя несколько векторов длинной <math>n</math> элементов). Однако так как в данном алгоритме доступ к матрице осуществляется построчно (в  [[#Вычислительное ядро алгоритма|вычислительном ядре алгоритма]] поэлементно ищутся  минимумы 2-х векторов, после чего все они складываются), то в этом случае из матрицы можно будет читать данные по строкам подряд, что снижает количество кеш-промахов. Использование только половины матрицы приводит к частым промахам в кеше в связи с необходимостью часть векторов загружать по строкам, а часть - по столбцам.
 +
  
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
  
−
=== Информационный граф ===
+
Для вычисления '''матрицы расстояний''' требуется:
 +
* <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> <br>
 +
<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> <br>
 +
<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> <br>
 +
<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> как по количеству найденных минимумов, так и по количеству операция сложения.
  
−
'''Входные данные''': ????? плотная квадратная матрица <math>A</math> (элементы <math>a_{ij}</math>).
+
Полученные оценки позволяют говорить о том, что PAM относится к алгоритмам '''с квадратичной сложностью'''.
  
−
'''Объём входных данных''':  ????? <math>n^2</math>.
+
=== Информационный граф ===
  
−
'''Выходные данные''':  ????? n вещественных собственных чисел <math> | l_{i} | </math> матрицы <math>A</math>
+
=== Ресурс параллелизма алгоритма ===
  
−
'''Объём выходных данных''':  ????? <math>n</math>.
+
=== Входные и выходные данные алгоритма ===
  
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
−
// TODO
 
  
 
== Программная реализация алгоритма ==
 
== Программная реализация алгоритма ==

Версия 10:59, 12 октября 2016

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

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

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

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


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

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

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


Модель [3], лежащая в основе алгоритма 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 относится к алгоритмам с квадратичной сложностью.

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

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

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

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

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

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

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

3 Литература

  1. Нейский И. М. Классификация и сравнение методов кластеризации //ББК 32.813 И 76 Составитель: ЮН Филиппович. – 2006. – С. 130.
    1. ↑ Rousseeuw P. J., Kaufman L. Finding Groups in Data. – Wiley Online Library, 1990.
    2. ↑ Речкалов Т.В. Параллельный алгоритм кластеризации для многоядерного сопроцессора Intel Xeon Phi // Суперкомпьютерные дни в России: Труды vеждународной конференции. – М.: Изд-во МГУ, 2015, c. 532.
    3. ↑ 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.
    Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Avasilenko/Partitioning_Around_Medoids_(PAM)&oldid=5594