Введение
Кластеризация является одной из наиболее важных задач Data Mining. В настоящее время разработано большое количество методов и алгоритмов кластеризации, но, к сожалению, не все они могут эффективно работать с большими массивами данных, поэтому дальнейшие исследования в этом направлении связаны с преодолением этой проблемы. Одним из широко известных в аналитическом сообществе алгоритмов кластеризации, позволяющих эффективно работать с большими объемами данных, является EM-алгоритм. Его название происходит от слов "expectation-maximization", что переводится как "ожидание-максимизация". Это связано с тем, что каждая итерация содержит два шага – вычисление математических ожиданий (expectation) и максимизацию (maximisation). Алгоритм основан на методике итеративного вычисления оценок максимального правдоподобия, предложенной в 1977 г. (A. P. Demster, N. M. Laird, D. B. Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm).
В основе идеи EM-алгоритма лежит предположение, что исследуемое множество данных может быть смоделировано с помощью линейной комбинации многомерных нормальных распределений, а целью является оценка параметров распределения, которые максимизируют логарифмическую функцию правдоподобия, используемую в качестве меры качества модели. Иными словами, предполагается, что данные в каждом кластере подчиняются определенному закону распределения, а именно, нормальному распределению (рисунок 1). С учетом этого предположения можно определить параметры - математическое ожидание и дисперсию, которые соответствуют закону распределения элементов в кластере, наилучшим образом "подходящему" к наблюдаемым данным.
Таким образом, мы предполагаем, что любое наблюдение принадлежит ко всем кластерам, но с разной вероятностью. Тогда задача будет заключаться в "подгонке" распределений смеси к данным, а затем в определении вероятностей принадлежности наблюдения к каждому кластеру. Очевидно, что наблюдение должно быть отнесено к тому кластеру, для которого данная вероятность выше.
Среди преимуществ EM-алгоритма можно выделить следующие:
- Мощная статистическая основа.
- Линейное увеличение сложности при росте объема данных.
- Устойчивость к шумам и пропускам в данных.
- Возможность построения желаемого числа кластеров.
- Быстрая сходимость при удачной инициализации.
Однако алгоритм имеет и ряд недостатков. Во-первых, предположение о нормальности всех измерений данных не всегда выполняется. Во-вторых, при неудачной инициализации сходимость алгоритма может оказаться медленной. Кроме этого, алгоритм может остановиться в локальном минимуме и дать квазиоптимальное решение.
Статистические основы алгоритма
Как отмечалось во введении, EM-алгоритм предполагает, что кластеризуемые данные подчиняются линейной комбинации (смеси) нормальных (гауссовых) распределений. Плотность вероятности нормального распределения имеет вид:
$$p(x) = \frac {1}{\sqrt {2\times\pi\times\sigma^2}}\times exp\left\{-\frac {{(x-\mu)}^2}{2\times\sigma^2}\right\}$$
где
$\mu = E(X)$ - математическое ожидание,
$\sigma^2 = E{(X-\mu)}^2$ - дисперсия.
Многомерное нормальное распределение для $q$-мерного пространства является обобщением предыдущего выражения. Многомерная нормальная плотность для $q$-мерного вектора $x=\biggl (x_1\mathbf{,}$ $ x_2\mathbf{,}$ $\ldots\mathbf{,}$ $x_q \biggr)$ может быть записана в виде:
$$p(x) = \frac {1}{{\Bigl(2\times\pi\Bigr)}^{\frac{q}{2}}\times\sqrt{\mid\sum\mid}}\times exp \biggl\{ -\frac{1}{2}\times{\Bigl(x-\mu\Bigr)}^T\times{\sum}^{-1}\times\Bigl(x-\mu\Bigr)\biggr\}$$
где
$\sum$ - ковариационная матрица размером $q \times q$, которая, как известно, является обобщением дисперсии для многомерной случайной величины,
$\mu$ представляет из себя $q$-мерный вектор математических ожиданий,
$\biggl | \sum \biggr |$ - определитель ковариационной матрицы,
$T$- оператор транспонирования.
Введем в рассмотрение функцию $\delta ^2 = {\Bigl(x-\mu\Bigr)}^T \times{\sum}^{-1}\times\Bigl(x-\mu\Bigr)$, называемую квадратичным расстоянием Махаланобиса.
Алгоритм предполагает, что данные подчиняются смеси многомерных нормальных распределений для $q$ переменных. Модель, представляющая собой смесь гауссовых распределений задается в виде:
$$p \Bigl(x\Bigr)=\sum_{t=1}^k w_i\times {p}\Bigl(x\mid i\Bigr)$$
где
$p\Bigl(x\mid i\Bigr)$ - нормальное распределение для $i$-го кластера,
$w_i\space$ - доля (вес) $i$-го кластера в исходной базе данных.
Существуют два подхода к решению задач кластеризации: основанный на расстоянии и основанный на плотности. Первый подход заключается в определении областей пространства признаков, внутри которых точки данных расположены ближе друг к другу, чем к точкам других областей, относительно некоторой функции расстояния (например, евклидовой). Второй - обнаруживает области, которые являются более "заселенными", чем другие. Алгоритмы кластеризации могут работать сверху вниз (иерархические) и снизу вверх (агломеративные). Агломаративные алгоритмы, как правило, являются более точными, хотя и работают медленнее.
Алгоритм EM основан на вычислении расстояний. Он может рассматриваться как обобщение кластеризации на основе анализа смеси вероятностных распределений. В процессе работы алгоритма происходит итеративное улучшение решения, а остановка осуществляется в момент, когда достигается требуемый уровень точности модели. Мерой в данном случае является монотонно увеличивающаяся статистическая величина, называемая логарифмическим правдоподобием. Целью алгоритма является оценка средних значений $C$, ковариаций $R$ и весов смеси $W$ для функции распределения вероятности, описанной выше. Параметры, оцененные алгоритмом, сохраняются в таблице вида:
Матрица | Размер | Содержит |
---|---|---|
$C$ |
$q\times k$ |
Математические ожидания, $\mu$ |
$R$ |
$q\times q$ |
Ковариации, $\sum$ |
$W$ |
$k\times 1$ |
Веса, $w_i$ |
Следует отметить, что один из популярных алгоритмов кластеризации k-means является частным случаем алгоритма EM, когда $W$ и $R$ постоянны: $w_ i = \frac {1}{k}$, $R=I$ ($I$ - единичная матрица).
Алгоритм начинает работу с инициализации, т.е. некоторого приближенного решения, которое может быть выбрано случайно или задано пользователем исходя из некоторых априорных сведений об исходных данных. Наиболее общим способом инициализации является присвоение элементам матрицы математических ожиданий случайных значений $C\leftarrow \mu \space$ $ \mbox {Random}$, начальная ковариационная матрица определяется как единчиная $r\leftarrow I$, веса кластеров задаются одинаковыми ($w_i\leftarrow \frac{1}{k}$).
Следует обратить внимание, что алгоритм может "застрять" в локальном оптимуме и дать квазиоптимальное решение при выборе неудачного начального приближения. Поэтому одним из его недостатков следует считать чувствительность к выбору начального состояния модели.
Реализация алгоритма EM может быть проиллюстрирована с помощью следующего псевдокода:
Вход: $k$ – число кластеров,
$Y = \biggl\{y_1,$ $y_2,$ $ \dots,y_n \biggr\}$– множество из $n$ наблюдений $q$-мерного пространства,
$\varepsilon$ - допустимое отклонение для логарифмического правдоподобия,
$Q$ - максимальное число итераций,
Выход: $C$, $R$, $W$ - матрицы, содержащие обновляемые параметры смеси.
$X$ - матрица с вероятностями членства в кластерах.
- Инициализация: установка начальных значений $C$, $R$, $W$, выбранных случайно или заданных пользователем.
- ПОКА изменение логарифмического правдоподобия $\Delta {llh}\geq\varepsilon$ и не достигнуто максимальное число итераций $Q$, выполнять шаги E и M.
Шаг E
$C^,=0$, $R^,=0$, $W^,=0$, $llh=0$
Для $i$, изменяющегося от $1$ до $n$
$sump_i=0$
Для $j$, изменяющегося от $1$ до $k$
$\delta_{ij}={(y_i - C_j)}^T\times R_j^{-1}\times(y_i - C_j)$
$p_{ij} =\frac{w_j}{(2\times\pi)^{\frac{q}{2}}\times |R_j|^\frac{1}{2}}\times exp \biggl\{-\frac{1}{2}\times\delta_{ij}\biggr\}$
$sump_i = sump_i$ $+$ $p_{ij}$
Конец цикла по $j$
$x_i = \frac{p_i}{sump_i}$, $llh = llh$ $+$ $\ln \biggl(sump_i\biggr)$
$C^, = C^, + y_i$ $ {x_i}^T$, $W^, = W^, + x_i$
Конец цикла по $i$
Шаг M
Для $j$, изменяющегося от $1$ до $k$
$C_j = \frac{{C_j}^,}{{W_j}^,}$
Для $i$, изменяющегося от $1$ до $n$
$R_j^, = R_j^, +(y_i-C_j)x_{ij}{(y_i - C_j)}^T$
Конец цикла по $i$
$R_j = \frac{R_j^,}{n}$, $W = \frac{W^,}{n}$
Конец цикла по $j$
Алгоритм содержит два шага: шаг ожидания (expectation) или Е-шаг и шаг максимизации (maximization) или M-шаг. Каждый из них повторяется до тех пор, пока изменение логарифмического правдоподобия $\Delta llh$ не станет меньше, чем $\varepsilon$, или пока не будет достигнуто максимальное число итераций.
Логарифмическое правдоподобие вычисляется как:
$$llh = \sum_{i=1}^n\ln\Bigl(sump_i\Bigr)$$
Переменные $\delta$, $R$, $P$ представляют собой матрицы, хранящие расстояния Махаланобиса, ковариации и вероятности членства в кластере для каждой из $n$ точек. $C^,$, $R^,$ и $W^,$ являются временными матрицами, используемыми только для вычислений. $\mid\mid W \mid\mid \,= 1$, т.е. $\sum \limits_{i=1}^{k} w_i = 1$. Обозначение вида $p_i$, использованное в псевдокоде, обозначает k-размерный вектор принадлежности $i$-го наблюдения к каждому из $k$ кластеров. Соответственно, $x_i$ – нормированная вероятность принадлежности к каждому из k кластеров. Столбец $C_j$ матрицы $C$ матрицы есть оценка математического ожидания по $j$-му кластеру, $R$ - диагональная матрица, т.е. $R_{ij} = 0$ для всех $i\neq j$. Со статистической точки зрения это означает, что ковариации являются независимыми.
Диагональность является ключевым предположением, которое делает алгоритм масштабируемым. В этом случае детерминант матрицы и его обращение может быть вычислено за время $O\bigl(p\bigr)$, а алгоритм имеет сложность $O\bigl(k\,p\,n\bigr)$. В случае недиагональной матрицы сложность алгоритма составит $O\bigl(k\,p^2n\bigr)$, т.е. будет квадратично возрастать с увеличением размерности данных.
Важнейшим действием, выполняемым на E-шаге, является вычисление расстояний Махаланобиса $\delta_{ij}$. Если матрица $R$ является диагональной, то расстояние Махаланобиса от точки $y$ до среднего значения кластера $C$, имеющего ковариацию $R$, будет:
$$\delta _{ij}={\bigl( y\,-\,C\bigr)}^TR^{-1}\bigl( y\,-\,C\bigr) = \sum_{k=1}^q \frac{{\bigl( y_k\,-\,C_k\bigr)}^2}{R_{kk}}$$
,поскольку ${R_{kk}}^{-1} = \frac{1}{R_{kk}}$, $\overline{k=1,q}$. Если матрица $R$ является диагональной, то ее обращение $R^{-1}$ легко вычисляется, т.к. $R^{-1}$ для любых $k\neq l$. Кроме этого, ускорению вычислений способствует то, что диагональная матрица $R$ может храниться в виде вектора ее диагональных элементов. Поскольку $R$ не изменяется в процессе Е-шага, ее детерминант вычисляется только единожды, что делает вычисление вероятностей $p_{ij}$ более быстрым. На M-шаге диагональность матрицы $R$ также упрощает вычисления, поскольку недиагональные элементы матрицы $\Bigl(y_i\,-\,C_j\Bigr)\,x_{ij}\,{\Bigl(y_i\,-\,C_j\Bigr)}^T$ равны нулю.
Для оптимизации используемого объема памяти, алгоритм может работать в двух режимах. В первом загружается только часть доступных данных и на их основе предпринимается попытка построения модели. Если она увенчалась успехом, то алгоритм завершает работу, в противном случае загружается следующая порция данных и т.д., пока не будут получены приемлемые результаты. Во втором режиме загружаются сразу все имеющиеся данные. Как правило, последний вариант обеспечивает более точную подгонку модели, но предъявляет более жесткие требования к объему доступной оперативной памяти.
Численный эксперимент
Для иллюстрации работы алгоритма EM и его сравнения с k-means рассмотрим результаты численного эксперимента, для проведения которого была взята выборка, представленная на рисунке 2.
Обратите внимание, что исходный набор данных не является простым с точки зрения задачи кластеризации, поскольку имеется явное перекрытие кластеров (области 1 и 2). В области 1 перекрываются кластеры 1 и 2, а в области 2 кластеры 4 и 5. Кластеры 3 и 6 расположены обособленно и, как ожидается, будут легко разделимы.
Для алгоритма k-means особые трудности должны возникнуть в местах перекрытия кластеров. Данное предположение подтверждается результатами, представленным на рисунке 3.
В местах перекрытия кластеров наблюдается наибольшее число ошибок. В то же время обособленные кластеры 3 и 6 были распознаны алгоритмом k-means без ошибок. Как можно увидеть на рисунке 4, алгоритм EM весьма успешно выявил перекрывающиеся кластеры, хотя и почти не распознал кластер 6.
Таким образом, можно сделать вывод, что алгоритм k-means может иметь преимущество при работе с обособленными (неперекрывающимися) кластерами, но полностью проигрывает алгоритму EM при наличии их перекрытия.
- Королёв В.Ю. ЕМ-алгоритм, его модификации и их применение к задаче разделения смесей вероятностных распределений. Теоретический обзор. - М.: ИПИРАН, 2007.
- Dempster, A., Laird, N., and Rubin, D. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39(1):1–38.
- McLachlan, G. and Krishnan, T. The EM algorithm and extensions. Wiley series in probability and statistics. John Wiley & Sons. (1997). Paul S. Bradley, Usama M. Fayyad, Cory A. Reina, Scaling EM (Expectation-Maximization) Clustering to Large Databases, Microsoft Research Technical Report MSR-TR-98-35. Redmond, WA 98052, 1999