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

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

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

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

1.1.1 EM-алгоритм

EM-алгоритм (англ. expectation-maximization) - алгоритм, используемый в математической статистике для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных. Каждая итерация алгоритма состоит из двух шагов. На E-шаге (expectation) вычисляется ожидаемое значение функции правдоподобия, при этом скрытые переменные рассматриваются как наблюдаемые. На M-шаге (maximization) вычисляется оценка максимального правдоподобия, таким образом увеличивается ожидаемое правдоподобие, вычисляемое на E-шаге. Затем это значение используется для E-шага на следующей итерации. Алгоритм выполняется до сходимости.

Как правило, ЕМ-алгоритм применяется для решения задач двух типов:

• К первому типу можно отнести задачи, связанные с анализом действительно неполных данных, когда некоторые статистические данные отсутствуют в силу каких-либо причин.
• Ко второму типу задач можно отнести те задачи, в которых функция правдоподобия имеет вид, не допускающий удобных аналитических методов исследования, но допускающий серьезные упрощения, если в задачу ввести дополнительные «ненаблюдаемые» (скрытые, латентные) переменные. Примерами прикладных задач второго типа являются задачи распознавания образов, реконструкции изображений. Математическую суть данных задач составляют задачи кластерного анализа, классификации и разделения смесей вероятностных распределений.

Алгоритм основан на методике итеративного вычисления оценок максимального правдоподобия, предложенной в 1977 г. (A. P. Demster, N. M. Laird, D. B. Rubin. Maximum Likelihood from Incomplete Data via the EM Algorithm).

1.1.2 Кластеризация с помощью EM-алгоритма

Кластеризация — задача разбиения заданной выборки объектов (ситуаций) на непересекающиеся подмножества, называемые кластерами, так, чтобы каждый кластер состоял из схожих объектов, а объекты разных кластеров существенно отличались. В основе идеи применения EM-алгоритма для кластеризации лежит предположение, что исследуемое множество данных может быть смоделировано с помощью линейной комбинации многомерных нормальных распределений, а целью является оценка параметров распределения, которые максимизируют логарифмическую функцию правдоподобия, используемую в качестве меры качества модели. К оцениваемым параметрам относятся математические ожидания и ковариационные матрицы. Предполагается, что любое наблюдение принадлежит ко всем кластерам, но с разной вероятностью. Задача заключается в "подгонке" распределений смеси к данным, а затем в определении вероятностей принадлежности наблюдения к каждому кластеру. Наблюдение должно быть отнесено к тому кластеру, для которого данная вероятность выше.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math] \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-алгоритм.

1.2.2 EM-алгоритм

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

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

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

     [math] \begin{align} & y_{ij} = \frac{w_j p_j(x_i)}{\sum_{l=1}^{k}w_lp_l(x_i)} \end{align} [/math]

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

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

Итерации происходят до сходимости (норма разности векторов скрытых переменных или изменение логарифмического правдоподобия на каждой итерации не будет превышать заданную константу) или достижения максимального числа итераций.

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

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

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

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

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

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

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

Шаг E

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

Шаг M

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

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

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

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

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

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

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

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

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

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

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

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

1.10.1 Достоинства и недостатки алгоритма

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

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

К недостаткам EM-алгоритма можно отнести следующее:

• Алгоритм неустойчив по начальным данным (то есть тем, которые инициализируют вектор параметров на первой итерации), так как он находит локальный экстремум, значение которого может оказаться гораздо ниже, чем глобальный максимум. В зависимости от выбора начального приближения алгоритм может сходиться к разным точкам. Также может сильно варьироваться скорость сходимости.
• Алгоритм не позволяет определять количество [math]k[/math] компонент смеси. Эта величина является структурным параметром алгоритма.
• Предположение о нормальности всех измерений данных не всегда выполняется.

3 Литература

1. A. Dempster, N. Laird and D. Rubin. Maximum likelihood estimation from incomplete data. – Journal of the Royal Statistical Society, Series B, 1977, vol. 39, p. 1-38.