Участник:Александр Куваев/Алгоритм кластеризации, основанный на максимизации ожидания

Авторы описания: Куваев А.С. и Щенявская Е.В.

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

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

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

Решение этой задачи принципиально неоднозначно по следующим причинам:

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

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

$p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \ \sum_{j=1}^{k}w_{j}=1, \ w_{j} \ge 0$, где $p_{j}$ - плотность распределения $j$-й компоненты смеси (кластера).

Везде в дальнейшем будем предполагать, что $p_{j}$ имеют вид многомерных нормальных плотностей с произвольной матрицей ковариации: смеси нормальных распределений позволяют аппроксимировать произвольные непрерывные функции плотности с наперед заданной точностью. Результатом работы EM-алгоритма являются оценки априорных вероятностей компонент смеси $w_{j}$, а также оценки векторов математических ожиданий и матриц ковариации для каждой компоненты. Зная параметры распределения смеси, каждому наблюдению будет сопоставляться тот кластер, вероятность принадлежности к которому будет максимальной.

EM-алгоритм позволяет значительно упростить задачу максимизации правдоподобия выборки путем искусственного введения вспомогательной матрицы скрытых переменных $G$. Алгоритм заключается в последовательном повторении шагов E(expectation) и M(maximization):

• На шаге E на основе текущего приближения параметров смеси по формуле Байеса вычисляются ожидаемые значения скрытых переменных $g_{ij}$ — апостериорные вероятности того, что $i$-й объект принадлежит кластеру $j$.
• На шаге M решается задача максимизации правдоподобия для нахождения следующего приближения параметров смеси на основе текущего приближения и матрицы скрытых переменных, при этом решение этой задачи выписывается в явном виде.

Рисунок 1. Визуализация компонент смеси после некоторых итераций алгоритма в двумерном пространстве признаков. Для каждой компоненты эллипсами выделены доверительные области с уровнями доверия 0.9, 0.95 и 0.99.

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

Пусть заданы $l$ объектов $x_{1},\dotsc,x_{l}$, каждый из которых описывается $n$ числовыми признаками. Таким образом, определена матрица объектов-признаков $X \in \R^{l \times n}.$ Предполагается, что объекты выбраны случайно и независимо из смеси $n$-мерных нормальных распределений с фиксированным числом компонент $k$. Плотность распределения на множестве объектов имеет вид:

$p(x)=\sum_{j=1}^{k}w_{j}p_{j}(x), \ \sum_{j=1}^{k}w_{j}=1, \ w_{j} \ge 0$, где $p_{j}$ - плотность распределения $j$-й компоненты смеси.

Каждая компонента смеси описывается $n$-мерным вектором математических ожиданий $\mu_{j}$ и матрицей ковариации $\Sigma_{j}$ порядка $n$, $j = 1,\dotsc,k$.

Исходные данные: матрица объектов-признаков $X$, число кластеров $k$, максимальное число итераций $imax$, минимальная величина изменения логарифма правдоподобия $\varepsilon$.

Вычисляемые данные: набор оценок весов, математических ожиданий и ковариационных матриц компонент смеси $\theta = (w_{1},...,w_{k}; \; \mu_{1},...,\mu_{k}; \; \Sigma_{1},...,\Sigma_{k})$, максимизирующий правдоподобие выборки.

Структура алгоритма:

• Инициализация параметров: существует большое число вариантов инициальзации параметров распределения. Один из возможных подходов — задание векторов математических ожиданий компонент случайными элементами выборки, задание матриц ковариации единичными матрицами и задание весов компонент равными $\frac{1}{k}.$
• Последовательное выполнение шагов E и M до тех пор, пока правдоподобие выборки не стабилизируется или не будет достигнуто максимальное число итераций:
  • Шаг E: вычисление значений скрытых переменных по формуле Байеса:

  $g_{ij} = \frac{w_{j} p_{j}(x_{i})}{\sum_{s=1}^{k} w_{s} p_{s}(x_{i})}$ , где $p_{j}(x_{i}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma_{j}|}} \exp \biggl( -\frac{1}{2}(x_{i} - \mu_{j})^{T} \Sigma_{j}^{-1} (x_{i} - \mu_{j}) \biggr), \ i = 1,\dotsc,l; \ j = 1,\dotsc,k$

  • Шаг M: перерасчет параметров смеси на основе текущего приближения и матрицы скрытых переменных:

  $w_{j} = \frac{1}{l} \sum_{i=1}^{l} g_{ij}, \ j = 1,\dotsc,k;$

  $\mu_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij} x_{i}, \ j = 1,\dotsc,k;$

  $\Sigma_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij}(x_{i} - \mu_{j})(x_{i} - \mu_{j})^T, \ j = 1,\dotsc,k.$

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

