一、概述

介绍

概率模型有时既包含观测变量（observed variable），又包含隐变量（latent variable）。当概率模型只包含观测变量时，那么给定观测数据，就可以直接使用极大似然估计法或者贝叶斯估计法进行模型参数的求解。然而如果模型包含隐变量，就不能直接使用这些简单的方法了。EM算法就是用来解决这种含有隐变量的概率模型参数的极大似然参数估计法。这里只讨论极大似然估计，极大后验估计与其类似。

算法

EM算法的输入如下：

$X$ :观测数据
$Z$ :未观测数据（隐变量）
$p(x,z|\theta )$ :联合分布
$p(z|x,\theta )$ :后验分布
$\theta$ :parameter

在算法运行开始时需要选择模型的初始化参数 $\theta ^{(0)}$ 。EM算法是一种迭代更新的算法，其计算公式为：

$\theta ^{t+1}=\underset{\theta }{argmax}E_{z|x,\theta^{t}}[log\; p(x,z|\theta )]\\ =\underset{\theta }{argmax}\int _{z}log\; p(x,z|\theta )\cdot p(z|x,\theta ^{t})\mathrm{d}z$

这个公式包含了迭代的两步：
①E step：计算 $p(x,z|\theta )$ 在概率分布 $p(z|x,\theta ^{t})$ 下的期望；
②M step：计算使这个期望最大化的参数得到下一个EM步骤的输入。

总结来说，EM算法包含以下步骤：
①选择初始化参数 $\theta ^{(0)}$ ；
②E step；
③M step；
④重复②③步直至收敛。

二、EM算法的收敛性

现在要证明迭代求得的 $\theta ^{t}$ 序列会使得对应的 $p(x|\theta ^{t})$ 是单调递增的（如果 $p(x|\theta ^{t})$ 是单调递增的，那么训练数据的似然就是单调递增的），也就是说要证明 $p(x|\theta ^{t})\leq p(x|\theta ^{t+1})$ 。首先我们有：

$log\; p(x|\theta )=log\; p(x,z|\theta )-log\; p(z|x,\theta )$

接下来等式两边同时求关于 $p(z|x,\theta ^{t})$ 的期望：

$左边=\int _{z}p(z|x,\theta ^{t})\cdot log\; p(x|\theta )\mathrm{d}z\\ =log\; p(x|\theta )\int _{z}p(z|x,\theta ^{t})\mathrm{d}z\\ =log\; p(x|\theta )\\ 右边=\underset{记作Q(\theta ,\theta ^{t})}{\underbrace{\int _{z}p(z|x,\theta ^{t})\cdot p(x,z|\theta )\mathrm{d}z}}-\underset{记作H(\theta ,\theta ^{t})}{\underbrace{\int _{z}p(z|x,\theta ^{t})\cdot log\; p(z|x,\theta )\mathrm{d}z}}$

因此有：

$log\; p(x|\theta )=\int _{z}p(z|x,\theta ^{t})\cdot p(x,z|\theta )\mathrm{d}z-\int _{z}p(z|x,\theta ^{t})\cdot log\; p(z|x,\theta )\mathrm{d}z$

这里定义了 $Q(\theta ,\theta ^{t})$ ，称为Q函数（Q function），这个函数也就是上面的概述中迭代公式里用到的函数，因此满足 $Q(\theta ^{t+1},\theta ^{t})\geq Q(\theta ^{t},\theta ^{t})$ 。

接下来将上面的等式两边 $\theta$ 分别取 $\theta ^{t+1}$ 和 $\theta ^{t}$ 并相减：

$log\; p(x|\theta ^{t+1})-log\; p(x|\theta ^{t})=[Q(\theta ^{t+1},\theta ^{t})-Q(\theta ^{t},\theta ^{t})]-[H(\theta ^{t+1},\theta ^{t})-H(\theta ^{t},\theta ^{t})]$

