EM算法

简介：本文主要介绍以下两个方面内容：

混合高斯(Mixtures of Gaussians)和EM算法
EM算法相关推导证明

混合高斯和EM算法

关于混合高斯可以看漫谈 Clustering (3): Gaussian Mixture Model 和混合高斯模型（Mixtures of Gaussians）和EM算法。观察样本数据X满足公式一概率密度函数：

$p(X\mid\Theta )=\sum_{i=1}^{m}\alpha _{i}p(X\mid\theta _{i})$ 公式1

其中 $\alpha_{i}$ 表示第i个Gaussian选中的概率，满足 $\sum_{i=1}^{k}\alpha _{i}=1$ ，数据满足第i个Guassian分布的概率为 $p_{i}(X\mid\theta_{i})=N(\mu_{i},\Sigma_{i})$ ，那么公式一中未知参数为

$\begin{align*} \Theta &=(\alpha _{1},\cdots ,\alpha _{k},\theta_{1},\cdots,\theta _{k}) \\ &= (\alpha _{1},\cdots ,\alpha _{k},(\mu_{1},\Sigma_{1}),\cdots,(\mu_{k},\Sigma_{k})) \end{align*}$

采用极大似然估计法maximum-likelihood estimation (MLE)求参数，似然函数如公式2所示：

$\begin{align*} ln\, p(X\mid \theta ) &=ln\,\prod_{i=1}^{n}\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j},\Sigma_{j}) \\ &= \sum_{i=1}^{n} ln\,\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j},\Sigma_{j}) \end{align*}$ 公式2

公式2对 $\mu_{r}$ 求导并令其等于0可得：

$\frac{\partial ln\, p(X\mid \theta )}{\partial \mu_{r}} =\sum_{i=1}^{n}\frac{\alpha _{r}N(x_{i}\mid\mu_{r},\Sigma_{r})}{\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j},\Sigma_{j})}\Sigma_{r}^{-1}(x_{i}-\mu_{r})=0$

设一个后验概率(Posterior probability )如公式3所示：

$p(\alpha_{r}\mid x_{i})=\frac{\alpha _{r}N(x_{i}\mid\mu_{r},\Sigma_{r})}{\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j},\Sigma_{j})}$ 公式3

可以得出 $\mu_{r}$ 的表达式如公式4所示：

$\mu_{r}=\frac{\sum_{i=1}^{n}p(\alpha_{r}\mid x_{i})x_{i}}{\sum_{i=1}^{n}p(\alpha_{r}\mid x_{i})}$ 公式4

用拉格朗日乘数法Lagrange multiplier求对 $\alpha _{r}$ 。构造拉格朗日乘子如下：

$\Lambda(\alpha)=\sum_{i=1}^{n} ln\,\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j},\Sigma_{j})+\beta (\sum_{j=1}^{k}\alpha _ {j}-1)$

求导可得

$\begin{align*} \frac{\partial \Lambda(\alpha)}{\partial \alpha _{r}} &= \sum_{i=1}^{n}\frac{N(x_{i}\mid\mu_{r},\Sigma_{r})}{\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j}, \Sigma_{j})}+\beta\\ &= \sum_{i=1}^{n}\frac{p(\alpha_{r}\mid x_{i})}{\alpha_{r}}+\beta\\ &= 0 \end{align*}$

易知 $\beta =\sum_{i=1}^{n}\sum_{j=1}^{k}p(\alpha_{j}\mid x_{i})=\sum_{i=1}^{n}1=n$ ，那么最终可以得到 $\alpha _{r}$ 如公式5所示：

$\alpha_{r}=\frac{1}{n}\sum_{i=1}^{n}p(\alpha_{r}\mid x_{i})$ 公式5

计算 $\Sigma_{r}$ 。这里用到单高斯分布极大似然估计的相关结论(详见Pattern Recognition and Machine Learning 练习2.34)。对于单高斯分布，有如公式6所示结论:

$\frac{\partial ln\, p(x_{i}\mid\mu_{r},\Sigma_{r})}{\partial \Sigma_{r}} =-\frac{1}{2}\Sigma_{r}^{-1}+\frac{1} {2}\Sigma_{r}^{-1}S\Sigma_{r}^{-1}$ 公式6

其中 $S=(x_i-\mu_{r})(x_i-\mu_{r})^{T}$ 。那么公式2对 $\Sigma_{r}$ 求导并令其等于0可得：

