EM-алгоритм

НЕТ ВОЙНЕ

24 февраля 2022 года российское руководство во главе с Владимиром Путиным развязало агрессивную войну против Украины. В глазах всего мира это военное преступление совершено от лица всей страны, всех россиян. Будучи гражданами Российской Федерации, мы против своей воли оказались ответственными за нарушение международного права, военное вторжение и массовую гибель людей. Чудовищность совершенного преступления не оставляет возможности промолчать или ограничиться пассивным несогласием. Мы убеждены в абсолютной ценности человеческой жизни, в незыблемости прав и свобод личности. Режим Путина — угроза этим ценностям. Наша задача — обьединить все силы для сопротивления ей. Эту войну начали не россияне, а обезумевший диктатор. И наш гражданский долг — сделать всё, чтобы её остановить. Антивоенный комитет России

Распространяйте правду о текущих событиях, оберегайте от пропаганды своих друзей и близких. Изменение общественного восприятия войны - ключ к её завершению.

meduza.io, Популярная политика, Новая газета, zona.media, Майкл Наки.

Определение

Алгоритм EM (англ. expectation-maximization) — итеративный алгоритм поиска оценок максимума правдоподобия модели, в ситуации, когда она зависит от скрытых (ненаблюдаемых) переменных.

Алгоритм ищет параметры модели итеративно, каждая итерация состоит из двух шагов:

E (Expectation) шаг — поиск наиболее вероятных значений скрытых переменных.

M (Maximization) шаг — поиск наиболее вероятных значений параметров, для полученных на шаге E значений скрытых переменных.

EM алгоритм подходит для решения задач двух типов:

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

Основной алгоритм

Постановка задачи

Плотность распределения смеси имеет вид:
[math]p(x) = \sum\limits_{j=1}^k w_j p_j(x)[/math].
Где [math] \sum\limits_{j=1}^k w_j = 1; w_j \geq 0; p_j(x) = \phi(x;\theta_j)[/math] — функция правдоподобия [math]j[/math]-ой компонеты смеси, [math]w_j[/math] — априорная вероятность [math]j[/math]-ой компоненты смеси.

Перед нами стоит две задачи:

1. По заданной выборке [math]X^m[/math] случайных и независимых наблюдений полученных из смеси [math]p(x)[/math], числу [math]k[/math] и функции [math]\phi[/math], оценить вектор параметров [math]\Theta = (w_1,..,w_k,\theta_1,..,\theta_k)[/math].
2. Найти [math]k[/math].

Проблема

Задачи подобного рода мы умеем решать, максимизируя логармиф правдоподобия:
[math] Q(\Theta) = ln \prod\limits_{i=1}^mp(x_i) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow \max\limits_{\Theta}[/math].

Но проблeма в том, что мы не знаем как аналитически посчитать логарифм суммы. Тут нам и поможет алгоритм EM.

Решение

Основная идея алгоритма EM заключается в том, что мы добавляем скрытые переменные такие, что:

1. Они могут быть выражены через [math]\Theta[/math].
2. Они помогают разбить сумму так: [math]p (X, H|\Theta) = \prod\limits_{i=1}^k p (X|H, \Theta) p(H|\Theta)[/math], где [math]H[/math] — матрица скрытых переменных.

Тогда алгоритм EM сводится к повторению шагов, указанных в Определении.

E-шаг

[math]p(x,\theta_j) = p(x)P(\theta_j | x) = w_jp_j(x)[/math].

Скрытые переменные представляют из себя матрицу [math]H = (h_{ij})_{m \times k}[/math],
где [math]h_{ij} = P(\theta_j | x_i)[/math] — вероятность того, что [math]x_i[/math] пренадлежит [math]j[/math]-ой компоненте.

По формуле Байеса справедливо равенство:
[math] h_{ij} = \frac{w_jp_j(x_i)}{p (x_i)} = \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^k w_s p_s(x_i)}[/math].

Также [math]\sum\limits_{j=1}^k h_{ij} = 1[/math].

Таким образом, зная значения вектора параметров [math]\Theta[/math], мы легко можем пересчитать значения скрытых переменных.

M-шаг

Теорема:

Если известны скрытые переменные, то задача минимизации [math]Q(\Theta)[/math] сводится к [math]k[/math] независимым подзадачам: [math]\theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta)[/math]. Оптимальные же веса считаются как: [math] w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}[/math].

Доказательство:

[math]\triangleright[/math]

