Участник:Антон Тодуа/Partitioning Around Medoids (PAM)

Основные авторы описания: Тодуа А.Р., Гусева Ю.В.

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

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

Алгоритм PAM (Partitioning Around Medoids) был предложен Чарли Кауфаманом (Charlie Kaufman) и Питером Русивом (Peter Rousseeuw) в 1987. Алгоритм является одним из распространённых вариантов реализации [math]k[/math]-medoids.

PAM используется для выполнения кластеризации[1] объектов — процесса автоматического разбиения некоторого множества объектов на кластеры (множества) на основе их схожести (наиболее схожие объекты, как правило, располагаются в одном кластере). В качестве исходных данных алгоритму необходима симметрическая матрица с нолями на главной диагонали (dissimilarity matrix), задающая расстояние (в смысле похожести) между объектами, для которых выполняется кластеризация и число кластеров, по которым следует распределить объекты. Стоит отметить, что не все алгоритмы кластеризации требует задания количества кластеров, например, методы иерархической[2] кластеризации свободны от данного ограничения.

Алгоритма PAM для каждого кластера в качестве его характерного представителя выбирает наиболее подходящий из объектов кластеризации. Выбранные объекты называются медоидами (medoids). Алгоритм состоит из двух стадий. Первая стадия – BUILD предназначена для выбора исходного множества медоидов. Вторая стадия – SWAP является итерационной. На каждой итерации стадии SWAP алгоритм пытается улучшить выбранное множество медоидов с помощью выполнения операции обмена, которая представляет собой замену одного медоида на один из объектов кластеризации. Алгоритм завершает свою работу, когда ни один из возможных обменов не позволяет улучшить разбиение на кластеры. Однако, часто на практике, условием завершения работы алгоритма является момент, когда наилучшее возможное улучшение кластеризации меньше заранее выбранной пороговой величины. Введение данной величины позволяет сократить число итераций, при сравнительно небольших потерях качества кластеризации.

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

В качестве исходных данных алгоритм принимает:

• Симметрическую матрицу [math]D[/math] порядка [math]N[/math] с нолями на главной диагонали для [math]N[/math] объектов кластеризации: [math]i[/math]-я строка матрицы [math]D[/math] определяет расстояние от [math]i[/math]-го объекта до каждого из остальных объектов
• Число кластеров [math]K[/math], при этом: [math]K \lt N[/math]

Для заданных входных данных, алгоритм выбирает [math]K[/math] объектов - медоидов (medoids) из [math]N[/math] объектов кластеризации. Все медоиды относится к разных кластерам и являются их характерными представителями. Кластеры остальных объектов (не медоидов) совпадают с кластером ближайшего к ним медоида. Ближайший к [math]i[/math]-му объекту медоид будем коротко назвать медоидом [math]i[/math]-го объекта. Выбор медоидов осуществляется исходя из минимизации среднего расстояния (average dissimilarity) между объектами внутри одного кластера. От минимизации среднего расстояния можно перейти к минимизации суммы расстояний от объектов до их медоидов. Такой переход не повлияет на выбор медоидов, но позволит не накапливать ошибки округления, возникающие в результате операции деления. Стоит отметить, что в общем случае расстояние между объектами следует понимать как меру похожести (similarity) объектов.

Для дальнейшего описания алгоритма введём следующие обозначения:

• [math]U[/math] – множество объектов кластеризации, не выбранных в качестве медоидов
• [math]S[/math] – множество объектов кластеризации, выбранных в качестве медоидов
• [math]C_{ij}[/math] – изменение расстояния от [math]j[/math]-го объекта до его медоида в результате выбора [math]i[/math]-го объекта в качестве медоида
• [math]T_{ijh}[/math] – изменение расстояния от [math]j[/math]-го объекта до его медоида, в результате обмена [math]i[/math]-го медоида на [math]h[/math]-й объект

Под обменом [math]i[/math]-го медоида на [math]h[/math]-й объект будем понимать операцию, в результате которой [math]i[/math]-й медоид перестаёт быть медоидом, а [math]h[/math]-й объект выбирается в качестве медоида:

[math] \begin{align} S & = S - \{i\} + \{h\} \\ U & = U - \{h\} + \{i\} \end{align} [/math]

Перед выполнением алгоритма множество [math]U[/math]содержит все объекты кластеризации, а множество [math]S[/math] – пусто. Описание значений [math]C_{ij}[/math] и [math]T_{ijh}[/math] и способа их вычисления будет дано при описании макроструктуры алгоритма.

Выполнение алгоритма включает в себя подготовительную стадию BUILD и итерационную стадию SWAP.