$\begin{align*} \frac{\partial ln\, p(X\mid \theta )}{\partial \Sigma_{r}} &= \sum_{i=1}^{n}\frac{\alpha _{r}N(x_{i}\mid\mu_{r},\Sigma_{r})}{\sum_{j=1}^{k}\alpha _{j}N(x_{i}\mid\mu_{j},\Sigma_{j})}(-\frac{1}{2}\Sigma_{r}^{-1}+\frac{1}{2}\Sigma_{r}^{-1}S\Sigma_{r}^{-1})\\ &= -\frac{1}{2}\Sigma_{r}^{-1}p(\alpha_{r}\mid x_{i})(\Sigma_{r}-S)\Sigma_{r}^{-1}\\ &= 0 \end{align*}$

最终整理可以得到 $\Sigma_{r}$ 的表达式如公式7所示：

$\Sigma_{r}=\frac{p(\alpha_{r}\mid x_{i})(x_i-\mu_{r})(x_i-\mu_{r})^{T}}{p(\alpha_{r}\mid x_{i})}$ 公式7

我们得到问题相关参数的表达式如公式4、公式5、公式7所示，但是无法根据样本数据直接获得，因为在计算后验概率 $p(\alpha_{r}\mid x_{i})$ 无法直接计算。EM算法给出这类问题的解决方法，即先初始化一组参数 $\Theta$ ，在E-Step计算后验概率 $p(\alpha_{r}\mid x_{i})$ ，在M步骤更新参数 $\Theta^{new}$ . 如算法一所示：

算法1:EM算法求解混合高斯

Initialize:
$\mu_{r}^{new}$ ， $\alpha _{r}^{new}$ ， $\Sigma_{r}^{new}$
E-Step:
$\mu_{r}\leftarrow \mu_{r}^{new}\, \alpha _{r}\leftarrow \alpha _{r}^{new}\, \Sigma_{r}\leftarrow \Sigma_{r}^{new}$

$p(\alpha _{r}\mid x_{i})\leftarrow \frac{\alpha _{r}^{new}N(x_{i}\mid \mu_{r}^{new},\Sigma_{r}^{new})}{\sum_{j=1}^{k}\alpha _{j}^{new}N(x_{i}\mid \mu_{j}^{new},\Sigma_{j}^{new})}$
M-Step:
$\mu _{r}^{new}\leftarrow \frac{\sum_{i=1}^{n}p(\alpha_{r}\mid x_{i})x_{i}}{\sum_{i=1}^{n}p(\alpha_{r}\mid x_{i})}$

$\alpha_{r}\leftarrow \frac{1}{n}\sum_{i=1}^{n}p(\alpha_{r}\mid x_{i})$

$\Sigma_{r}^{new}\leftarrow \frac{p(\alpha_{r}\mid x_{i})(x_i-\mu_{r})(x_i-\mu_{r})^{T}}{p(\alpha_{r}\mid x_{i})}$

这部分主要是利用混合高斯引出EM算法，下部分内容详细讲解EM算法推导及相关证明。

EM算法

EM算法是用极大似然估计法求解存在隐含变量(如混合高斯中的 $\alpha_{i}$ )问题时的一种有效的方法。下面假定观察样本为X，隐含变量为H，联合概率密度 $p(X,H\mid \Theta)$ 由参数 $\Theta$ 决定：

$p(X\mid \Theta)= \int_{H}\, p(X,H\mid \Theta)\mathrm{d}H$

我们的目标就是极大化公式8来求参数 $\Theta$ :

$L(\Theta ) &= \sum_{i=1}^{n}ln\, p(x\mid \theta)= \sum_{i=1}^{n}ln\int_{h}\, p(x_{i},h\mid \theta)\mathrm{d}h$ 公式8

直接求解公式8有难度，通常如果隐含变量H知道，那么公式8的求解将会容易一些。

EM算法对此类问题提供了一个有效的解决方法： 重复构造 $L(\Theta )$ 的一个下界(E-Step)，优化(提高)这个下界(M-Step)。如果构造和优化下界的任务比较容易，那么我们就间接的极大化公式8求得参数 $\Theta$ 。

对于每一个i，我们定义隐含变量满足某种分布的概率 $q(h_{i})$ (即满足 $\int_{h}q(h_{i})=1$ 和 $q(h_{i})\geq 0$ )，那么我们可以得到公式9：