Посчитаем логарифм правдоподобия: [math] Q(\Theta) = \sum\limits_{i=1}^m ln\sum\limits_{j=1}^k w_j p_j(x_i; \theta_j) \longrightarrow \max\limits_{\Theta}[/math]. При условии, что[math] \sum\limits_{j=1}^k w_j = 1; w_j \geq 0[/math] имеет смысл рассматривать Лагранжиан задачи: [math] L(\Theta, X^m) = \sum\limits_{i=1}^m ln \biggl( \sum\limits_{j=1}^k w_j p_j(x_i) \biggr) - \lambda \biggl(\sum\limits_{j=1}^k w_j - 1 \biggr) [/math]. Приравняв нулю производную Лагранжиана по [math]w_j[/math], получим: [math]\frac{\partial L} {\partial w_j} = \sum\limits_{i=1}^m \frac{p_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} - \lambda = 0, j = 1..k[/math]. Умножим на [math]w_j[/math] и просуммируем уравнения для всех [math]j[/math]: [math]\sum\limits_{j=1}^k \sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_s p_s(x_i)} = \lambda \sum\limits_{j=1}^kw_j[/math]. А так как [math]\sum\limits_{j=1}^k \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = 1[/math] и [math]\sum\limits_{j=1}^kw_j = 1[/math], из чего следует [math]\lambda = m[/math]. [math]w_j = \frac{1}{m}\sum\limits_{i=1}^m \frac{w_jp_j(x_i)}{\sum\limits_{s=1}^kw_sp_s(x_i)} = \frac{1}{m}\sum\limits_{i=1}^m h_{ij}[/math]. Приравняв к нулю производную Лагранжиана по [math]\theta_j[/math], схожим способом найдем: [math] \theta_j = \arg\max\limits_{\theta}\sum\limits_{i=1}^m h_{ij}*\ln\phi(x_i;\theta).[/math]

[math]\triangleleft[/math]

Критерий остановки

Алгоритм EM выполняется до сходимости, но как нам определить, что сходимость наступила? Мы можем останавливаться, когда либо [math]Q(\Theta)[/math], либо [math]H[/math] перестают сильно меняться. Но, обычно, удобней контролировать изменения значений скрытых переменных, так как они имеют смысл вероятностей и принимают значения из отрезка [math][0,1][/math]. Поэтому один из возможных критериев остановки будет выглядеть так: [math]\max\limits_{i,j} |h_{ij} - h_{ij}^{(0)}| \gt \delta[/math].

Псевдокод

Input:[math]X^m, k, \Theta^{(0)}[/math]
Repeat
   E-step: for all i = 1..m; j = 1..k:
       [math]h_{ij} = \frac{w_j \phi(x_i; \theta_j)}{\sum\limits_{s=1}^k w_s \phi(x_i; \theta_j)}[/math]
   M-step: for all j = 1..k:
       [math]\theta_j = \arg\max\limits_{\theta} \sum\limits_{i=1}^m h_{ij}*ln \phi (x_i, \theta)[/math]
       [math]w_j = \frac {1} {m} \sum\limits_{i=1}^m h_{ij}[/math]
Until a stopping criterion is satisfied
Return [math]\Theta = (\theta_j, w_j)_{j=1}^k[/math]

Плюсы и минусы

Плюсы:

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

Минусы:

• Чуствителен к начальному приближению. Могут быть ситуации, когда сойдемся к локальному экстремуму.
• Число компонент [math]k[/math] является гиперпараметром.

Модификации

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

Generalized EM-algorithm

Осоновная идея этой модификации заключается в том, что на шаге M мы не будем пытаться найти наилучшее решение. Это применимо в случаях, когда максимизация [math]Q(\Theta)[/math] является сликшом дорогой, поэтому нам достаточно сделать лишь несколько итераций, для того, чтобы сместиться в сторону максимума значения [math]Q(\Theta)[/math]. Эта модификация имеет неплохую сходимость.

Stochastic EM-algorithm

Как уже было отмечено в Плюсы и минусы, базовый алгоритм чувствителен к начальному приближению и могут быть ситуации, когда алгоритм "застрянет" в локальном экстремуме. Для того, чтобы предотвратить это, будем на каждой итерации алгоритма случайно "встряхивать" выборку. В этой модификации у нас добавляется шаг S, на котором мы и будем "встряхивать" выборку. И на шаге M мы будем решать уже задачу максимуму невзвешенного правдоподобия. Эта модификация хороша тем, что нечуствиетльная к начальном приблежению.

Пример. Разделение смеси Гауссиана

Несколько итераций алгоритма

Каноническим примером использования EM алгоритма является задача разделения смеси гауссиана. Данные у нас получены из нормального распределения. В этом случае параметрами функций ялвяются матожидание и дисперсия.

