EM算法

机器学习

---

常见使用场景

机器学习中，频率派经常会通过最大似然函数来求模型中的参数，但是有些模型中包含有不可见的隐变量Z，这样再通过最大似然求解就非常困难。一般如果模型中包含隐变量，我们就会使用加强版的最大似然-EM算法来求参数。常见包含隐参数的模型有高斯混合模型，隐马尔科夫模型，PLSA模型等等。下面我们以高斯混合模型为例来阐述为何有隐变量我们无法使用最大似然来估计参数。
首先我们先看看GMM的对数似然函数，单个高斯模型的参数我们用$\theta=(u,\Sigma)$表示，而在GMM模型中$\pi$是隐变量：

$$lnp(X|\pi,\theta)=\sum_{n=1}^Nln\sum_{\pi} p(x_n,\pi|\theta)=\sum_{n=1}^Nln\sum_{k=1}^K \pi_kN(x_n|u_k,\Sigma_k)$$
由于对数里面有加和，直接对$p(X|\theta)$求最值很困难，我们可以试想下，如果隐变量$\pi$是可以观测的，也就是知道每个数据点是由哪几个分布生成的，那么我们求解就会很方便，但关键的是$\pi$是隐变量，我们没法观测到。但是我们可以用它的期望来表示。因此，只要将隐变量用它的期望表示就得到了大名鼎鼎的EM算法。具体过程如下:
$lnp(X|\theta)$是我们的原始似然函数，加入隐变量Z后可以写成$ln\sum_Zp(X,Z|\theta)$，我们先初始化模型的参数$\theta_{old}$，由于隐变量无法观测到，我们用参数来得到它的后验$p(Z|X,\theta_{old})$，然后我们通过隐藏变量的期望得到新的complete-data likelihood：
$$Q(\theta,\theta_{old})=\int lnp(X,Z|\theta)p(Z|X,\theta_{old})dZ$$
以上就是E 步，当我们得到期望后就可以求出Q函数的最优解，也就是新的参数$\theta=argmaxQ(\theta,\theta_{old})$，注意这个Q函数是包含隐变量的complete-data的似然函数，并不是我们一开始的似然函数，所以通过complete-data likelihood是可解析出新的参数。当我们将新的参数$\theta$ 再带入Q 函数不断就开始了第二轮迭代。

EM算法

下面我们将来看看EM算法是怎么推导的：
先回顾优化理论中的一些概念。设f是定义域为实数的函数，如果对于所有的实数x，$f''(x)\geq 0$，那么f是凸函数。当x是向量时，如果其Hessian矩阵H是半正定的$H\succeq 0$，那么f是凸函数。如果$f''(x)>0$或者$H\succ 0$，那么称f是严格凸函数。

如果f是凸函数， 通过Jensen不等式我们知道：

$$f(\sum_{i=1}^n \lambda_ix_i)\leq \sum_{i=1}^n \lambda_if(x_i)\quad \lambda_i\geq 0\quad\sum_{i=1}^n\lambda_i=1$$
如果将$x_i$看成随机变量，而将$\lambda_i$认为是$x_i$的出现概率则有：
$$f(E(x)) \leq E[f(x)]$$
特别地，如果f是严格凸函数，那么当且仅当$x=E[x]$恒成立时，也就是说X是常量，有
$$E[f(x)]=f(E(x))$$
如果用图表示会很清晰：

图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到

$$f(E(x)) \leq E[f(x)]$$ 成立

回到正题，我们现在已经知道了Jensen不等式，假设给定的训练样本是${x_1,x_2,...,x_n}$，样例间独立，我们想找到每个样例隐含的类别$z_i$，能使得似然估计最大：

$$\begin{align} \ell(\theta) &= \sum_{i=1}^mlogp(x|\theta)\\ &= \sum_{i=1}^mlog\sum_zp(x,z|\theta) \end{align}$$

第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求$\theta$一般比较困难，因为有隐变量z存在，但是一般确定了z后，求解就容易了。EM是一种解决存在隐变量优化问题的有效方法。既然不能直接最大化$\ell(\theta)$，我们可以不断地建立对数似然$\ell(\theta)$的下界（E步），然后优化下界（M步）。这句话比较抽象，继续看下面：

对于每一个样例i，让$Q_i$表示该样例隐变量z的某种分布，$Q_i$满足的条件是$\sum_z Q_i(z)=1,Q_i(z)>0$。（如果z是连续性的，那么$Q_i$是概率密度函数，需要将求和符号换做积分符号）。比如要将班上学生聚类，假设隐藏变量z是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布了。可以由前面阐述的内容得到下面的公式：

