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

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

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

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

Самый распространённый вариант реализации $k$-medoids называется PAM (Partitioning Around Medoids). Он представляет собой «жадный» алгоритм с очень неэффективной эвристикой. Данный метод был предложен Кауфманом в 1987 году.

Алгоритм начинается с выбора набора данных, состоящих из медоидов и остальных объектов. После выбора $k$ произвольных объектов в качестве медоидов, необходимо для каждой пары ${x_i}$ и ${m_k}$ вычислить значение $S({x_{i}}{m_{k}})$. По ходу выполнения алгоритма происходит итеративная замена одного из медоидов $x_i$ одним из остальных объектов $m_k$, если посчитанное ранее значение $S$ меньше нуля. Процесс будет повторятся, пока медоиды не стабилизируются.

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

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

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

Для заданных входных данных, алгоритм выбирает $K$ объектов - медоидов (medoids) из $N$ объектов кластеризации. При этом выбор медоидов осуществляется исходя из минимизации средней непохожести (average dissimilarity) объектов и ближайшего к ним медоида. От минимизации средней непохожести можно перейти к минимизации суммы непохожестей объектов на ближайший к ним медоид, т.к. это не повлияет на выбор медоидов, но позволит не накапливать ошибки округления, возникающие в результате операции деления.

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

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

Перед выполнением алгоритма множество $U$ содержит все объекты кластеризации, а множество $S$ – пусто.

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

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

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

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

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

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

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

В результате работы алгоритма множество $S$ будет содержать ровно $K$ объектов – медоидов, принадлежащих разным кластерам. Кластер объекта $j$ из множества $U$ совпадает с кластером ближайшего к нему медоида: $arg\min_{i \in S} d_{ij}$

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

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

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

Т.е. выполняется выбор такой пары медоида $i$ и объекта $h$ операция обмена которых приведёт к улучшению кластеризации или же делается вывод о невозможности такого обмена и работа алгоритма завершается.

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

• $D_{i}$ – непохожесть $i$-го объекта на ближайший к нему медоид
• $E_{i}$ – непохожесть $i$-го объекта на второй ближайший к нему медоид

Числа $D_{i}$ и $E_{i}$ могут измениться при изменениях в множествах $U$ и $S$.

2. Для каждого объекта $i \in U$ выполнить шаги 3-5

3. Для каждого объекта $j \in U - \{i\}$ вычислить $D_{j}$ и выполнить шаг 4

4. Положить $C_{ij} = max\{ D_{j} - d_{ij}, 0 \}$, мера улучшения кластеризации

5. Положить $g_{i} = \sum_{j \in U - \{i\}} C_{ij}$, улучшение от выбора $i$-го объекта в качестве медоида

6. Выбрать объект, максимизирующий величину $g_{i}$ в качестве медоида:

$\begin{align} k & = arg\min_{i \in U} \sum_{j \in U} d_{ij} \\ U & = U - \{k\} \\ S & = S + \{k\} \\ \end{align}$

7. Если множество $S$ содержит менее $K$ объектов (медоидов) перейти к шагу 2

8. КОНЕЦ стадии BUILD

Стадия SWAP (итерационное улучшение кластеризации):

1. Для всех возможных пар объектов $i \in S$ и $h \in U$ выполнить шаги 2-4

2. Вычислить эффект $T_{ih}$ от выполнения обмена объектов $i$ и $h$ (т.е. объект $i$ больше не выбран в качестве медоида, а объект $h$, наоборот, выбран в качестве медоида). Способ вычисления $T_{ih}$ описан ниже.

3. Положим $(i', h') = arg\min_{i \in S, h \in U} T_{ih}$

4. Если $T_{i'h'}$ больше 0 (т.е. кластеризацию можно улучшить), то выполняем обмен объектов $i$ и $h$ ($U = U - \{k\}, S = S + \{k\}$) и переходим к шагу 1

5. КОНЕЦ стадии SWAP

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

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

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

При описании последовательной сложности будем считать, что выполняется распределение $N$ объектов по $K$ кластерам. В соответствии с приведённым математическим описанием алгоритм включает в себя подготовительную стадию 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]

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

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] (кластеры объектов кластеризации)

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

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

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