Стадия BUILD предназначена для построения начального множества медоидов и состоит в последовательном выполнении следующих шагов:

1. Добавить в множество [math]S[/math] объект, сумма расстояний от которого до всех остальных объектов минимальна (самый центральный объект):

[math]S = S + \{ arg\min_{i \in U} \sum_{j \in U} d_{ij} \}[/math]

2. Выбрать ещё [math]K-1[/math] медоидов следующим образом

[math]S = S + \{ arg\max_{i \in U} \sum_{j \in U - \{i\}} C_{ij} \}[/math]

Стадия SWAP является итерационной, одна итерация состоит в попытке улучшить кластеризацию путём выполнения обмена. На каждой итерации выполняется операция [math]arg\min_{i \in S, h \in U} \sum_{j \in U - \{h\}} T_{ijh}[/math]. Если для найденной пары [math]i \in S[/math] и [math]h \in U[/math] значение [math]\sum_{j \in U - \{h\}} T_{ijh}[/math] меньше ноля (это означает, что кластеризацию можно улучшить), то выполняется обмен [math]i[/math]-го медоида и [math]h[/math]-го объекта. В противном случае дальнейшее улучшение кластеризации с помощью обмена невозможно и алгоритм завершает свою работу.

В результате работы алгоритма множество [math]S[/math] будет содержать ровно [math]K[/math] объектов – медоидов. Как было сказано [math]j[/math]-й объект лежит в том же кластере, что и его медоид: [math]arg\min_{i \in S} d_{ij}[/math]

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

Из математического описания алгоритма следует, что вычислительное ядро алгоритма сосредоточено в итерационной стадии SWAP. Одна итерация состоит в вычислении:

[math]arg\min_{i \in S, h \in U} \sum_{j \in U - \{h\}} T_{ijh}[/math]

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

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

Для описания значений [math]C_{ij}[/math] и [math]T_{ijh}[/math] и способа их вычисления введём следующие обозначения (в дополнение к обозначениям, введённым в математическом описании алгоритма):

• [math]D_{i}[/math] – расстояние от [math]i[/math]-го объекта до его медоида
• [math]E_{i}[/math] – расстояние от [math]i[/math]-го объекта до второго ближайшего к нем медоида

Очевидно, что для любых объектов верно неравенство [math]D_{i} \lt = E_{i}[/math]. Отметим также, что значения [math]D_{i}[/math] и [math]E_{i}[/math] могут меняться при изменениях в множествах [math]U[/math] и [math]S[/math], следовательно они должны пересчитываться после добавления нового медоида или выполнения операции обмена.

Значение [math]C_{ij}[/math] представляет собой число, на которое уменьшается значение [math]D_{j}[/math] в результате выбора [math]i[/math]-го объекта в качестве медоида. Очевидно, что данное значение всегда неотрицательно. Для вычисления [math]C_{ij}[/math] следует рассмотреть два случая. В первом случае, когда [math]D_{j} \lt = d_{ij}[/math], величина [math]D_{j}[/math] не изменится, т.к. ближайший к [math]j[/math]-му объекту медоид останется на том же расстоянии от него, следовательно будет справедливо равенство [math]C_{ij} = 0[/math]. В втором случае, когда [math]D_{j} \gt d_{ij}[/math], ближайшим к [math]j[/math]-му объекту медоидом станет новый выбранный медоид ([math]i[/math]-ый объект), следовательно будет справедливо равенство [math]C_{ij} = D_{j} - d_{ij}[/math].

Вычисление значения [math]C_{ij}[/math]

Значение [math]T_{ijh}[/math] представляет собой число, которое добавляется к значению [math]D_{j}[/math] в результате выполнения операции обмена [math]i[/math]-го медоида и [math]h[/math]-го объекта. Для вычисления [math]C_{ij}[/math] следует рассмотреть два случая. В первом случае, когда [math]D_{j} = d_{ij}[/math] ([math]i[/math]-й объект является медоидом [math]j[/math]-го объекта) может быть два варианта: [math]E_{j} \lt = d_{jh}[/math] (медоидом [math]j[/math]-го объекта станет второго ближайший к нему медоид), тогда справедливо равенство [math]T_{ijh} = E_{j} - D_{j}[/math] и [math]E_{j} \gt d_{jh}[/math] (медоидом [math]j[/math]-го объекта станет [math]h[/math]-й объект), тогда справедливо равенство [math]T_{ijh} = d_{jh} - D_{j}[/math]. Во втором случае, когда [math]D_{j} != d_{ij}[/math] ([math]i[/math]-й объект НЕ является медоидом [math]j[/math]-го объекта) также может быть два варианта: [math]D_{j} \lt = d_{jh}[/math] (медоид [math]j[/math]-го объекта не изменится), тогда справедливо равенство [math]T_{ijh} = 0[/math] и [math]D_{j} \gt d_{jh}[/math] (медоидом [math]j[/math]-го объекта станет [math]h[/math]-й объект), тогда справедливо равенство [math]T_{ijh} = d_{jh} - D_{j}[/math]. Таким образом, в первом случае ([math]i[/math]-й объект является медоидом [math]j[/math]-го объекта) значение [math]T_{ijh}[/math] может быть вычислено по формуле [math]T_{ijh} = min(E_{j}, d_{jh}) - D_{j}[/math], а во втором ([math]i[/math]-й объект НЕ является медоидом [math]j[/math]-го объекта) – по формуле [math]T_{ijh} = min(d_{jh} - D_{j}, 0)[/math].

