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

Авторы: Зайцева М.Ф., Ковальчук А.А..

PAM (Partitioning Around Medoids) –– алгоритм разбиения множества объектов на кластеры, базирующийся на использовании наиболее представительных объектов каждого кластера (медоидов) для их описания. Был разработан в 1987 году Леонардом Кауфманом и Питером Русивом . Показывает хорошие результаты разбиения на небольших наборах данных, однако имеет высокую вычислительную сложность, ввиду чего неэффективно работает на больших наборах данных.

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

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

Данный алгоритм используется для решения задачи кластеризации множества объектов на заранее известное число кластеров [math]k[/math]. Каждый кластер характеризуется наиболее представительным (центральным ) объектом –– медоидом. Принадлежность объекта к соответствующему кластеру определяется индексом наиболее близкого к нему медоида.

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

Алгоритм состоит из двух последовательных этапов:

1. Шаг BUILD. Построение начального набора [math]k[/math] медоидов путем их последовательного выбора. Первый из них выбирается как объект, сумма дистанций которого до всех других объектов минимальна (центральный объект). Каждый следующий выбирается из оставшихся объектов как лучше всего минимизирующий сумму расстояний от объектов до ближайших к ним медоидов. Выбранный набор медоидов характеризует "конфигурацию" кластеров, для которой может быть вычислена ее стоимость.

2. Шаг SWAP. На данном шаге производится попытка улучшить стоимость текущей конфигурации путем перебора всех возможных замен (медоид-немедоид). Среди рассматриваемых конфигураций ищется наиболее оптимальная. В случае, если данная конфигурация по стоимости лучше текущей, она принимается в качестве текущей и шаг повторяется с начала, иначе алгоритм завершается.

Преимущества данного алгоритма по сравнению с подобными ему:

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

Недостатком данного алгоритма является неэффективность при больших наборах данных ввиду необходимости полного перебора.

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

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

• [math]o = \{1,\dots,n\} \subset \mathbb{Z}[/math] –– множество индексов объектов
• [math]O = \{O_1,\dots,O_n\}[/math] –– множество объектов
• [math]k \leq n[/math] –– количество кластеров для разбиения
• [math]d:O\times O\rightarrow \mathbb{R}^{1}[/math] –– метрика расчета расстояний между объектами, [math]d_{ij} = d(O_i,O_j)[/math]
• [math]D\in \mathbb{R}^{n\times n}[/math] –– симметрическая матрица дистанций объектов, где [math]D=(d_{ij}), i \in o, j \in o[/math]

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

• [math]m = \{m_1,\dots,m_k\}\subset o[/math] –– множество индексов объектов, принятых в качестве медоидов

Принадлежность объекта к определенному кластеру определяется как индекс ближайшего к нему медоида:

[math]g_m(i)=\arg \min_{m_j\in m} d_{im_j}[/math].

Основная задача состоит в выборе такого множества индексов медоидов [math]m[/math], наилучшим образом минимизирующего функционал:

[math]F_o(m) = \sum_{j \in o} d_{jg_m(j)}[/math].

Введем обозначения:

• [math]D_j = d_{jg_m(j)}[/math] –– индекс первого по близости к j медоида
• [math]E_j = d_{jg_{m \backslash g_m(j)}(j)}[/math] –– индекс второго по близости к j медоида

Решение задачи состоит из следующих шагов:

1. Шаг BUILD

 1.1. Выбор первого медоида

[math]m=\{ \arg \min_{i \in o} \sum_{j \in o} d_{ij} \}[/math]

 1.2 Последовательный выбор оставшихся медоидов до тех пор, пока [math]\left\vert{m}\right\vert \neq k[/math]

[math]m = m \cup \{ \arg \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0) \}[/math]

2. Шаг SWAP

 2.1. Производится перебор всех возможных пар медоид-немедоид [math](i, h): i \in m, h \in o \backslash m[/math], для каждой из которых рассчитывается стоимость смены текущей конфигурации [math]T_m(i,h)[/math]. При смене конфигурации множество медоидов изменяется [math]m = m \cup \{h\} \backslash \{i\}[/math].

