Участник:Sergey Lavrushkin/EM-алгоритм кластеризации

Автор статьи: Сергей Лаврушкин (группа 620)

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

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

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-алгоритма лежит предположение, что исследуемое множество данных может быть смоделировано с помощью линейной комбинации многомерных нормальных распределений, а целью является оценка параметров распределения, которые максимизируют логарифмическую функцию правдоподобия, используемую в качестве меры качества модели. Иными словами, предполагается, что данные в каждом кластере подчиняются определенному закону распределения, а именно, нормальному распределению. С учетом этого предположения можно определить параметры - математическое ожидание и дисперсию, которые соответствуют закону распределения элементов в кластере, наилучшим образом "подходящему" к наблюдаемым данным.

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

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

Преимущества и недостатки алгоритма

Среди преимуществ EM-алгоритма можно выделить следующие:

• Мощная статистическая основа.
• Линейное увеличение сложности при росте объема данных.
• Устойчивость к шумам и пропускам в данных.
• Возможность построения желаемого числа кластеров.
• Быстрая сходимость при удачной инициализации.

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

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

Исходные данные:

• [math]\displaystyle{ k }[/math] — число кластеров,
• [math]\displaystyle{ X =\{x_{1}, x_{2}, ..., x_{n}\} }[/math] — множество из [math]\displaystyle{ n }[/math] наблюдений [math]\displaystyle{ q }[/math]-мерного пространства,
• [math]\displaystyle{ \epsilon }[/math] — допустимое отклонение для логарифмического правдоподобия,
• [math]\displaystyle{ m }[/math] — максимальное число итераций

Вычисляемые данные:

• [math]\displaystyle{ (w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k) }[/math] — параметры смеси гауссовых распределений
• [math]\displaystyle{ Y }[/math] — матрица с вероятностями членства в кластерах

Постановка задачи разделения смеси гауссовых распределений

Пусть [math]\displaystyle{ w_1, ..., w_k }[/math] — априорные вероятности кластеров, [math]\displaystyle{ p_1(x), ..., p_k(x) }[/math] — плотности распределения кластеров, тогда плотность распределения вектора признаков [math]\displaystyle{ x }[/math] сразу по всем кластерам равна:

[math]\displaystyle{ \begin{align} p(x) = \sum_{j=1}^{k}w_jp_j(x) \end{align} }[/math]

Необходимо на основе выборки оценить параметры модели [math]\displaystyle{ w_1, ..., w_k, p_1(x), ..., p_k(x) }[/math]. Это позволит оценивать вероятность принадлежности к кластеру и, таким образом, решить задачу кластеризации. Такая задача называется задачей разделения смеси распределений:

[math]\displaystyle{ \begin{align} p(x) = \sum_{j=1}^{k}w_jp_j(x),\quad p_j(x) = \phi(\theta_j; x), \end{align} }[/math]

где [math]\displaystyle{ \theta_j }[/math] — параметры распределения [math]\displaystyle{ p_j(x) }[/math]. В случае смести гауссовых распределений предполагается, что все компоненты имеют многомерное нормальное распределение. То есть в случае [math]\displaystyle{ q }[/math]-мерного пространства признаков [math]\displaystyle{ p_j(x) }[/math] имеют следующий вид:

[math]\displaystyle{ \begin{align} & p_j(x) = \frac{1}{\Big(2\pi\Big)^{\frac{q}{2}}\sqrt{|\Sigma_j|}}\exp{\bigg\{-\frac{\delta^2}{2}\bigg\} },\\ & \delta^2 = \Big( x - \mu_j \Big)^T \Sigma_j^{-1} \Big( x - \mu_j \Big), \end{align} }[/math]

где:

• [math]\displaystyle{ \Sigma_j }[/math] — ковариационная матрица размером </math>q \times q</math>,
• [math]\displaystyle{ \mu_j }[/math] — [math]\displaystyle{ q }[/math]-мерный вектор математических ожиданий,
• [math]\displaystyle{ \delta^2 }[/math] — квадратичное расстояние Махаланобиса.

