blumental May 12 2020 at 13:38

EM-алгоритм для кластеризации

3 min

7.1K

Technotext 2020

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

Задача

Множество точек $X= \{ x_i, i\in1..N \}$ нужно разбить на $inline$ кластеров.

Идея решения

Составим вероятностную модель распределения точек по кластерам. Найдём параметры модели, при которых вероятность наблюдать множество $inline$ максимальна. С этими параметрами мы сможем определять, к какому кластеру наиболее вероятно относится данная точка $inline$ .

Модель данных

Введем ряд обозначений, заимствованных из курса.

$inline$ — вероятность наблюдать точку $inline$ .

$p(X) = \prod_{i=1}^{N}p(x_i)$ — вероятность наблюдать множество $inline$ .

$p_j (x) = \varphi(x; \theta_j)$ — вероятность встретить точку $inline$ в кластере $inline$ . Это распределение параметризовано параметром (или вектором параметров) $\theta_j$ , индивидуальным для кластера $inline$ .

$inline$ — вероятность кластера $inline$ , т.е. вероятность того, что случайно выбранная точка относится к кластеру $inline$ . Случайно выбранная точка точно относится к какому-то кластеру, поэтому $\sum_{j=1}^K w_j = 1$ .

Из определений выше следует, что $p(x) = \sum_{j=1}^K w_j p_j(x) = \sum_{j=1}^K w_j \varphi(x; \theta_j)$ , т.е. распределение точек моделируется как смесь распределений кластеров.

В итоге, вероятностная модель множества точек $inline$ :

$p(X) = \prod_{i=1}^{N}\left(\sum_{j=1}^K w_j \varphi(x_i; \theta_j)\right)$

Поиск параметров

Параметры модели $inline$ и $\theta$ , как и обсуждалось выше, должны обеспечить максимальную вероятность нашим данным:

$w, \theta = \textrm{argmax} \ p(X) = \textrm{argmax} \ \log p(X) = \textrm{argmax}_{w, \theta} \sum_{i=1}^{N} \log \left(\sum_{j=1}^K w_j \varphi(x_i; \theta_j)\right)$

Сумма под знаком логарифма мешает решить задачу аналитически. Ограничение $\sum_{j=1}^K w_j = 1$ не дает нам просто применить автоматическое дифференцирование (например, TensorFlow или PyTorch).

Придется использовать трюк:

L := нижняя граница log p(X)
while log p(X) увеличивается заметно:
    приближаем L к log p(X)
    w, theta = argmax L

Иначе говоря, мы не будем оптимизировать сам $\log p(X)$ , раз это сложно. Мы возьмем его нижнюю границу $\mathcal{L}$ и будем итеративно делать два шага:

Уточнять $\mathcal{L}$ : "подтягивать" её вверх, ближе к $\log p(X)$ .
Искать $inline$ и $\theta$ , максимизирующие $\mathcal{L}$ .

Мы надеемся, что если полученные параметры максимизируют "близкую" нижнюю границу функции, то они неплохо максимизируют и саму функцию.

Нижняя граница $\mathcal{L}$

Ограничим снизу выражение:

$\log p(X) = \sum_{i=1}^{N} \log \left(\sum_{j=1}^K w_j \varphi(x_i; \theta_j)\right)$

Сначала дополним нашу вероятностную модель. Введем распределение вероятностей $inline$ данной точки $inline$ быть в том или ином кластере:

$g_i(j) \equiv p(\textrm{быть в кластере} \ j| \textrm{это точка} \ i)$

Домножим и поделим на это распределение слагаемые под логарифмом:

$\sum_{i=1}^{N} \log \left(\sum_{j=1}^K w_j \varphi(x_i; \theta_j)\right) =\sum_{i=1}^{N} \log \left(\sum_{j=1}^K \frac{ g_i(j) }{ g_i(j) } w_j \varphi(x_i; \theta_j)\right)$

Теперь применим неравенство Йенсена. Оно позволяет логарифм взвешенной суммы ограничить снизу взвешенной суммой логарифмов:

$\log \left(\sum_i q_i x_i \right) \geq \sum_i q_i \log x_i$

Чтобы неравентсво выполнялось, нужно чтобы веса $inline$ были неотрицательны и их сумма была $inline$ .

В нашем случае $inline$ будет весом: его значения неотрицательны и $\sum_{j=1}^K g_i(j) = 1$ . Применим неравенство:

$\sum_{i=1}^{N} \log \left(\sum_{j=1}^K \frac{ g_i(j) }{ g_i(j) } w_j \varphi(x_i; \theta_j)\right) \geq \sum_{i=1}^{N} \sum_{j=1}^K g_i(j) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g_i(j) }\right)$

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

$\mathcal{L}(g, w, \theta) \equiv \sum_{i=1}^{N} \sum_{j=1}^K g_i(j) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g_i(j) }\right)$