$$\begin{align} \ell(\theta) &= \sum_{i=1}^mlog\sum_zp(x,z|\theta)\\ &= \sum_{i=1}^mlog\sum_zQ_i \frac{p(x,z|\theta)}{Q_i}\\ &\geq \sum_{i=1}^m\sum_z Q_i log \frac{p(x,z|\theta)}{Q_i} \end{align}$$
上式中第2步到第3步利用了Jensen不等式，考虑到$log$是凹函数（二阶导数小于0），而且$\sum_zQ_i \frac{p(x,z|\theta)}{Q_i}$是函数$\frac{p(x,z|\theta)}{Q_i}$在$Q_i$分布下的期望，所以$logE[ \frac{p(x,z|\theta)}{Q_i}]\geq E[log \frac{p(x,z|\theta)}{Q_i}]$。

首先我们令$\ell'(\theta)=\sum_{i=1}^m\sum_z Q_i log \frac{p(x,z|\theta)}{Q_i}$。上述过程可以看作是$\ell'(\theta)$是$\ell(\theta)$的下界。对于$Q_i$的选择，有多种可能，哪种更好呢？假设$\theta$已经给定，那么$\ell'(\theta)$的值就决定于$Q_i$和$p(x,z)$了。我们可以通过调整这两个概率使下界不断上升，以逼近$\ell(\theta)$的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够使$\ell'(\theta)=\ell(\theta)$了。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：

$$\frac{p(x,z|\theta)}{Q_i}=c \qquad (1)$$
c为常数，不依赖于$z_i$。对此式子做进一步推导，我们知道$\sum_z Q_i(z)=1$，那么也就有$\sum_z p(x,z|\theta)=c$（多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），将$\sum_z p(x,z|\theta)=c$带入(1)式：
$$\begin{align} Q_i(z) &= \frac{p(x,z|\theta)}{c}\\ &= \frac{p(x,z|\theta)}{\sum p(x,z|\theta)}\\ &= p(z|x,\theta) \qquad (2) \end{align}$$
至此，我们推出了在固定参数$\theta$后，$Q_i(z)$的计算公式就是后验概率$p(z|x,\theta)$，解决了$Q_i(z)$如何选择的问题。这一步就是E步，建立$\ell(\theta)$的下界。接下来的M步，就是在给定$Q_i(z)$后，调整参数$\theta$，取极大化$\ell(\theta)$的下界（在固定$Q_i(z)$后，调整$\theta$，使下界还可以调整的更大）。

上图中$F(\theta)$就是原始的极大似然函数$\ell(\theta)$，而$G_{\theta_t}(\theta_t)$和$G_{\theta_t}(\theta_{t+1})$分别是第t次和第t+1次迭代的下界函数$\ell'(\theta)$，从图中可以看出通过上次一次迭代或者最开始的初始化一个参数$\theta_t$，在E步我们可以求出在参数$G_{\theta_t}(\theta_t)$与原始极大似然函数$F(\theta)$想等的下界函数$G_{\theta_t}(\theta_t)$，然后在M步我们下界函数求的极大值，也就是下次迭代的参数$\theta_{t+1}$。
将(2)式带入$\ell'(\theta)$继续推导至简单形式：
$$\begin{align} \theta_{t+1} &= argmax_\theta \ell'(\theta)\\ &= argmax_\theta \sum_{i=1}^m\sum_z Q_i log \frac{p(x,z|\theta)}{Q_i}\\ &= argmax_\theta \sum_{i=1}^m\sum_zp(z|x,\theta_{t}) log \frac{p(x,z|\theta)}{p(z|x,\theta_{t})}\\ &= argmax_\theta \sum_{i=1}^m\sum_zp(z|x,\theta_{t})logp(x,z|\theta)-p(z|x,\theta_{t})logp(z|x,\theta_{t}) \end{align}$$

因为$p(z|x,\theta_{t})logp(z|x,\theta_{t})$中的$\theta_{t}$都是定值，所以对$\theta$求偏导时，会消掉，因此M步最大化的只有第一项$sum_z p(z|x,\theta_{t})logp(x,z|\theta)$，即complete-data关于隐变量z的期望$E_{z|x,\theta_t}[logp(x,z|\theta)]$。
那么一般的EM算法的步骤如下：
1. 初始化$\theta_{old}$
2. E步：对于每一个$i$，计算$p(z_i|x_i,\theta)$，求complete-data liklihood $Q(\theta,\theta_{old})=\sum_z logp(x,z|\theta)p(z|x,\theta_{old})dz$
3. M步：更新参数$\theta=argmax_\theta Q(\theta,\theta_{old})$
4. 跳转到第2步

应用案例

1.GMM

参考资料

1. http://blog.csdn.net/abcjennifer/article/details/8170378
2. http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html
3. http://blog.csdn.net/zouxy09/article/details/8537620
4. PRML, chapter 9
5. Notes on Pattern Recognition and Machine Learning (Jian Xiao)
6. Pattern Recognition And Machine Learning 读书会, chapter 9