[math]\theta = (w_1,..,w_k;\;\mu_1,..,\mu_k;\;\sigma_1,..,\sigma_k)[/math] — вектор параметров,
[math]p_j(x) = N(x;\mu_j, \sigma_j) = \frac1{\sqrt{2\pi}\sigma_j} \exp \biggl(-\frac{(x - \mu_j)^2}{2\sigma_j^2}\biggr) [/math] — плотность распределения.

Посчитаем значения для каждого шага.

E-шаг:

[math] h_{ij} = \frac{w_j N(x_i, \mu_j, \sigma_j)}{\sum\limits_{s=1}^k w_s N(x_i, \mu_s, \sigma_s)}.[/math]

M-шаг:

[math]w_j = \frac{1}{m} \sum\limits_{i=1}^m h_{ij}.[/math]

[math] \mu_j = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}x_i.[/math]

[math] \sigma_j^2 = \frac {1} {mw_j} \sum\limits_{i=1}^m h_{ij}(x_i - \mu_j)^2, j = 1..k.[/math]

Использование в задаче кластеризации

Пример работы k-means

Как уже упоминалось в Определении, алгоритм EM подходит для решения задачи кластеризации. И одной из его имплементаций для этой задачи является алгоритм [math]k[/math]-Means. В этом алгоритме в качестве скрытых переменных выступают метки классов объектов. Параметрами же являются центроиды искомых классов. Тогда на шаге E мы относим объекты к какому-то одному классу на основе расстояний до центроид. А на шаге M мы пересчитываем центроиды кластеров, исходя из полученной на шаге E разметке.

Также стоит упомянуть алгоритм [math]c[/math]-means. В нем качестве скрытых переменных выступают вероятности принадлежности объекта к классам. На шаге E мы пересчитывем вероятности принадлежности объектов, иходя из расстояния до центроид. Шаг M, идейно, остается без изменений.

Реализация на python

В пакете sklearn алгоритм EM представлен объектом GaussianMixture. Проиллюстрируем его работу на примере задачи кластеризации и сравним его с алгоритмом [math]k[/math]-means:

Результат выполнения программы

import numpy as np
import matplotlib.pyplot as plt
from sklearn import cluster, datasets, mixture
from sklearn.preprocessing import StandardScaler
from itertools import cycle, islice
np.random.seed(12)

# Создаем datasets с использованием стандартных sklearn.datasets
n_samples = 2000
random_state = 170
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5, noise=.05)
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)
varied = datasets.make_blobs(n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state)

# Создаем анизатропно разделенные данные
X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
transformation = [[0.6, -0.6], [-0.4, 0.8]]
X_aniso = np.dot(X, transformation)
aniso = (X_aniso, y)

# Выставляем параметры для matplotlib.pyplot
plt.figure(figsize=(9 * 2 + 3, 12.5))
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
plot_num = 1
defaul_n = 3

# Варьируем значение количества классов в зависимости от данных, ведь для нас это гиперпараметр
datasets = [
   (varied, defaul_n),
   (aniso, defaul_n),
   (blobs, defaul_n),
   (noisy_circles, 2)]
for i_dataset, (dataset, n_cluster) in enumerate(datasets):
   X, y = dataset

   # Нормализация данных
   X = StandardScaler().fit_transform(X)

   # Непосредственно наш алгоритм - Gaussian Mixture
   gmm = mixture.GaussianMixture(n_components=n_cluster, covariance_type='full')
   
   # Для сравнения берем алгоритм - k-means
   two_means = cluster.KMeans(n_clusters=n_cluster)
   clustering_algorithms = {
       'Real distribution': None,
       'Gaussian Mixture': gmm,
       'k-Means': two_means
   }
   for name, algorithm in clustering_algorithms:

       # Этап обучения
       if algorithm is not None:
           algorithm.fit(X)
       
       # Применяем алгоритм
        y_pred = y if algorithm is None else algorithm.predict(X)
       
       # Рисуем результаты
       plt.subplot(len(datasets), len(clustering_algorithms), plot_num)
       if i_dataset == 0:
           plt.title(name, size=18)
       colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a']), int(max(y_pred) + 1))))
       plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])
       plt.xlim(-2.5, 2.5)
       plt.ylim(-2.5, 2.5)
       plt.xticks(())
       plt.yticks(())
       plot_num += 1
plt.show()

Как и следовало ожидать, алгоритм EM работает на некоторых данных лучше чем k-means, однако есть данные, с которыми он не справляется без дополнительных преобразований.

См. также

• Кластеризация
• Алгоритм k-Means

Источники информации

1. Материалы лекции про кластеризацию курса "Машинное обучение" университета ИТМО, 2019 год
2. Математические методы обучения по прецедентам К. В. Воронцов
3. Статья про EM-алгоритм на machinelearning.ru
4. A Gentle Introduction to Expectation-Maximization
5. k-means