Уточняем $\mathcal{L}$ (E-шаг)

Попробуем максимально приблизить $\mathcal{L}(g, w, \theta)$ к $\log p(X)$ . Параметры $inline$ и $\theta$ будем считать фиксированными, а приближать $\mathcal{L}$ будем через выбор распределения $inline$ .

Запишем разность $\log p(X)$ и $\mathcal{L}$ , а потом посмотрим, как её уменьшить:

$\log p(X) - \mathcal{L}(g, w, \theta) = \sum_{i=1}^N \log p(x_i) - \sum_{i=1}^{N} \sum_{j=1}^K g_i(j) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g_i(j) }\right)=$

$= \sum_{i=1}^N \left(\log p(x_i) \sum_{j=1}^K g_i(j) - \sum_{j=1}^K g_i(j) \log \frac{w_j \varphi(x_i; \theta_j)}{g_i(j)} \right) = \sum_{i=1}^N \sum_{j=1}^K g_i(j) \log \frac{p(x_i) g_i(j)}{w_j \varphi(x_i; \theta_j)}$

Заметим, что под знаком логарифма можно выделить апостериорную вероятность кластера $inline$ :

$p(j|x_i) = \frac{\varphi(x_i; \theta_j) w_j}{p(x_i)}$

Тогда:

$\sum_{i=1}^N \sum_{j=1}^K g_i(j) \log \frac{p(x_i) g_i(j)}{w_j \varphi(x_i; \theta_j)} = \sum_{i=1}^N \sum_{j=1}^K g_i(j) \log \frac{g_i(j)}{p(j|x_i)}= \sum_{i=1}^N \mathbb{E}_{g_i} \frac{g_{i}}{p(j|x_i)}$

Мы получили интересное выражение: матожидание отношения двух распределений по первому из них. Оно называется дивергеницией Кульбака-Лейблера (или кратко KL-дивергенцией) и служит "расстоянием" между вероятностными распределениями.

Таким образом, разность $\log p(X)$ и $\mathcal{L}$ — это сумма KL-дивергенций:

$\log p(X) - \mathcal{L}(g, w, \theta) = \sum_{i=1}^N KL(g_i || p(j|x_i))$

KL-дивергенция неотрицательна, поэтому лучшее, что мы можем сделать для приближения нижней границы — это сделать KL-дивергенцию нулевой. А этого несложно добиться: KL-дивергенция будет нулём, если её аргументы — это одинаковые распределения. Вот и сделаем распределение $inline$ тождественным $inline$ :

$g_i(j) = p(j|x_i) = \frac{w_j \varphi(x_i; \theta_j)}{p(x_i)}$

Вычисление $inline$ по данной формуле и позволит нам приблизить нижнюю границу $\mathcal{L}$ к $\log p(X)$ .

Максимизируем $\mathcal{L}$ по параметрам (M-шаг)

Теперь вторая часть итерации: поиск параметров по нижней границе. В этой части наши предположения будут противоположными:

распределение $inline$ фиксированно;
параметры $inline$ и $\theta$ подлежат оптимизации.

Перед оптимизацией немного упростим $\mathcal{L}$ :

$\mathcal{L}(g, \theta) = \sum_{i=1}^N\left( \sum_{j=1}^K g_i(j) \log \frac{w_j p(x_i; \theta_j)}{g_i(j)} \right) =$

$= \sum_{i=1}^N \sum_{j=1}^K g_i(j) \log \left( w_j p(x_i; \theta_j) \right) -\sum_{i=1}^N \sum_{j=1}^K g_i(j) \log g_i(j)$

Второе слагаемое не зависит от параметров $inline$ и $\theta$ , поэтому далее будем оптимизировать только первое слагаемое:

$w, \theta = \textrm{argmax}_{w, \theta}\sum_{i=1}^N \sum_{j=1}^K g_i(j) \log \left( w_j \varphi(x_i; \theta_j) \right)$

Разложим логарифм произведения на сумму логарифмов и получим:

$w = \textrm{argmax}_{w}\sum_{i=1}^N \sum_{j=1}^K g_i(j) \log w_j, \textrm{ при условии }\sum_{j=1} w_j = 1$

$\theta_j = \textrm{argmax} \sum_{i=1}^N g_i(j) \log \varphi (x_i; \theta_j)$

Первая задача решается методом множителей Лагранжа. Результат:

$w_j = \frac{1}{N} \sum_{i=1}^N g_i(j)$

Решение второй задачи зависит от конкретного вида распределения кластера $\varphi (x_i; \theta_j)$ . Как видно, для её решение больше не придётся иметь дело с суммой под знаком логарифма, поэтому, например, для гауссовых распределений решение может быть выписано аналитически.

Итог

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

Tags:

Hubs:

Machine learning