[math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}[/math], где [math]C_{jih}[/math] –– вклад объекта [math]j[/math] в перестановку медоида [math]i[/math] и не медоида [math]h[/math].

 При этом [math]C_{jih}[/math] рассчитывается по-разному в зависимости от положения [math]j[/math] относительно [math]i[/math] и [math]h[/math] ([math]j[/math] дальше от [math]i[/math] и от [math]h[/math], чем от своего ближайшего медоида; [math]i[/math] оказался ближайшим медоидом [math]j[/math]; [math]j[/math] оказался дальше от [math]i[/math], но ближе к [math]h[/math], чем к своему медоиду).

[math]C_{jih} = \left\{\begin{matrix} 0, & d_{ji} \gt D_j, d_{jh} \gt D_j \\ d_{jh} - d_{ji}, & d_{ji} = D_j, d_{jh} \lt E_j \\ E_j - D_j, & d_{ji} = D_j, d_{jh} \geqslant E_j \\ d_{jh} - D_j, & d_{ji} \gt D_j, d_{jh} \lt D_j \end{matrix}\right.[/math]

 Среди всех пар находится такая пара [math](i_s, h_s)[/math], что [math]T_m(i_s, h_s)[/math] минимально.

 2.2 В случае, если [math]T_m(i_s, h_s) \geq 0[/math], смена конфигурации медоидов не производится и решение задачи окончено. Иначе производится смена конфигурации и повторяется шаг 2.1.

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

Вычислительное ядро алгоритма на этапе BUILD состоит последовательного выбора такого объекта из числа немедоидов, добавление которого к множеству медоидов максимизирует следующую функцию:

[math]F_{m} = \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0)[/math]

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

Вычислительное ядро алгоритма на этапе SWAP состоит в расчете функции минимизации стоимости изменения конфигурации [math]T_m(i, h)[/math] при замене медоида [math]i[/math] на немедоид [math]h[/math] в текущей конфигурации для всех возможных [math](i, h): i \in m, h \in o \backslash m[/math]:

[math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}[/math].

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

Как было указано в описании ядра алгоритма, основную часть метода занимают:

• вычисление [math](k - 1)[/math] медоида [math]\arg \max_{i \in o \backslash m} \sum_{j \in o \backslash \{m \cup i\} } \max (D_j - d_{ji}, 0)[/math];
• поиск оптимальной замены замены медоида на немедоид среди [math]k(n-k)[/math] кандидатов путем вычисления стоимости замены [math]T_m(i, h) = \sum_{j \in o \backslash \{m \cup h\}} C_{jih}[/math] относительно текущей конфигурации [math]m[/math] для каждого кандидата.

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

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

Вход: матрица дистанций D; количество кластеров k; максимальное количество итераций алгоритма maxIter
Выход: множество медоидов m
1 # Фаза BUILD
2 Выбор стартового объекта в качестве начального медоида: m = {первый медоид}
3 Пока |m| != k:
4     Выбор следующего медоида i_m
5     Добавление найденного медоида в множество: m = m + {i_m}
6 # Фаза SWAP
7 Пока (maxIter--):
8     Лучшая стоимость замены: Tbest = MAXINT
9     Кандидаты на замену: (i_s, h_s) = (-1, -1)
10     Для всех i из множества медоидов:
11         Для всех j из множества немедоидов:
12             Расчет стоимости замены Tm(i,h)
13             Если Tm(i, h) < Tbest:
14                 Tbest = Tm(i, j)
15                 Запомнить найденные кандидаты: (i_s, j_s) = (i, h)
16     Если Tbest >= 0:
17         Выход
18     Иначе:
19         m = m + {h} - {i}

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

Для осуществления шага BUILD потребуется:

• [math]n^2[/math] операций сложения и [math]n[/math] операций сравнения для выбора начального медоида;
• [math]\sum_{l=1}^{k-1} 2(n-l) \cdot (n-l-1)=\frac{2}{3}(k-1)(k^2 - 3kn + k + 3(n-1)n)[/math] операций сложения и вычитания для всех итераций выбора медоидов;
• [math]\sum_{l=1}^{k-1} n-l = \frac{1}{2}(k-1)(2n-k)[/math] операций сравнения для всех итераций выбора медоидов.

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

Для осуществления каждого шага SWAP потребуется:

• [math]2k(n-k)(n-k-1)[/math] операций вычитания и сложения;
• [math]k(n-k)[/math] операций сравнения;

Таким образом, последовательная сложность каждого шага SWAP оценивается как [math]O(kn^2)[/math]. Заранее число итераций шага SWAP неизвестно, поэтому принимаем его равным [math]T[/math]. Отсюда, общая последовательная сложность алгоритма PAM оценивается как [math]O(Tkn^2)[/math].

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

Это очень важный раздел описания. Именно здесь можно показать (увидеть) как устроена параллельная структура алгоритма, для чего приводится описание и изображение его информационного графа (графа алгоритма [1]). Для рисунков с изображением графа будут составлены рекомендации по их формированию, чтобы все информационные графы, внесенные в энциклопедию, можно было бы воспринимать и интерпретировать одинаково. Дополнительно можно привести полное параметрическое описание графа в терминах покрывающих функций [1].

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

В целом, задача изображения графа алгоритма весьма нетривиальна. Начнем с того, что это потенциально бесконечный граф, число вершин и дуг которого определяется значениями внешних переменных, а они могут быть весьма и весьма велики. В такой ситуации, как правило, спасают упомянутые выше соображения подобия, делающие графы для разных значений внешних переменных "похожими": почти всегда достаточно привести лишь один граф небольшого размера, добавив, что графы для остальных значений будут устроены "точно также". На практике, увы, не всегда все так просто, и здесь нужно быть аккуратным.

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

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

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

Рис.1. Информационная структура алгоритма умножения матриц

Рис.2. Информационная структура одного из вариантов алгоритма решения систем линейных алгебраических уравнений с блочно-двухдиагональной матрицей

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

Для осуществления шага BUILD потребуется:

• [math]n[/math] параллельных операций сложения и [math]n[/math] последовательных операций сравнения для выбора начального медоида;
• [math]\sum_{l=1}^{k-1} 2(n-l) = (k-1)(2n-k)[/math] параллельных операций сложения и вычитания для всех итераций выбора медоидов;
• [math]\sum_{l=1}^{k-1} n-l = \frac{1}{2}(k-1)(2n-k)[/math] последовательных операций сравнения для всех итераций выбора медоидов.

Таким образом, параллельная сложность шага BUILD оценивается как [math]O(kn)[/math].

Для осуществления каждого шага SWAP потребуется:

• [math]k(n-k) [/math] параллельных операций вычитания и сложения;
• [math]k(n-k)[/math] последовательных операций сравнения.

Таким образом, параллельная сложность каждого шага SWAP оценивается как [math]O(kn)[/math]. Заранее число итераций шага SWAP неизвестно, поэтому принимаем его равным [math]T[/math]. Отсюда, общая последовательная сложность алгоритма PAM оценивается как [math]O(Tkn)[/math].

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

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

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

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

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

Вычислительная мощность алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. Она показывает, сколько операций приходится на единицу переданных данных. Несмотря на простоту данного понятия, это значение исключительно полезно на практике: чем выше вычислительная мощность, тем меньше накладных расходов вызывает перемещение данных для их обработки, например, на сопроцессоре, ускорителе или другом узле кластера. Например, вычислительная мощность скалярного произведения двух векторов равна всего лишь [math]1[/math], а вычислительная мощность алгоритма умножения двух квадратных матриц равна [math]2n/3[/math].

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

