Contents

Лекция

Для чего нужен EM-алгоритм?

Хотя вариации этого алгоритма используются, например, в байесовском подходе, сейчас он интересует нас в первую очередь как способ максимизировать сложную функцию правдоподобия \(l(x | \theta)\). Здесь \(x\) - наблюдения, \(\theta\) - вектор параметров, которые мы хотим оценить. EM-алгоритм заменяет сложную функцию \(l(x | \theta)\) (которая называется неполным правдоподобием) на легко максимизируемую \(l(x, z | \theta)\) (которая называется полным правдоподобием) и максимизирует её. Здесь \(z\) - некая латентная, ненаблюдаемая переменная. Несмотря на то, что интуитивно первая функция \(l(x | \theta)\) кажется менее громоздкой, часто бывает, что итеративно работать проще с функцией, зависящей также от латентной переменной. Что именно мы будем использовать как \(z\), зависит от конкретной задачи. Но поэтому EM-алгоритм особенно хорошо проявляет себя в тех случаях, когда в задаче интуитивно и логично выделяется такая ненаблюдаемая переменная.

Устройство EM-алгоритма

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

  • E-шага (expectation-step), шага выведения функции \(l(x, z | \theta)\).

  • M-шага (maximization-step), шага ее максимизации.

Условием останова алгоритма может быть как достижение сходимости, так, например, и обычное превышение максимально допустимого количества итераций. Опишем теперь работу отдельных частей алгоритма подробнее.

  • Инициализация

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

\[ \theta_{old} = [\text{случайные/эвристика/}0] \]

  • E-шаг, часть 1

На этом этапе мы хотим найти условное распределение \(p(z | x, \theta_{old})\), где \(\theta_{old}\) - последние посчитанные значения вектора параметров (или, на стартовой итерации, значения, которыми он был проинициализирован).

  • E-шаг, часть 2

Определим функцию, по которой будем оптимизировать новый вектор параметров \(\theta\) относительно старого \(\theta_{old}\):

\[ Q(\theta, \theta_{old}) = \mathbb{E}\left[ l(x, z, |\theta) | x, \theta_{old} \right] \]

  • M-шаг

Максимизируем описанную выше функцию \(Q(\theta, \theta_{old})\) по \(\theta\). В результате получаем новый вектор параметров \(\theta_{new}\), который и будем использовать как \( \theta_{old}\) на следующей итерации. Или, если достигли условия останова, прекратим работу алгоритма и возвратим \(\theta_{new}\) как ответ.

Пример работы

Чтобы лучше понять, как работает EM-алгоритм, решим задачу кластеризации с известным числом кластеров. Задача формулируется так: пусть есть набор независимых наблюдений, вектор \(x\), и все наблюдения принадлежат одному из двух кластеров:

  • Кластер 1, в котором наблюдения \(x_i \sim \mathcal{N}(\mu_1,\sigma_1^2)\). Вероятность того, что случайное наблюдение из смеси относится к первому кластеру равна \(p_1\).

  • Кластер 2, в котором \(x_i \sim \mathcal{N}(\mu_2,\sigma_2^2)\). Вероятность того, что случайное наблюдение из смеси относится ко второму равна \(p_2=1-p_1\).

Другими словами, если мы хотим получить смесь из \(N\) наблюдений, то мы \(N\) раз с вероятностью \(p_1\) генерируем наблюдение согласно распределению первого кластера, иначе - согласно распределению второго. После этого мы хотим, получив как данные только саму смесь и число кластеров (в нашем примере 2), научиться восстанавливать неизвестные параметры \(\theta = (\mu_1, \sigma^2_1, \mu_2, \sigma^2_2, p_1)\).

Попробуем разобраться, почему трудно решать эту задачу методом максимального правдоподобия. Вспомним функцию плотности нормального распределения:

\[ f(x) = \dfrac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\dfrac{(x-\mu)^2}{2\sigma^2}\right) \]

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

\[ p(x_i | \theta) = p_1 \cdot \dfrac{1}{\sqrt{2\pi\sigma_1^2}}\exp\left(-\dfrac{(x_i-\mu_1)^2}{2\sigma_1^2}\right) + (1 - p_1) \cdot \dfrac{1}{\sqrt{2\pi\sigma_2^2}}\exp\left(-\dfrac{(x_i-\mu_2)^2}{2\sigma_2^2}\right) \]
\[ l(x|\theta) = \ln p(x | \theta) = \sum_{i = 1}^n \ln p(x_i | \theta) = \sum_{i = 1}^n \ln [ p_1 \cdot f_1(x_i) + (1 - p_1) \cdot f_2(x_i) ] \]

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