В случае смеси гауссовых распределений получаем, что [math]\displaystyle{ \theta_j = (\mu_j, \Sigma_j) }[/math]. То есть для решения задачи разделения смеси гауссовых распределений необходимо оценить вектор параметров [math]\displaystyle{ (w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k) }[/math]. Согласно принципу максимизации правдоподобия:

[math]\displaystyle{ \begin{align} w, \theta = \underset{w, \theta}{\operatorname{argmax}} \sum_{i=1}^{n}\ln{p(x_i)} = \underset{w, \theta}{\operatorname{argmax}} \sum_{i=1}^{n}\ln{\sum_{j=1}^{k}w_jp_j(x_i)} \end{align} }[/math]

Таким образом, имеет место задача максимизации суммы логарифмов сумм, решение которой представляет большую трудность. В таком случае полезным оказывается итеративный метод решения — EM-алгоритм.

EM-алгоритм

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

EM-алгоритм заключается в следующем:

1. В начале работы алгоритма задаются параметры начального приближения. Наиболее общим способом инициализации является присвоение элементам матрицы математических ожиданий случайных значений [math]\displaystyle{ \mu_j \leftarrow Random }[/math], начальные ковариационные матрицы определяются как единчиные [math]\displaystyle{ \Sigma_j \leftarrow I }[/math], веса кластеров задаются одинаковыми [math]\displaystyle{ w_i \leftarrow \frac{1}{k} }[/math]. Также в качестве начальных параметров можно использовать результат работы алгоритма K-means. (Данная эвристика применяется, так как K-means требуется намного меньше итераций до достижения стабилизации, в то время как каждый шаг EM требует больших вычислительных затрат).
2. Далее итеративно выполняется следующая пара процедур:
   • E-шаг: используя текущее значение вектора параметров [math]\displaystyle{ (w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k) }[/math], вычисляем значение вектора скрытых переменных [math]\displaystyle{ g }[/math]:

     [math]\displaystyle{ \begin{align} & g_{ij} = \frac{w_j p_j(x_i)}{\sum_{l=1}^{k}w_lp_l(x_i)}, \\ & n_j = \sum_{i = 1}^{n}g_{ij}. \end{align} }[/math]

   • М-шаг: переоценка вектора параметров, используя текущее значение вектора скрытых переменных:

     [math]\displaystyle{ \begin{align} & \mu_j^{new} = \frac{1}{n_j} \sum_{i = 1}^{n}g_{ij}x_i, \\ & \Sigma_j^{new} = \frac{1}{n_j} \sum_{i = 1}^{n}g_{ij}(x_i - \mu_j^{new})(x_i - \mu_j^{new})^T,\\ & w_j^{new} = \frac{n_j}{n}. \end{align} }[/math]

Итерации происходят до сходимости или достижения максимального числа итераций.

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

Важнейшим действием вычислительного ядра EM-алгоритма, выполняемым на E-шаге, является вычисление скрытых переменных [math]\displaystyle{ g_{ij} }[/math] (всего их [math]\displaystyle{ nk }[/math]). Для M-шага таким действием является вычисление [math]\displaystyle{ g_{ij}(x_i - \mu_j^{new})(x_i - \mu_j^{new})^T }[/math] (всего их [math]\displaystyle{ nk }[/math]).

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

Как записано в описании ядра алгоритма, основную часть метода составляют множественные вычисления скрытых переменных [math]\displaystyle{ g_{ij} }[/math] на E-шаге и выражений вида [math]\displaystyle{ g_{ij}(x_i - \mu_j^{new})(x_i - \mu_j^{new})^T }[/math] на M-шаге. Сложность вычисления этих выражений в случае недиагональных матриц ковариаций — [math]\displaystyle{ O(q^2) }[/math], а в случае диагональных — [math]\displaystyle{ O(q) }[/math].

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

Последовательная реализация алгоритма EM может быть проиллюстрирована с помощью следующего псевдокода:

1. Инициализация: установка начальных значений [math]\displaystyle{ \mu, \, \Sigma_j \, (j = \overline{1, k}), \, w }[/math].
2. Пока изменение логарифмического правдоподобия [math]\displaystyle{ \Delta llh \geq \epsilon }[/math] и не достигнуто максимальное число итераций [math]\displaystyle{ m }[/math], выполнять шаги E и M.

Шаг E

def E-step:
    Инициализировать нулевыми значениями [math]\displaystyle{ \mu^{'}, \, \Sigma_j^{'} \, (j = \overline{1, k}), \, w^{'}, \, llh }[/math]
    for [math]\displaystyle{ j = \overline{1, k} }[/math]:
        Вычислить определитель и обратную матрицу матрицы [math]\displaystyle{ \Sigma_j }[/math]
    for [math]\displaystyle{ i = \overline{1, n} }[/math]:
        [math]\displaystyle{ sump_i = 0 }[/math]
        for [math]\displaystyle{ j = \overline{1, k} }[/math]:
            [math]\displaystyle{ \delta_{ij} = \big( x_i - \mu_j \big)^T \Sigma_j^{-1} \big( x_i - \mu_j \big) }[/math]
            [math]\displaystyle{ gij = \frac{w_j}{{2 \pi}^\frac{q}{2}|\Sigma_j|}\exp{\big\{-\frac{\delta^2}{2}\big\} } }[/math]
            [math]\displaystyle{ sump_i \mathrel{+}= g_{ij} }[/math]
        [math]\displaystyle{ y_i = \frac{g_i}{sump_i} }[/math]
        [math]\displaystyle{ llh \mathrel{+}= \ln(sump_i) }[/math]
        [math]\displaystyle{ \mu^{'} \mathrel{+}= x_iy_i^T }[/math]
        [math]\displaystyle{ w^{'} \mathrel{+}= y_i }[/math]

Шаг M

def M-step:
    for [math]\displaystyle{ j = \overline{1, k} }[/math]:
        [math]\displaystyle{ \mu_j = \frac{\mu_j^{'}}{w_j^{'}} }[/math]
        for [math]\displaystyle{ i = \overline{1, n} }[/math]:
            [math]\displaystyle{ \Sigma_j^{'} \mathrel{+}= (x_i - \mu_j)y_{ij}(x_i - \mu_j)^T }[/math]
        [math]\displaystyle{ \Sigma_j = \frac{\Sigma_j^{'}}{w_j^{'}} }[/math]
    [math]\displaystyle{ w = \frac{w^{'}}{n} }[/math]

Логаририфмическое правдоподобие вычисляется как:

[math]\displaystyle{ \begin{align} llh=\sum_{i=1}^{n}\ln{(sump_i)} \end{align} }[/math]

Матрицы [math]\displaystyle{ \mu^{'}, \, \Sigma_j^{'} \, (j = \overline{1, k}), \, w^{'} }[/math] являются временными и используются только для вычислений.

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

В случае диагональности матриц ковариации определитель матрицы и ее обращение может быть вычислено за время [math]\displaystyle{ O(q) }[/math], а алгоритм имеет сложность [math]\displaystyle{ O(mkqn) }[/math]. В случае недиагональной матрицы сложность алгоритма составит [math]\displaystyle{ O(mkq^2n) }[/math], то есть будет квадратично возрастать с увеличением размерности данных.

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

Входные данные: [math]\displaystyle{ X =\{x_{1}, x_{2}, ..., x_{n}\} }[/math] — множество из [math]\displaystyle{ n }[/math] наблюдений [math]\displaystyle{ q }[/math]-мерного пространства.

Объем входных данных: [math]\displaystyle{ nq }[/math]

Выходные данные: [math]\displaystyle{ (w, \theta) = (w_1, ..., w_k; \mu_1, ..., \mu_k; \Sigma_1, ..., \Sigma_k) }[/math] — параметры смеси гауссовых распределений; [math]\displaystyle{ Y }[/math] — матрица с вероятностями членства в кластерах.

Объем выходных данных: [math]\displaystyle{ k(q^2 + q + 1) + nq }[/math]