Сбалансированность вычислительного процесса можно рассматривать с разных сторон. Здесь и сбалансированность типов операций, в частности, арифметических операций между собой (сложение, умножение, деление) или же арифметических операций по отношению к операциям обращения к памяти (чтение/запись). Здесь и сбалансированность операций между параллельными ветвями алгоритма. С одной стороны, балансировка нагрузки является необходимым условием эффективной реализации алгоритма. Вместе с этим, это очень непростая задача, и в описании должно быть отмечено явно, насколько алгоритм обладает этой особенностью. Если обеспечение сбалансированности не очевидно, желательно описать возможные пути решения этой задачи.

На практике важна детерминированность алгоритмов, под которой будем понимать постоянство структуры вычислительного процесса. С этой точки зрения, классическое умножение плотных матриц является детерминированным алгоритмом, поскольку его структура при фиксированном размере матриц никак не зависит от элементов входных матриц. Умножение разреженных матриц, когда матрица хранятся в одном из специальных форматов, свойством детерминированности уже не обладает: его свойства, например, степень локальности данных зависит от структуры разреженности входных матриц. Итерационный алгоритм с выходом по точности также не является детерминированным: число итераций, а значит и число операций, меняется в зависимости от входных данных. В этом же ряду стоит использование датчиков случайных чисел, меняющих вычислительный процесс для различных запусков программы. Причина выделения свойства детерминированности понятна: работать с детерминированным алгоритмом проще, поскольку один раз найденная структура и будет определять качество его реализации. Если детерминированность нарушается, то это должно быть здесь описано вместе с описанием того, как недетерминированность влияет на структуру вычислительного процесса.

Серьезной причиной недетерминированности работы параллельных программ является изменение порядка выполнения ассоциативных операций. Типичный пример - это использование глобальных MPI-операций на множестве параллельных процессов, например, суммирование элементов распределенного массива. Система времени исполнения MPI сама выбирает порядок выполнения операций, предполагая выполнение свойства ассоциативности, из-за чего ошибки округления меняются от запуска программы к запуску, внося изменения в конечный результат ее работы. Это очень серьезная проблема, которая сегодня встречается часто на системах с массовым параллелизмом и определяет отсутствие повторяемости результатов работы параллельных программ. Данная особенность характерна для второй части AlgoWiki, посвященной реализации алгоритмов, но вопрос очень важный, и соответствующие соображения, по возможности, должны быть отмечены и здесь.

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

Степень исхода вершины информационного графа показывает, в скольких операциях ее результат будет использоваться в качестве аргумента. Если степень исхода вершины велика, то на этапе реализации алгоритма нужно позаботиться об эффективном доступе к результату ее работы. В этом смысле, особый интерес представляют рассылки данных, когда результат выполнения одной операции используется во многих других вершинах графа, причем число таких вершин растет с увеличением значения внешних переменных.

"Длинные" дуги в информационном графе [1] говорят о потенциальных сложностях с размещением данных в иерархии памяти компьютера на этапе выполнения программы. С одной стороны, длина дуги зависит от выбора конкретной системы координат, в которой расположены вершины графа, а потому в другой системе координат они попросту могут исчезнуть (но не появится ли одновременно других длинных дуг?). А с другой стороны, вне зависимости от системы координат их присутствие может быть сигналом о необходимости длительного хранения данных на определенном уровне иерархии, что накладывает дополнительные ограничения на эффективность реализации алгоритма. Одной из причин возникновения длинных дуг являются рассылки скалярных величин по всем итерациям какого-либо цикла: в таком виде длинные дуги не вызывают каких-либо серьезных проблем на практике.

Для проектирования специализированных процессоров или реализации алгоритма на ПЛИС представляют интерес компактные укладки информационного графа [1], которые также имеет смысл привести в данном разделе.

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

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

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

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

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

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

Рис.5 Масштабируемость классического алгоритма умножения плотных матриц в зависимости от числа процессоров и размера задачи

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

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

1. Последовательные реализации

• В MATLAB (Commercial software)
• ELKI на языке JAVA (Свободное распространение (AGPL))
• В пакете Cluster на языке R (Свободное распространение (GPL))
• В пакете ClusterR на языке R (Свободное распространение (MIT))

2. Параллельные реализации

• В пакете SPRINT на языке R (Свободное распространение (GPL))