Инициализация

Зададим разумный стартовый набор параметров с помощью простейших эвристик. Самое простое, что мы можем сказать о смеси двух распределений на прямой, что данные слева принадлежат одному кластеру, данные справа - второму. Тогда в качестве параметров матожиданий можно взять минимум и максимум по всей выборке, а в качестве дисперсий - среднее минимума и максимума по всей выборке. Кластеры на данном шаге можно считать равновероятными, то есть вероятность попадания в тот или иной кластер равна \(\dfrac{1}{K}\), где \(K\)- число кластеров. В нашей задаче получаем:

\[\begin{split} \theta_{old} = \begin{cases} \mu_1 = \min(x) \\ \mu_2 = \max(x) \\ \sigma_1 = \sigma_2 = \dfrac{\max(x)-\min(x)}{2} \\ p_1 = \dfrac{1}{2} \end{cases} \end{split}\]

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

Главное - как-то разумно инициализировать стартовый набор параметров, а ЕМ-алгоритм дальше будет его улучшать. Чуть позже мы докажем, что ЕМ-алгоритм точно не будет его ухудшать: в отличие от градиентного спуска, который может перескочить экстремум и уменьшить значение целевой функции, EM-алгоритм никогда не уменьшает ее значение.

Но понятно, что и EM-алгоритм не является панацеей от всех бед, так как ни один алгоритм численной оптимизации не гарантирует сходимость к глобальному экстремуму.

E-шаг, часть 1

Напомним, что главная идея ЕМ-алгоритма: максимизировать не функцию правдоподобия \(l(x | \theta)\), а функцию \(l(x, z | \theta)\) с латентной переменной \(z\). В этой задаче разумно объявить \(z\) как переменную, отвечающую за номер кластера. \(z = (z_1, z_2, .., z_n)\) - вектор всех \(z_i\): \(z_i = 1\), если наблюдение \(x_i\) принадлежит первому кластеру, \(z_i = 2\) - если второму. В реальных данных мы не знаем значения этих переменных (иначе самой задачи классификации бы не было), но мы можем посчитать условные вероятности \(p(z | x, \theta_{old})\).

Что это значит? Это значит, что хоть мы и не можем восстановить сами значения \(z_i\) для каждого наблюдения \(x_i\), мы можем посчитать вероятности \(p(z_i = k | x_i, \theta_{old})\) принадлежности данного наблюдения к кластеру \(k\). Заметим, что для \(K\) кластеров необходимо для каждого наблюдения \(x_i\) посчитать \(K-1\) вероятность, так как вероятность принадлежности последнему кластеру равна \(p(z_i = K | x_i, \theta_{old}) = 1 - \sum_{k = 1}^{K-1} p(z_i = k | x_i, \theta_{old})\).

Обобщая, посчитать закон распределения \(p(z| x, \theta_{old})\) в случае дискретного \(z\) означает посчитать в явном виде \(N*(K-1)\) вероятностей \(p(z_i = k | x_i, \theta_{old})\), где \(N\) - число наблюдений, \(K\) - число кластеров. Что делать в случае, когда латентная переменная является непрерывной? Необходимо посчитать некоторую аппроксимацию непрерывного распределения, так как оно тяжело представляется на компьютере, в отличие от дискретного, которое представляется вектором вероятностей.

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

Итак, вернемся к задаче кластеризации. Нам достаточно посчитать вероятности принадлежности первому кластеру для каждого наблюдения \(p(z_i = 1 | x_i, \theta_{old})\). Для этого воспользуемся формулой условной вероятности:

\[ P(A | B) = \dfrac{P(A \cap B)}{P(B)} \]

В нашем случае:

\[ p(z| x, \theta_{old}) = \dfrac{p(z, x | \theta_{old})}{p(x | \theta_{old})} \]

Тогда:

\[ p(z_i = 1 | x_i, \theta_{old}) = \dfrac{p(z_i = 1, x_i | \theta_{old})}{p(x_i | \theta_{old})} \]

Заметим, что \(\theta_{old}\) нам известны! То есть мы знаем параметры обоих распределений и вероятности принадлежности кластерам. Тогда \(p(z_i = 1, x_i | \theta_{old}) = p_1 \cdot p(x_i |z_i = 1, \theta_{old})\), а \(p(x_i | \theta_{old})\) можно расписать по формуле полной вероятности. Получаем:

\[ p(z_i = 1 | x_i, \theta_{old}) = \dfrac{p_1 \cdot p(x_i |z_i = 1, \theta_{old})}{p_1 \cdot p(x_i |z_i = 1, \theta_{old}) + (1- p_1) \cdot p(x_i |z_i = 2, \theta_{old})} \]

В данной формуле нам известны все значения: \(p_1\) берется из \(\theta_{old}\), а \(p(x_i |z_i = 1, \theta_{old})\) и \(p(x_i |z_i = 2, \theta_{old})\) - это значения функции плотности нормального распределения в точке \(x_i\) с параметрами \(\mu_1\), \(\sigma_1\) и \(\mu_2\), \(\sigma_2\) соответственно, которые также берутся из вектора \(\theta_{old}\). Итого, мы можем в явном виде посчитать \(p(z_i = 1 | x_i, \theta_{old})\)!

E-шаг, часть 2

Определим для нашей задачи функцию, по которой будем оптимизировать новый вектор параметров \(\theta\) относительно старого \(\theta_{old}\):

\[ Q(\theta, \theta_{old}) = \mathbb{E}\left[ l(x, z, |\theta) | x, \theta_{old} \right] = \mathbb{E}\left[ \ln p(x, z |\theta) | x, \theta_{old} \right] \]

На этом шаге мы вычисляем математическое ожидание, используя вероятности \(p(x, z |\theta_{old})\), вычисленные в предыдущем пункте. Обратим внимание, что мы ничего не «предрассчитываем», а лишь определяем функцию для её дальнейшей максимизации.

Вспомним, что логарифм функции плотности для нормального распределения равен

\[ \ln f(x) = -\dfrac{1}{2}\ln(2\pi\sigma^2) - \dfrac{(x-\mu)^2}{2\sigma^2} \]

Теперь у нас есть все необходимое, чтобы расписать формулу \(Q(\theta, \theta_{old})\):

(1)\[\begin{multline} Q(\theta, \theta_{old}) = \mathbb{E}\left[ \ln p(x, z |\theta) | x, \theta_{old} \right] = \\ \sum_{i = 1}^n P(z_i = 1 | x, \theta_{old}) \cdot \ln[p_1 \cdot f_1(x)] + (1 - P(z_i = 1 | x, \theta_{old})) \cdot \ln[p_2 f_2(x)] = \\ \sum\limits_{i = 1}^n \underbrace{P(z_i = 1 | x, \theta_{old})}_{\text{посчитано ранее }} \Big[-\dfrac{1}{2} \ln(2\pi\sigma^2_1) - \dfrac{(x_i-\mu_1)^2}{2\sigma^2_1} + \ln p_1 \Big] + \\ + \underbrace{(1 - P(z_i = 1 | x, \theta_{old}))}_{\text{посчитано ранее }} \Big[ -\dfrac{1}{2}\ln(2\pi\sigma^2_2) - \dfrac{(x_i-\mu_2)^2}{2\sigma^2_2} + \ln (1 - p_1) \Big] \end{multline}\]

Параметры для \(\theta_{new}\) будем знать после максимизации, пока что только определили функцию.

Для общего случая, можем определить \(Q(\theta, \theta_{old})\) через интеграл:

\[ Q(\theta, \theta_{old}) = \int\limits_{R^d} \ln p(x, z | \theta) \cdot p(z | x. \theta_{Old}) \,dz \]

Эта формула корректна для зависимых наблюдений, но в большинстве задач наблюдения независимы, поэтому логарифм по выборке \(\ln p(x, z | \theta)\) будет распадаться в сумму по отдельным наблюдениям.

М-шаг

Наконец, на этом шаге, алгоритмически максимизируем по \(\theta\) функцию, полученную ранее, и находим оптимальную \(\theta_{new}\):

\[ \theta_{new} = \arg\max_{\theta} Q(\theta, \theta_{old}) \]

Заметим, что благодаря нашим преобразованиям, максимумы отдельных переменных из вектора \(\theta_{new}\) легко вычисляются. Для демонcтрации, выведем \(\mu_1^{new}\):

\[ \dfrac{\partial Q}{\partial \mu_1} = \sum\limits_{i = 1}^n P(z_i = 1 | x, \theta_{old}) \dfrac{(x_i-\mu_1)^2}{\sigma^2_1} \]
\[ \mu_1^{new} = \dfrac{\sum P(z_i = 1 | x, \theta_{old}) x_i}{\sum P(z_i = 1 | x, \theta_{old})} \]

Дополнительные материалы

Благодарности

  • Григорьев Пётр

  • Сибагатова Софья

  • Терехова Юлия

Задачи к квизу Семинар