Оффлайн
- Регистрация
- 12.04.17
- Сообщения
- 19.095
- Реакции
- 107
- Репутация
- 0
EM-алгоритм – полезный инструмент моделирования данных, когда максимизация правдоподобия "в лоб", через дифференцирование, невозможна. Кластеризация – одна из задач, где этот алгоритм приходит на помощь. В статье приведен общий вывод EM-алгоритма для кластеризации.
Задача
Множество точек $inline$X= \{ x_i, i\in1..N \}$inline$ нужно разбить на $inline$K$inline$ кластеров.
Идея решения
Составим вероятностную модель распределения точек по кластерам. Найдём параметры модели, при которых вероятность наблюдать множество $inline$X$inline$ максимальна. С этими параметрами мы сможем определять, к какому кластеру наиболее вероятно относится данная точка $inline$x$inline$.
Модель данных
Введем ряд обозначений, заимствованных из
$inline$p(x)$inline$ — вероятность наблюдать точку $inline$x$inline$.
$inline$p(X) = \prod_{i=1}^{N}p(x_i)$inline$ — вероятность наблюдать множество $inline$X$inline$.
$inline$p_j (x) = \varphi(x; \theta_j)$inline$ — вероятность встретить точку $inline$x$inline$ в кластере $inline$j$inline$. Это распределение параметризовано параметром (или вектором параметров) $inline$\theta_j$inline$, индивидуальным для кластера $inline$j$inline$.
$inline$w_j$inline$ — вероятность кластера $inline$j$inline$, т.е. вероятность того, что случайно выбранная точка относится к кластеру $inline$j$inline$. Случайно выбранная точка точно относится к какому-то кластеру, поэтому $inline$\sum_j w_j = 1$inline$.
Из определений выше следует, что $inline$p(x) = \sum_j w_j p_j(x) = \sum_j w_j \varphi(x; \theta_j)$inline$, т.е. распределение точек моделируется как смесь распределений кластеров.
В итоге, вероятностная модель множества точек $inline$X$inline$:
$$display$$p(X) = \prod_i^{N}\left(\sum_j^K w_j \varphi(x_i; \theta_j)\right)$$display$$
Поиск параметров
Параметры модели $inline$w$inline$ и $inline$\theta$inline$, как и обсуждалось выше, должны обеспечить максимальную вероятность нашим данным:
$$display$$w, \theta = \textrm{argmax} \ p(X) = \textrm{argmax} \ \log p(X) = \textrm{argmax}_{w, \theta} \sum_i^{N} \log \left(\sum_j^K w_j \varphi(x_i; \theta_j)\right)$$display$$
Сумма под знаком логарифма мешает решить задачу аналитически. Ограничение $inline$\sum_j w_j = 1$inline$ не дает нам просто применить автоматическое дифференцирование (например, TensorFlow или PyTorch).
Придется использовать трюк:
L := нижняя граница log p(X)
while log p(X) увеличивается заметно:
приближаем L к log p(X)
w, theta = argmax L
Иначе говоря, мы не будем оптимизировать сам $inline$\log p(X)$inline$, раз это сложно. Мы возьмем его нижнюю границу $inline$\mathcal{L}$inline$ и будем итеративно делать два шага:
Мы надеемся, что если полученные параметры максимизируют "близкую" нижнюю границу функции, то они неплохо максимизируют и саму функцию.
Нижняя граница $inline$\mathcal{L}$inline$
Ограничим снизу выражение:
$$display$$ \log p(X) = \sum_i^{N} \log \left(\sum_j^K w_j \varphi(x_i; \theta_j)\right)$$display$$
Сначала дополним нашу вероятностную модель. Введем распределение вероятностей данной точки $inline$x_i$inline$ быть в том или ином кластере:
$$display$$g_{ij} = p(\textrm{быть в кластере} \ j| \textrm{это точка} \ i)$$display$$
Домножим и поделим на это распределение слагаемые под логарифмом:
$$display$$\sum_i^{N} \log \left(\sum_j^K w_j \varphi(x_i; \theta_j)\right) =\sum_i^{N} \log \left(\sum_j^K \frac{ g(ij) }{ g(ij) } w_j \varphi(x_i; \theta_j)\right)$$display$$
Теперь применим
$$display$$\log \left(\sum_i q_i x_i \right) \geq \sum_i q_i \log x_i$$display$$
Чтобы неравентсво выполнялось, нужно чтобы веса $inline$q_i$inline$ были неотрицательны и их сумма была $inline$1$inline$.
В нашем случае $inline$g_{ij}$inline$ будет весом: его значения неотрицательны и $inline$\sum_j^K g(ij) = 1$inline$. Применим неравенство:
$$display$$\sum_i^{N} \log \left(\sum_j^K \frac{ g(ij) }{ g(ij) } w_j \varphi(x_i; \theta_j)\right) \geq \sum_i^{N} \sum_j^K g(ij) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g(ij) }\right)$$display$$
Полученное выражение и будем использовать в качестве нижней границы:
$$display$$\mathcal{L}(g, w, \theta) \equiv \sum_i^{N} \sum_j^K g(ij) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g(ij) }\right)$$display$$
Уточняем $inline$\mathcal{L}$inline$ (E-шаг)
Попробуем максимально приблизить $inline$\mathcal{L}(g, w, \theta)$inline$ к $inline$\log p(X)$inline$. Параметры $inline$w$inline$ и $inline$\theta$inline$ будем считать фиксированными, а приближать $inline$\mathcal{L}$inline$ будем через выбор распределения $inline$g$inline$.
Запишем разность $inline$\log p(X)$inline$ и $inline$\mathcal{L}$inline$, а потом посмотрим, как её уменьшить:
$$display$$\log p(X) - \mathcal{L}(g, w, \theta) = \sum_{i=1}^N \log p(x_i) - \sum_i^{N} \sum_j^K g(ij) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g(ij) }\right)=$$display$$
$$display$$= \sum_i^N \left(\log p(x_i) \sum_j^K g(ij) - \sum_j^K g(ij) \log \frac{w_j \varphi(x_i; \theta_j)}{g(ij)} \right) = \sum_i^N \sum_j^K g(ij) \log \frac{p(x_i) g(ij)}{w_j \varphi(x_i; \theta_j)}$$display$$
Заметим, что под знаком логарифма можно выделить апостериорную вероятность кластера $inline$j$inline$:
$$display$$p(j|x_i) = \frac{\varphi(x_i; \theta_j) w_j}{p(x_i)}$$display$$
Тогда:
$$display$$\sum_i^N \sum_j^K g(ij) \log \frac{p(x_i) g(ij)}{w_j \varphi(x_i; \theta_j)} = \sum_i^N \sum_j^K g(ij) \log \frac{g(ij)}{p(j|x_i)}= \mathbb{E}_g \frac{g(ij)}{p(j|x_i)}$$display$$
Мы получили интересное выражение: матожидание отношения двух распределений по первому из них. Оно называется
Таким образом, разность $inline$\log p(X)$inline$ и $inline$\mathcal{L}$inline$ — это KL-дивергенция между двумя распределениями:
$$display$$\log p(X) - \mathcal{L}(g, w, \theta) = KL(g_{ij} || p(j|x_i)) $$display$$
KL-дивергенция неотрицательна, поэтому лучшее, что мы можем сделать для приближения нижней границы — это сделать KL-дивергенцию нулевой. А этого несложно добиться: KL-дивергенция будет нулём, если её аргументы — это одинаковые распределения. Вот и сделаем распределение $inline$g_{ij}$inline$ тождественным $inline$p(j|x_i)$inline$:
$$display$$g_{ij} = p(j|x_i) = \frac{w_j \varphi(x_i; \theta_j)}{p(x_i)}$$display$$
Вычисление $inline$g_{ij}$inline$ по данной формуле и позволит нам приблизить нижнюю границу $inline$\mathcal{L}$inline$ к $inline$\log p(X)$inline$.
Максимизируем $inline$\mathcal{L}$inline$ по параметрам (M-шаг)
Теперь вторая часть итерации: поиск параметров по нижней границе. В этой части наши предположения будут противоположными:
Перед оптимизацией немного упростим $inline$\mathcal{L}$inline$:
$$display$$\mathcal{L}(g, \theta) = \sum_i^N\left( \sum_j^K g(ij) \log \frac{w_j p(x_i; \theta_j)}{g(ij)} \right) =$$display$$
$$display$$ = \sum_i^N \sum_j^K g(ij) \log \left( w_j p(x_i; \theta_j) \right) -\sum_i^N \sum_j^K g(ij) \log g(ij)$$display$$
Второе слагаемое не зависит от параметров $inline$w$inline$ и $inline$\theta$inline$, поэтому далее будем оптимизировать только первое слагаемое:
$$display$$w, \theta = \textrm{argmax}_{w, \theta}\sum_i^N \sum_j^K g(ij) \log \left( w_j \varphi(x_i; \theta_j) \right)$$display$$
Разложим логарифм произведения на сумму логарифмов и получим:
$$display$$w = \textrm{argmax}_{w}\sum_i^N \sum_j^K g(ij) \log w_j, \textrm{ при условии }\sum_j w_j = 1$$display$$
$$display$$\theta_j = \textrm{argmax} \sum_i^N g(ij) \log \varphi (x_i; \theta_j)$$display$$
Первая задача решается методом множителей Лагранжа. Результат:
$$display$$w_j = \frac{1}{N} \sum_i^N g(ij)$$display$$
Решение второй задачи зависит от конкретного вида распределения кластера $inline$\varphi (x_i; \theta_j)$inline$. Как видно, для её решение больше не придётся иметь дело с суммой под знаком логарифма, поэтому, например, для гауссовых распределений решение может быть выписано аналитически.
Итог
Мы выяснили, в чем заключается суть итераций EM-алгоритма для кластеризации, и увидели, как в общем виде выводятся их формулы.
Задача
Множество точек $inline$X= \{ x_i, i\in1..N \}$inline$ нужно разбить на $inline$K$inline$ кластеров.
Идея решения
Составим вероятностную модель распределения точек по кластерам. Найдём параметры модели, при которых вероятность наблюдать множество $inline$X$inline$ максимальна. С этими параметрами мы сможем определять, к какому кластеру наиболее вероятно относится данная точка $inline$x$inline$.
Модель данных
Введем ряд обозначений, заимствованных из
You must be registered for see links
.$inline$p(x)$inline$ — вероятность наблюдать точку $inline$x$inline$.
$inline$p(X) = \prod_{i=1}^{N}p(x_i)$inline$ — вероятность наблюдать множество $inline$X$inline$.
$inline$p_j (x) = \varphi(x; \theta_j)$inline$ — вероятность встретить точку $inline$x$inline$ в кластере $inline$j$inline$. Это распределение параметризовано параметром (или вектором параметров) $inline$\theta_j$inline$, индивидуальным для кластера $inline$j$inline$.
$inline$w_j$inline$ — вероятность кластера $inline$j$inline$, т.е. вероятность того, что случайно выбранная точка относится к кластеру $inline$j$inline$. Случайно выбранная точка точно относится к какому-то кластеру, поэтому $inline$\sum_j w_j = 1$inline$.
Из определений выше следует, что $inline$p(x) = \sum_j w_j p_j(x) = \sum_j w_j \varphi(x; \theta_j)$inline$, т.е. распределение точек моделируется как смесь распределений кластеров.
В итоге, вероятностная модель множества точек $inline$X$inline$:
$$display$$p(X) = \prod_i^{N}\left(\sum_j^K w_j \varphi(x_i; \theta_j)\right)$$display$$
Поиск параметров
Параметры модели $inline$w$inline$ и $inline$\theta$inline$, как и обсуждалось выше, должны обеспечить максимальную вероятность нашим данным:
$$display$$w, \theta = \textrm{argmax} \ p(X) = \textrm{argmax} \ \log p(X) = \textrm{argmax}_{w, \theta} \sum_i^{N} \log \left(\sum_j^K w_j \varphi(x_i; \theta_j)\right)$$display$$
Сумма под знаком логарифма мешает решить задачу аналитически. Ограничение $inline$\sum_j w_j = 1$inline$ не дает нам просто применить автоматическое дифференцирование (например, TensorFlow или PyTorch).
Придется использовать трюк:
L := нижняя граница log p(X)
while log p(X) увеличивается заметно:
приближаем L к log p(X)
w, theta = argmax L
Иначе говоря, мы не будем оптимизировать сам $inline$\log p(X)$inline$, раз это сложно. Мы возьмем его нижнюю границу $inline$\mathcal{L}$inline$ и будем итеративно делать два шага:
- Уточнять $inline$\mathcal{L}$inline$: "подтягивать" её вверх, ближе к $inline$\log p(X)$inline$.
- Искать $inline$w$inline$ и $inline$\theta$inline$, максимизирующие $inline$\mathcal{L}$inline$.
Мы надеемся, что если полученные параметры максимизируют "близкую" нижнюю границу функции, то они неплохо максимизируют и саму функцию.
Нижняя граница $inline$\mathcal{L}$inline$
Ограничим снизу выражение:
$$display$$ \log p(X) = \sum_i^{N} \log \left(\sum_j^K w_j \varphi(x_i; \theta_j)\right)$$display$$
Сначала дополним нашу вероятностную модель. Введем распределение вероятностей данной точки $inline$x_i$inline$ быть в том или ином кластере:
$$display$$g_{ij} = p(\textrm{быть в кластере} \ j| \textrm{это точка} \ i)$$display$$
Домножим и поделим на это распределение слагаемые под логарифмом:
$$display$$\sum_i^{N} \log \left(\sum_j^K w_j \varphi(x_i; \theta_j)\right) =\sum_i^{N} \log \left(\sum_j^K \frac{ g(ij) }{ g(ij) } w_j \varphi(x_i; \theta_j)\right)$$display$$
Теперь применим
You must be registered for see links
. Оно позволяет логарифм взвешенной суммы ограничить снизу взвешенной суммой логарифмов:$$display$$\log \left(\sum_i q_i x_i \right) \geq \sum_i q_i \log x_i$$display$$
Чтобы неравентсво выполнялось, нужно чтобы веса $inline$q_i$inline$ были неотрицательны и их сумма была $inline$1$inline$.
В нашем случае $inline$g_{ij}$inline$ будет весом: его значения неотрицательны и $inline$\sum_j^K g(ij) = 1$inline$. Применим неравенство:
$$display$$\sum_i^{N} \log \left(\sum_j^K \frac{ g(ij) }{ g(ij) } w_j \varphi(x_i; \theta_j)\right) \geq \sum_i^{N} \sum_j^K g(ij) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g(ij) }\right)$$display$$
Полученное выражение и будем использовать в качестве нижней границы:
$$display$$\mathcal{L}(g, w, \theta) \equiv \sum_i^{N} \sum_j^K g(ij) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g(ij) }\right)$$display$$
Уточняем $inline$\mathcal{L}$inline$ (E-шаг)
Попробуем максимально приблизить $inline$\mathcal{L}(g, w, \theta)$inline$ к $inline$\log p(X)$inline$. Параметры $inline$w$inline$ и $inline$\theta$inline$ будем считать фиксированными, а приближать $inline$\mathcal{L}$inline$ будем через выбор распределения $inline$g$inline$.
Запишем разность $inline$\log p(X)$inline$ и $inline$\mathcal{L}$inline$, а потом посмотрим, как её уменьшить:
$$display$$\log p(X) - \mathcal{L}(g, w, \theta) = \sum_{i=1}^N \log p(x_i) - \sum_i^{N} \sum_j^K g(ij) \log \left(\frac{ w_j \varphi(x_i; \theta_j) }{ g(ij) }\right)=$$display$$
$$display$$= \sum_i^N \left(\log p(x_i) \sum_j^K g(ij) - \sum_j^K g(ij) \log \frac{w_j \varphi(x_i; \theta_j)}{g(ij)} \right) = \sum_i^N \sum_j^K g(ij) \log \frac{p(x_i) g(ij)}{w_j \varphi(x_i; \theta_j)}$$display$$
Заметим, что под знаком логарифма можно выделить апостериорную вероятность кластера $inline$j$inline$:
$$display$$p(j|x_i) = \frac{\varphi(x_i; \theta_j) w_j}{p(x_i)}$$display$$
Тогда:
$$display$$\sum_i^N \sum_j^K g(ij) \log \frac{p(x_i) g(ij)}{w_j \varphi(x_i; \theta_j)} = \sum_i^N \sum_j^K g(ij) \log \frac{g(ij)}{p(j|x_i)}= \mathbb{E}_g \frac{g(ij)}{p(j|x_i)}$$display$$
Мы получили интересное выражение: матожидание отношения двух распределений по первому из них. Оно называется
You must be registered for see links
(или кратко KL-дивергенцией) и служит "расстоянием" между вероятностными распределениями.Таким образом, разность $inline$\log p(X)$inline$ и $inline$\mathcal{L}$inline$ — это KL-дивергенция между двумя распределениями:
$$display$$\log p(X) - \mathcal{L}(g, w, \theta) = KL(g_{ij} || p(j|x_i)) $$display$$
KL-дивергенция неотрицательна, поэтому лучшее, что мы можем сделать для приближения нижней границы — это сделать KL-дивергенцию нулевой. А этого несложно добиться: KL-дивергенция будет нулём, если её аргументы — это одинаковые распределения. Вот и сделаем распределение $inline$g_{ij}$inline$ тождественным $inline$p(j|x_i)$inline$:
$$display$$g_{ij} = p(j|x_i) = \frac{w_j \varphi(x_i; \theta_j)}{p(x_i)}$$display$$
Вычисление $inline$g_{ij}$inline$ по данной формуле и позволит нам приблизить нижнюю границу $inline$\mathcal{L}$inline$ к $inline$\log p(X)$inline$.
Максимизируем $inline$\mathcal{L}$inline$ по параметрам (M-шаг)
Теперь вторая часть итерации: поиск параметров по нижней границе. В этой части наши предположения будут противоположными:
- распределение $inline$g$inline$ фиксированно;
- параметры $inline$w$inline$ и $inline$\theta$inline$ подлежат оптимизации.
Перед оптимизацией немного упростим $inline$\mathcal{L}$inline$:
$$display$$\mathcal{L}(g, \theta) = \sum_i^N\left( \sum_j^K g(ij) \log \frac{w_j p(x_i; \theta_j)}{g(ij)} \right) =$$display$$
$$display$$ = \sum_i^N \sum_j^K g(ij) \log \left( w_j p(x_i; \theta_j) \right) -\sum_i^N \sum_j^K g(ij) \log g(ij)$$display$$
Второе слагаемое не зависит от параметров $inline$w$inline$ и $inline$\theta$inline$, поэтому далее будем оптимизировать только первое слагаемое:
$$display$$w, \theta = \textrm{argmax}_{w, \theta}\sum_i^N \sum_j^K g(ij) \log \left( w_j \varphi(x_i; \theta_j) \right)$$display$$
Разложим логарифм произведения на сумму логарифмов и получим:
$$display$$w = \textrm{argmax}_{w}\sum_i^N \sum_j^K g(ij) \log w_j, \textrm{ при условии }\sum_j w_j = 1$$display$$
$$display$$\theta_j = \textrm{argmax} \sum_i^N g(ij) \log \varphi (x_i; \theta_j)$$display$$
Первая задача решается методом множителей Лагранжа. Результат:
$$display$$w_j = \frac{1}{N} \sum_i^N g(ij)$$display$$
Решение второй задачи зависит от конкретного вида распределения кластера $inline$\varphi (x_i; \theta_j)$inline$. Как видно, для её решение больше не придётся иметь дело с суммой под знаком логарифма, поэтому, например, для гауссовых распределений решение может быть выписано аналитически.
Итог
Мы выяснили, в чем заключается суть итераций EM-алгоритма для кластеризации, и увидели, как в общем виде выводятся их формулы.