我们需要证明 $log\; p(x|\theta ^{t+1})-log\; p(x|\theta ^{t})\geq 0$ ，同时已知 $Q(\theta ^{t+1},\theta ^{t})-Q(\theta ^{t},\theta ^{t})\geq 0$ ，现在来观察 $H(\theta ^{t+1},\theta ^{t})-H(\theta ^{t},\theta ^{t})$ ：

$H(\theta ^{t+1},\theta ^{t})-H(\theta ^{t},\theta ^{t})\\ =\int _{z}p(z|x,\theta ^{t})\cdot log\; p(z|x,\theta ^{t+1})\mathrm{d}z-\int _{z}p(z|x,\theta ^{t})\cdot log\; p(z|x,\theta ^{t})\mathrm{d}z\\ =\int _{z}p(z|x,\theta ^{t})\cdot log\frac{p(z|x,\theta ^{t+1})}{p(z|x,\theta ^{t})}\mathrm{d}z\\ \leq log\int _{z}p(z|x,\theta ^{t})\frac{p(z|x,\theta ^{t+1})}{p(z|x,\theta ^{t})}\mathrm{d}z\\ =log\int _{z}p(z|x,\theta ^{t+1})\mathrm{d}z\\ =log\; 1\\ =0$

这里的不等号应用了Jensen不等式：

$log\sum _{j}\lambda _{j}y_{j}\geq \sum _{j}\lambda _{j}log\; y_{j},其中\lambda _{j}\geq 0，\sum _{j}\lambda _{j}=1$

也可以使用KL散度来证明 $\int _{z}p(z|x,\theta ^{t})\cdot log\frac{p(z|x,\theta ^{t+1})}{p(z|x,\theta ^{t})}\mathrm{d}z\leq 0$ ，两个概率分布 $P(x)$ 和 $Q(x)$ 的KL散度是恒 $\geq 0$ 的，定义为：

$D_{KL}(P||Q)=E_{x\sim P}[log\frac{P(x)}{Q(x)}]$

因此有：

$\int _{z}p(z|x,\theta ^{t})\cdot log\frac {p(z|x,\theta ^{t+1})}{p(z|x,\theta ^{t})}\mathrm{d}z=-KL(p(z|x,\theta ^{t})||p(z|x,\theta ^{t+1}))\leq 0$

因此得证 $log\; p(x|\theta ^{t+1})-log\; p(x|\theta ^{t})\geq 0$ 。这说明使用EM算法迭代更新参数可以使得 $log\; p(x|\theta)$ 逐步增大。