$\begin{align*} L(\Theta ) &= \sum_{i=1}^{n}ln\, p(x\mid \theta)= \sum_{i=1}^{n}ln\int_{h}\, p(x_{i},h\mid \theta)\mathrm{d}h\\ &= \sum_{i=1}^{n}ln\int_{h}\, q(h)\frac{p(x_{i},h\mid \theta)}{q(h)}\mathrm{d}h\\ &\geq \sum_{i=1}^{n}\int_{h}\, q(h)ln\,\frac{p(x_{i},h\mid \theta)}{q(h)}\mathrm{d}h \end{align*}$ 公式9

这里公式9最后一步用到了 Jensen不等式，如过f是一个凹函数(如开口向下的抛物线)，X是随机变量，有如下结论：

$E[f(X)]\leq f(E[X])$

等式成立的条件是X=E[X]即X是常量。易知lnx是一个凹函数，那么最后一步推导如下:

$ln\,(E_{h}[\frac{p(x_{i},h\mid \theta)}{q(h)}]) \geq E_{h}ln\,([\frac{p(x_{i},h\mid \theta)}{q(h)}])$

给定一组隐含变量的分布 $q(h_{i})$ ，公式9给出了 $L(\Theta )$ 的一个下界。对于给定一组参数 $\theta$ 等号成立的条件是 $\frac{p(x_{i},h_{i}\mid \theta)}{q(h_{i})}=C$ ，那么我们可以得到 $q(h_{i})$ 如公式10所示：

$\begin{align*} q(h_{i})&= \frac{p(x_{i},h_{i}\mid \theta)}{C}=\frac{p(x_{i},h_{i}\mid \theta)}{\sum_{h}p(x_{i},h\mid \theta)}=\frac{p(x_{i},h_{i}\mid \theta)}{p(x_{i}\mid \theta)}\\ &=p(h_{i}\mid x_{i},\theta) \end{align*}$ 公式10

即 $q(h_{i})$ 定义为给定样本x和参数 $\theta$ 时隐含变量h的后验概率。

给出这个下界后我们可以通过 $\theta$ 来优化这个下界：

$\begin{align*} \widehat{\theta} &= \underset{\theta}{arg\, max}L(\Theta)\\ &= \underset{\theta}{arg\, max}\sum_{i=1}^{n}\int_{h}\, q(h)ln\,\frac{p(x_{i},h\mid \theta)}{q(h)}\mathrm{d}h\\ &= \underset{\theta}{arg\, max}\sum_{i=1}^{n}(\int_{h}\, q(h)ln\,p(x_{i},h\mid \theta)\mathrm{d}h-\int_{h}\, q(h)ln\,q(h)\mathrm{d}h)\\ &= \underset{\theta}{arg\, max}\sum_{i=1}^{n}\int_{h}\, q(h)ln\,p(x_{i},h\mid \theta)\mathrm{d}h \end{align*}$

最终我们可以得到EM算法如算法2所示：

算法2：EM算法过程

Initialize:
$q^{[0]}(h),\theta^{[0]}$
E-Step:
$q^{[t]}(h)\leftarrow p(h\mid x,\theta^{[t-1]})$
M-Step:
$\theta^{[t]}\leftarrow \underset{\theta}{arg\, max}\sum_{i=1}^{n}\int_{h}\, q^{[t]}(h)ln\,p(x_{i},h\mid \theta)\mathrm{d}h$

剩余的最后一个问题是如何证明EM算法会收敛，也就是证明 $L(\Theta^{[t+1]})\geq L(\Theta^{[t]})$ ，证明如公式11所示:

$\begin{align*} L(\Theta^{[t+1]}) &\geq \sum_{i=1}^{n}\int_{h}\, q^{[t]}(h)ln\,p(x_{i},h\mid \theta^{[t+1]})\mathrm{d}h\\ &\geq \sum_{i=1}^{n}\int_{h}\, q^{[t]}(h)ln\,p(x_{i},h\mid \theta^{[t]})\mathrm{d}h\\ &= L(\Theta^{[t]}) \end{align*}$ 公式11

其中第一个不等式利用就是公式9结论，第二个不等式是第t+1次迭代式M-Step的过程，最后一个等号就是t次迭代完的结果。至此，关于EM算法的相关推导证明全部结束。