Вычисление значения [math]T_{ijh}[/math]

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

В соответствии с приведённым выше математическим описанием последовательное выполнение алгоритма представляет собой выполнение стадии BUILD, а затем последовательное повторение итераций стадии SWAP до тех пор пока существует пара объектов (один из которых медоид, а другой – нет), которые можно обменять, чтобы улучшить кластер.

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

При описании последовательной сложности будем считать, что выполняется распределение [math]N[/math] объектов по [math]K[/math] кластерам. В соответствии с приведённым математическим описанием алгоритм включает в себя подготовительную стадию BUILD и итерационную стадию SWAP.

Стадия BUILD сводится к вычислению:

Одной операции [math]arg\min_{i} \sum_{j} d_{ij}[/math]
[math]K-1[/math] операций [math]arg\max_{i \in U} \sum_{j \in U - \{i\}} C_{ij}[/math]

Т.е. стадия BUILD требует:

[math]K * N^2[/math] операций сложения
[math]K * N[/math] операций сравнения (определения максимума/минимума)
Около [math]K * N^2[/math] вычислений [math]C_{ij}[/math]

Одна итерация стадии SWAP сводится к вычислению:

[math]arg\min_{i \in S, h \in U} \sum_{j \in U - \{h\}} T_{ijh}[/math]

Т.е. одна итерации стадии SWAP требует порядка:

[math]K * N^2[/math] операций сложения
[math]K * N[/math] операций сравнения (определения максимума/минимума)
[math]K * (N - K) * (N - K - 1)[/math] вычислений [math]T_{ijh}[/math]

Вычисления [math]C_{ij}[/math] и [math]T_{ijh}[/math] имеют константную сложность для любых [math]i[/math], [math]j[/math] и [math]h[/math], поэтому сложность одной итерации алгоритма определяется как [math]O(K * N^2)[/math]. Таким образом можно отнести алгоритм PAM к алгоритмам с квадратичной сложностью по количеству [math]N[/math] объектов.

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

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

1. Матрица расстояния (в смысле непохожести) между объектами кластеризации (dissimilarity matrix) [math]D[/math] (элементы [math]d_{ij}[/math]).
Дополнительные ограничения:
   [math]D[/math] – симметрическая матрица порядка [math]N[/math], т.е. [math]d_{ij}= d_{ji}, i, j = 1, \ldots, N[/math]
   [math]D[/math] – содержит нули на главной диагонали, т.е. [math]d_{ii}= 0, i = 1, \ldots, N[/math]
2. [math]K[/math] - число кластеров, на которое следует разбить объекты кластеризации

Объём входных данных:

[math]\frac{N (N - 1)}{2}[/math] (в силу симметричности и нулевой главной диагонали достаточно хранить только  над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. 

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

[math]N[/math] чисел [math]k_{1}, k_{2}, \ldots, k_{N}[/math], где [math]k_{i}[/math] - целое число, соответствующее кластеру [math]i[/math]-го объекта.

Выходными данными алгоритма также можно считать [math]K[/math] чисел, задающих номера объектов кластеризации, выделенных в качестве медоидов (medoids). Однако в большинстве случаев требуется именно определение кластеров объектов, а не поиск соответствующих медоидов.

Объём выходных данных:

[math]N[/math] (кластеры объектов кластеризации)

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

Ограничения: небольшой объем данных. Достоинства: простота использования; быстрота использования; понятность и прозрачность алгоритма, алгоритм менее чувствителен к выбросам в сравнении с k-means. Недостатки: необходимо задавать количество кластеров; медленная работа на больших базах данных.

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

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

PAM работает эффективно для небольших наборов данных, но не очень хорошо масштабируется для больших наборов данных