Вычислительное ядро EM-алгоритма — процедура последовательного выполнения шагов E и M:

• На шаге E на основе текущего приближения параметров смеси вычисляются ожидаемые значения скрытых переменных
• На шаге M вычисляется следующее приближение параметров смеси на основе текущего приближения и матрицы скрытых переменных

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

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

Алгоритма состоит из двух итерационно проводимых макроопераций: E-шага и M-шага.

В ходе E-шага используется:

• Обращение ковариационной матрицы и вычисление ее определителя
• Умножение вектор-строки на матрицу и вектор-строки на вектор-столбец

В ходе M-шага используется:

• Умножение вектор-столбца на вектор-строку
• Взвешенное суммирование векторов и матриц

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

1. Для всех $j = 1,\dotsc,k$:

   Инициализировать $w_{j}, \ \mu_{j}, \ \Sigma_{j}$

2. Для всех $iter = 1,\dotsc,imax$:
   1. Шаг E:

      Для всех $j = 1,\dotsc,k$:

      Вычислить $|\Sigma_{j}|, \ \Sigma_{j}^{-1}$

      Для всех $i = 1,\dotsc,l$:

      $sum = 0$

      Для всех $j = 1,\dotsc,k$:

      Вычислить $p_{j}(x_{i}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma_{j}|}} \exp \biggl( -\frac{1}{2}(x_{i} - \mu_{j})^{T} \Sigma_{j}^{-1} (x_{i} - \mu_{j}) \biggr)$

      $sum = sum + w_{j}p_{j}(x_{i})$

      Для всех $j = 1,\dotsc,k$:

      Вычислить $g_{ij} = \frac{w_{j} p_{j}(x_{i})}{sum}$

   2. Шаг M:

      Для всех $j = 1,\dotsc,k$:

      Вычислить $w_{j} = \frac{1}{l} \sum_{i=1}^{l} g_{ij}$

      Вычислить $\mu_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij} x_{i}$

      Вычислить $\Sigma_{j} = \frac{1}{l w_{j}} \sum_{i=1}^{l} g_{ij}(x_{i} - \mu_{j})(x_{i} - \mu_{j})^T$

   3. Вычислить изменения логарифма правдоподобия $\Delta$
   4. Если $\Delta \lt \varepsilon$, то досрочно выйти из цикла
3. Вернуть $w_{j}, \mu_{j}, \Sigma_{j}, \ j = 1,\dotsc,k.$

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

Рассмотрим мультипликативную сложность одной итерации алгоритма:

• E-шаг:
  • Сложность обращения матрицы ковариаций и вычисления ее определителя — $O(n^{3})$ (для простоты будем рассматривать метод Гаусса)
  • Сложность вычисления расстояния Махаланобиса (показателя экспоненты) при вычислении значения каждой скрытой переменной — $O(n^{2})$
  • Общая сложность E-шага — $O(k * n^{3} + k * l * n^{2})$
• M-шаг:
  • Сложность пересчета веса кластера — $O(l)$
  • Сложность пересчета центра кластера — $O(l * n)$
  • Сложность пересчета матрицы ковариации кластера — $O(l * n)$
  • Общая сложность M-шага — $O(k * l * n)$

Таким образом, общая сложность одной итерации — $O(k * n^{2} * (l + n))$. При учете ограничения на максимальное число итераций получим оценку общей сложности алгоритма — $O(imax * k * n^{2} * (l + n))$.

Большое влияние на сложность E-шага оказывает необходимость обращать ковариационные матрицы и вычислять их определители. Помимо того, что это трудоемкая операция, ковариационные матрицы могут оказаться плохо обусловленными, что может привести к неустойчивости выборочных оценок параметров смеси. Обращения матриц можно избежать, если принять гипотезу о том, что в каждой компоненте смеси признаки некоррелированы, то есть ковариационные матрицы диагональные. В таком случае общая сложность алгоритма составит $O(imax * k * n * l)$.

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

Рисунок 2. Граф шага E с отображением входных и выходных данных.

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

Входные данные: плотная матрица объектов-признаков $X \in \R^{l \times n}$, число кластеров $k \in \N$, максимальное число итераций $imax \in \N$, минимальная величина изменения логарифма правдоподобия $\varepsilon \in \R$.

Объём входных данных: $ln + 3$

Выходные данные: вещественный вектор весов $w \in \R^{k}$, $k$ вещественных векторов математических ожиданий $\mu_{j} \in \R^{n}$ и $k$ вещественных ковариационных матриц $\Sigma_{j} \in \R^{n \times n}$.

Объём выходных данных: $k(n^{2}+n+1)$