EM算法的推导

算法用途

在统计计算中，最大期望（EM）算法是在 概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐藏变量（Latent Variable）。

如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法或者贝叶斯估计法来估计模型参数。当概率模型含有隐藏变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法。

符号定义

• Y：观测的随机变量数据，
• Z：隐藏的随机变量数据
• 完全数据（complete data）：Y和Z
• 不完全数据（incomplete data）：Y

问题描述

• 输入：观测数据 $Y$
• 求解：模型参数 $\theta$
• 方法：最大化观测数据的似然函数
  $$\theta^*=\arg_\theta {\max{L(\theta;Y)} }$$

其中，似然函数定义为：

$$L(\theta;Y) = \log P(Y;\theta) = \log \sum_{Z}{P(Y|Z;\theta)\cdot P(Z|\theta)}$$

• 困难：上式往往包含对未观测数据及和或者积分的对数，无解析解，不可直接求。

我们的解决方案是：通过迭代，不断优化似然函数 $L(\theta;Y)$ 的下界，来逐步近似极大化 $L(\theta;Y)$ 。这一下界，就是著名的Q函数。

Q函数

$$\begin{align*} Q(\theta, \theta^i) & = E_Z[\log P(Y,Z\ |\ \theta)\ |\ Y, \theta^i]\\ &= \sum_Z{ P(Z\ |\ Y,\theta^i)} \cdot \log P(Y,Z\ |\ \theta) \end{align*}$$

Q函数的物理意义为：

完全数据 ${Y, Z}$ 的 对数似然函数 $\log P(Y,Z\ |\ \theta)$，关于未观测数据 $Z$ 的 在给定 观测数据 $Y$ 和 当前参数估计 $\theta^i$ 的情况下的 后验概率分布 $P(Z\ |\ Y,\theta^i)$ 的条件期望。

EM算法的核心就在于Q函数的最大化。

极大似然估计的求解方法是，最大化观测变量关于模型参数的似然函数；而在EM方法的试用场景中含有隐藏变量，无法直接最大化观测变量的似然函数，因此我们转而通过迭代来不断优化观测变量似然函数的下界，即包括观测变量和隐藏变量在内的所有变量的似然函数的条件期望。

推导思路一

假设第 $i$ 轮迭代后，$\theta$ 的估计值是 $\theta^i$。

我们希望在新的一轮迭代中（即第 $i+1$ 轮），对 $\theta$ 的新的估计能够使得 $Y$ 的似然函数值增加，即

$$L(\theta;Y) \geq L(\theta^i;Y)$$

于是，我们考虑这二者的差：

$$L(\theta;Y) - L(\theta^i;Y)$$

我们注意到：

观测数据的似然函数为：

$$L(\theta;Y) = \log P(Y | \theta)$$

完全数据的似然函数为：

$$\begin{align*} L_C(\theta;Y, Z) &= \log P(Y,Z | \theta) \\ &= \log P(Y | \theta)\cdot P(Z| Y; \theta) \\ &= \log P(Y | \theta) +\log P(Z| Y; \theta) \\ &= L(\theta;Y)+\log P(Z| Y; \theta) \end{align*}$$

也就是说：

$$L(\theta;Y) = L_C(\theta;Y, Z) - \log P(Z| Y; \theta)$$

所以有

$$\begin{align*} L(\theta) - L(\theta^i) &= (L_C(\theta) - \log P(Z| Y; \theta)) - (L_C(\theta^i) - \log P(Z| Y; \theta^i))\\ &= L_C(\theta) - L_C(\theta^i) + \log \frac{P(Z| Y; \theta)}{P(Z| Y; \theta^i)} \end{align*}$$

对该式左右两边取关于Z的期望（$Y, \theta, i$ 已给定）。左式两项与Z无关，所以结果不变；右式三项均与Z相关。结果如下：