另外还有其他定理保证了EM的算法收敛性。首先对于 $\theta ^{i}(i=1,2,\cdots )$ 序列和其对应的对数似然序列 $L(\theta ^{t})=log\; p(x|\theta ^{t})(t=1,2,\cdots )$ 有如下定理：
①如果 $p(x|\theta )$ 有上界，则 $L(\theta ^{t})=log\; p(x|\theta ^{t})$ 收敛到某一值 $L^*$ ；
②在函数 $Q(\theta,\theta^{'})$ 与 $L(\theta )$ 满足一定条件下，由EM算法得到的参数估计序列 $\theta ^{t}$ 的收敛值 $\theta ^{*}$ 是 $L(\theta )$ 的稳定点。

三、EM算法的导出

ELBO+KL散度的方法

对于前面用过的式子，首先引入一个新的概率分布 $q(z)$ ：

$log\; p(x|\theta )=log\; p(x,z|\theta )-log\; p(z|x,\theta )\\ =log\; \frac {p(x,z|\theta )}{q(z)}-log\; \frac{p(z|x,\theta )}{q(z)}\; \; q(z)\neq 0$

以上引入一个关于 $z$ 的概率分布 $q(z)$ ，然后式子两边同时求对 $q(z)$ 的期望：

$左边=\int _{z}q(z)\cdot log\; p(x|\theta )\mathrm{d}z=log\; p(x|\theta )\int _{z}q(z)\mathrm{d}z=log\; p(x|\theta )\\ 右边=\underset{ELBO(evidence\; lower\; bound)}{\underbrace{\int _{z}q(z)log\; \frac{p(x,z|\theta )}{q(z)}\mathrm{d}z}}\underset{KL(q(z)||p(z|x,\theta ))}{\underbrace{-\int _{z}q(z)log\; \frac{p(z|x,\theta )}{q(z)}\mathrm{d}z}}$

因此我们得出 $log\; p(x|\theta )=ELBO+KL(q||p)$ ，由于KL散度恒 $\geq0$ ，因此 $log\; p(x|\theta )\geq ELBO$ ，则 $ELBO$ 就是似然函数 $log\; p(x|\theta )$ 的下界。使得 $log\; p(x|\theta )=ELBO$ 时，就必须有 $KL(q||p)=0$ ，也就是 $q(z)=p(z|x,\theta )$ 时。在每次迭代中我们取 $q(z)=p(z|x,\theta ^{t})$ ，就可以保证 $log\; p(x|\theta ^{t})$ 与 $ELBO$ 相等，也就是：

$log\; p(x|\theta )=\underset{ELBO}{\underbrace{\int _{z}p(z|x,\theta ^{t})log\; \frac {p(x,z|\theta )}{p(z|x,\theta ^{t})}\mathrm{d}z}}\underset{KL(p(z|x,\theta ^{t})||p(z|x,\theta ))}{\underbrace{-\int _{z}p(z|x,\theta ^{t})log\; \frac{p(z|x,\theta )}{p(z|x,\theta ^{t})}\mathrm{d}z}}$

当 $\theta =\theta ^{t}$ 时， $log\; p(x|\theta ^{t})$ 取ELBO，即：

$log\; p(x|\theta ^{t})=\underset{ELBO}{\underbrace{\int _{z}p(z|x,\theta ^{t})log\; \frac{p(x,z|\theta ^{t})}{p(z|x,\theta ^{t})}\mathrm{d}z}}\underset{=0}{\underbrace{-\int _{z}p(z|x,\theta ^{t})log\; \frac{p(z|x,\theta ^{t})}{p(z|x,\theta ^{t})}\mathrm{d}z}}=ELBO$

也就是说 $log\; p(x|\theta )$ 与 $ELBO$ 都是关于 $\theta$ 的函数，且满足 $log\; p(x|\theta )\geq ELBO$ ，也就是说 $log\; p(x|\theta )$ 的图像总是在 $ELBO$ 的图像的上面。对于 $q(z)$ ，我们取 $q(z)=p(z|x,\theta ^{t})$ ，这也就保证了只有在 $\theta =\theta ^t$ 时 $log\; p(x|\theta )$ 与 $ELBO$ 才会相等，因此使 $ELBO$ 取极大值的 $\theta ^{t+1}$ 一定能使得 $log\; p(x|\theta ^{t+1})\geq log\; p(x|\theta ^{t})$ 。该过程如下图所示：

ELBO

然后我们观察一下 $ELBO$ 取极大值的过程：

$\theta ^{t+1}=\underset{\theta }{argmax}ELBO \\ =\underset{\theta }{argmax}\int _{z}p(z|x,\theta ^{t})log\; \frac{p(x,z|\theta )}{p(z|x,\theta ^{t})}\mathrm{d}z\\ =\underset{\theta }{argmax}\int _{z}p(z|x,\theta ^{t})log\; p(x,z|\theta )\mathrm{d}z-\underset{与\theta 无关}{\underbrace{\underset{\theta }{argmax}\int _{z}p(z|x,\theta ^{t})p(z|x,\theta ^{t})\mathrm{d}z}}\\ {\color{Red}{=\underset{\theta }{argmax}\int _{z}p(z|x,\theta ^{t})log\; p(x,z|\theta )\mathrm{d}z}} \\ {\color{Red}{=\underset{\theta }{argmax}E_{z|x,\theta ^{t}}[log\; p(x,z|\theta )]}}$

由此我们就导出了EM算法的迭代公式。

ELBO+Jensen不等式的方法

首先要具体介绍一下Jensen不等式：对于一个凹函数 $f(x)$ （国内外对凹凸函数的定义恰好相反，这里的凹函数指的是国外定义的凹函数），我们查看其图像如下：

Jensen不等式

$t\in [0,1]\\ c=ta+(1-t)b\\ \phi =tf(a)+(1-t)f(b)$

凹函数恒有 $f(c)\geq \phi$ ，也就是 $f(ta+(1-t)b)\geq tf(a)+(1-t)f(b)$ ，当 $t=\frac{1}{2}$ 时有 $f(\frac{a}{2}+\frac{b}{2})\geq \frac{f(a)}{2}+\frac{f(b)}{2}$ ，可以理解为对于凹函数来说 先求期望再求函数值 恒 $\geq$ 先求函数值再求期望 ，即 $f(E)\geq E[f]$ 。

上面的说明只是对Jensen不等式的一个形象的描述，而非严谨的证明。接下来应用Jensen不等式来导出EM算法：

$log\; p(x|\theta )=log\int _{z}p(x,z|\theta )\mathrm{d}z\\ =log\int _{z}\frac{p(x,z|\theta )}{q(z)}\cdot q(z)\mathrm{d}z\\ =log\; E_{q(z)}[\frac{p(x,z|\theta )}{q(z)}]\\ \geq \underset{ELBO}{\underbrace{E_{q(z)}[log\frac{p(x,z|\theta )}{q(z)}]}}$

这里应用了Jensen不等式得到了上面出现过的 $ELBO$ ，这里的 $f(x)$ 函数也就是 $log$ 函数，显然这是一个凹函数。当 $log\frac {P(x,z|\theta )}{q(z)}$ 这个函数是一个常数时会取得等号，利用这一点我们也同样可以得到 $q(z)=p(z|x,\theta )$ 时能够使得 $log\; p(x|\theta )=ELBO$ 的结论：

$\frac{p(x,z|\theta )}{q(z)}=C\\ \Rightarrow q(z)=\frac{p(x,z|\theta )}{C}\\ \Rightarrow \int _{z}q(z)\mathrm{d}z=\int _{z}\frac{1}{C}p(x,z|\theta )\mathrm{d}z\\ \Rightarrow 1=\frac{1}{C}\int _{z}p(x,z|\theta )\mathrm{d}z\\ \Rightarrow C=p(x|\theta )\\ 将C代入q(z)=\frac{p(x,z|\theta )}{C}得\\ {\color{Red}{q(z)=\frac{p(x,z|\theta )}{p(x|\theta )}=p(z|x,\theta )}}$

这种方法到这里就和上面的方法一样了，总结来说就是：

$log\; p(x|\theta )\geq \underset{ELBO}{\underbrace{E_{q(z)}[log\frac{p(x,z|\theta )}{q(z)}]}}$
上面的不等式在 $q(z)=p(z|x|\theta )$ 时取等号，因此在迭代更新过程中取 $q(z)=p(z|x,\theta ^{t})$ 。接下来的推导过程就和第1种方法一样了。

四、广义EM算法

上面介绍的EM算法属于狭义的EM算法，它是广义EM的一个特例。在上面介绍的EM算法的E步中我们假定 $q(z)=p(z|x,\theta ^{t})$ ，但是如果这个后验 $p(z|x,\theta ^{t})$ 无法求解，那么必须使⽤采样（MCMC）或者变分推断等⽅法来近似推断这个后验。前面我们得出了以下关系：

$log\; p(x|\theta )=\int _{z}q(z)log\; \frac{p(x,z|\theta )}{q(z)}\mathrm{d}z-\int _{z}q(z)log\; \frac{p(z|x,\theta )}{q(z)}\mathrm{d}z=ELBO+KL(q||p)$