$$\begin{align*} L(\theta) - L(\theta^i) &= \sum_Z{L_C(\theta) \cdot P(Z|Y,\theta^i)} - \sum_Z{L_C(\theta^i)\cdot P(Z|Y,\theta^i)} + \sum_Z{ P(Z|Y,\theta^i)\cdot \log \frac{P(Z| Y; \theta)}{P(Z| Y; \theta^i)}}\\ &\geq \sum_Z{L_C(\theta) \cdot P(Z|Y,\theta^i)} - \sum_Z{L_C(\theta^i)\cdot P(Z|Y,\theta^i)} \end{align*}$$

不等号成立的原因是右式第三项实质上是两个概率分布 $P(Z| Y; \theta)$ 和 $P(Z| Y; \theta^i)$ 的KL-divergence，不为负。

上式又可转化为

$$L(\theta) \geq L(\theta^i) +\sum_Z{L_C(\theta) \cdot P(Z|Y,\theta^i)} - \sum_Z{L_C(\theta^i)\cdot P(Z|Y,\theta^i)}$$

此时，我们得到了 $L(\theta)$ 的下界，即右式所示。我们将该下界记为函数 $B(\theta, \theta^i)$。任何可以使 $B(\theta, \theta^i)$ 增大的 $\theta$，也可以使 $L(\theta)$ 增大。为了使 $L(\theta)$ 尽可能地变大，我们选择可以最大化函数 $B(\theta, \theta^i)$ 的 $\theta$ 作为对第 $i+1$ 轮的估计，即：

$$\begin{align*} \theta^{i+1} &= \arg_{\theta} \max B(\theta, \theta^i) \\ &= \arg_{\theta} \max L(\theta^i) +\sum_Z{L_C(\theta) \cdot P(Z|Y,\theta^i)} - \sum_Z{L_C(\theta^i)\cdot P(Z|Y,\theta^i)} \end{align*}$$

我们注意到：右式第一项和第三项均为常数，因此，上述问题等价于：

$$\begin{align*} \theta^{i+1} &= \arg_{\theta} \sum_Z{L_C(\theta) \cdot P(Z|Y,\theta^i)} \\ &= \arg_{\theta} \sum_Z{P(Z|Y,\theta^i) \log P(Y,Z\ |\ \theta)} \\ &= Q(\theta, \theta^i) \end{align*}$$

这一下界就是EM算法中的Q函数：

$$\begin{align*} Q(\theta, \theta^i) &= \sum_Z{L_C(\theta) \cdot P(Z\ |\ Y,\theta^i)} \\ &= \sum_Z{ P(Z\ |\ Y,\theta^i)} \cdot \log P(Y,Z\ |\ \theta)\\ & = E_Z[\log P(Y,Z\ |\ \theta)\ |\ Y, \theta^i] \end{align*}$$

推导思路二

同样，假设第 $i$ 轮迭代后，$\theta$ 的估计值是 $\theta^i$。

我们希望在新的一轮迭代中（即第 $i+1$ 轮），对 $\theta$ 的新的估计能够使得 $Y$ 的似然函数值增加，即

$$L(\theta;Y) \geq L(\theta^i;Y)$$

于是，我们考虑这二者的差：

$$\begin{align*} L(\theta;Y) - L(\theta^i;Y) &= \log P(Y\ |\ \theta) - \log P(Y\ |\ \theta^i)\\ &= \log \sum_Z P(Y,Z\ |\ \theta) - \log P(Y\ |\ \theta^i)\\ &= \log \sum_Z P(Z\ |\ Y,\theta^i)\frac{P(Y,Z\ |\ \theta)}{P(Z\ |\ Y,\theta^i)} - \log P(Y\ |\ \theta^i)\\ &\geq \sum_Z P(Z\ |\ Y,\theta^i)\log\frac{P(Y,Z\ |\ \theta)}{P(Z\ |\ Y,\theta^i)} - \log P(Y\ |\ \theta^i)\\ &= \sum_Z P(Z\ |\ Y,\theta^i)\log\frac{P(Y,Z\ |\ \theta)}{P(Z\ |\ Y,\theta^i)} -\sum_Z P(Z\ |\ Y,\theta^i)\log P(Y\ |\ \theta^i)\\ &= \sum_Z P(Z\ |\ Y,\theta^i)\log\frac{P(Y,Z\ |\ \theta)}{P(Z\ |\ Y,\theta^i)P(Y\ |\ \theta^i)}\\ &= \sum_Z P(Z\ |\ Y,\theta^i)\log\frac{P(Y,Z\ |\ \theta)}{P(Y,Z\ |\ \theta^i)}\\ \end{align*}$$

上式不等号是利用了Jensen不等式：

对于凸函数 $f(x)$：

$$f\Big(\sum_{\lambda_i}x_i\Big) \geq \sum_{\lambda_i}f\Big(x_i\Big)$$

其中 $\lambda_i\geq 0, \sum \lambda_i = 1 .$

而对数函数 $log(\cdot)$ （以 $e$ 为底）为凹函数。

将上式的形式进一步转换可以得到：

$$\begin{align*} L(\theta;Y) &= L(\theta^i;Y) + \sum_Z P(Z\ |\ Y,\theta^i)\log\frac{P(Y,Z\ |\ \theta)}{P(Y,Z\ |\ \theta^i)}\\ &= L(\theta^i;Y) + \sum_Z P(Z\ |\ Y,\theta^i)\log P(Y,Z\ |\ \theta) - \sum_Z P(Z\ |\ Y,\theta^i) P(Y,Z\ |\ \theta^i)\\ \end{align*}$$

此时，我们得到了 $L(\theta)$ 的下界，和第一种推导思路一致。

EM算法框架

至此，我们可以给出EM算法框架如下：

输入：观测数据 Y
设定：模型及其参数 $\theta$，以及$P(Y\ |\ Z)$ 和 $P(Z\ |\ \theta)$
求解：模型参数 $\theta$

(1) 初始化：选择参数初值 $\theta_0$，开始迭代

(2) E-step：记 $\theta_i$ 为第 $i$ 次迭代的参数值，则在第 $i+1$ 次迭代的E-step，计算Q函数，即完全数据在当前参数下的条件期望：

$$\begin{align*} Q(\theta, \theta^i) & = E_Z[\log P(Y,Z\ |\ \theta)\ |\ Y, \theta^i]\\ &= \sum_Z{ P(Z\ |\ Y,\theta^i)} \cdot \log P(Y,Z\ |\ \theta) \end{align*}$$

其中Q函数的计算涉及两项概率分布：

• 完全数据 ${Y, Z}$ 的 对数似然函数 $\log P(Y,Z\ |\ \theta)$：
  $$\begin{align*} \log P(Y,Z\ |\ \theta) &= \log P(Y\ |\ Z) \ P(Z\ |\ \theta) \end{align*}$$
• 未观测数据 $Z$ 的 后验概率分布 $P(Z\ |\ Y,\theta^i)$：
  $$\begin{align*} P(Z\ |\ Y,\theta^i) = \frac{P(Y,Z\ | \theta^i)}{P(Y\ |\ \theta^i)} =\frac{P(Y,Z\ | \theta^i)}{\sum_{Z'}P(Y,Z'\ |\ \theta^i)} =\frac{P(Y\ |\ Z) \ P(Z\ | \theta^i)}{\sum_{Z'}P(Y\ |\ Z')\ P(Z'\ |\ \theta^i)} \end{align*}$$

(3) M-step：求使得上述期望最大的 $\theta$，确定第 $i+1$ 次迭代的参数估计值 $\theta_{i+1}$：

$$\theta_{i+1} = \arg_{\theta} \max Q(\theta, \theta^i)$$

(4) 重复第2、3步，直到收敛。