EM算法|机器学习推导系列（十二）

一、概述

假设有如下数据：

$X$X:observed data
$Z$Z:unobserved data(latent variable)
$(X,Z)$(X,Z):complete data
$\theta$θ:parameter

EM算法的目的是解决具有隐变量的参数估计（MLE、MAP）问题。EM算法是一种迭代更新的算法，其计算公式为：

$\theta ^{t+1}=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$θt+1=Ez∣x,θt​[logp(x,z∣θ)]=θargmax​∫z​logp(x,z∣θ)⋅p(z∣x,θt)dz

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

二、EM的算法收敛性

现在要证明迭代求得的$\theta ^{t}$θt序列会使得对应的$p(x|\theta ^{t})$p(x∣θt)是单调递增的，也就是说要证明$p(x|\theta ^{t})\leq p(x|\theta ^{t+1})$p(x∣θt)≤p(x∣θt+1)。首先我们有：

$log\; p(x|\theta )=log\; p(x,z|\theta )-log\; p(z|x,\theta )$logp(x∣θ)=logp(x,z∣θ)−logp(z∣x,θ)

接下来等式两边同时求关于$p(z|x,\theta ^{t})$p(z∣x,θ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$左边=∫z​p(z∣x,θt)⋅logp(x∣θ)dz=logp(x∣θ)∫z​p(z∣x,θt)dz=logp(x∣θ)右边=Q(θ,θt)∫z​p(z∣x,θt)⋅p(x,z∣θ)dz​​−H(θ,θt)∫z​p(z∣x,θt)⋅logp(z∣x,θ)dz​​因此有logp(x∣θ)=∫z​p(z∣x,θt)⋅p(x,z∣θ)dz−∫z​p(z∣x,θt)⋅logp(z∣x,θ)dz

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

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

$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})]$logp(x∣θt+1)−logp(x∣θt)=[Q(θt+1,θt)−Q(θt,θt)]−[H(θt+1,θt)−H(θt,θt)]

我们需要证明$log\; p(x|\theta ^{t+1})-log\; p(x|\theta ^{t})\geq 0$logp(x∣θt+1)−logp(x∣θt)≥0，同时已知$Q(\theta ^{t+1},\theta ^{t})-Q(\theta ^{t},\theta ^{t})\geq 0$Q(θt+1,θt)−Q(θt,θt)≥0，现在来观察$H(\theta ^{t+1},\theta ^{t})-H(\theta ^{t},\theta ^{t})$H(θt+1,θt)−H(θt,θ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不等式：{\color{Red} {log\sum _{j}\lambda _{j}y_{j}\geq \sum _{j}\lambda _{j}logy_{j},其中\lambda _{j}\geq 0，log\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散度的定义为\\ D_{KL}(P||Q)=E_{x\sim P}[log\frac{P(x)}{Q(x)}]，KL散度是恒\geq 0的。\\ 因此\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$H(θt+1,θt)−H(θt,θt)=∫z​p(z∣x,θt)⋅logp(z∣x,θt+1)dz−∫z​p(z∣x,θt)⋅logp(z∣x,θt)dz=∫z​p(z∣x,θt)⋅logp(z∣x,θt)p(z∣x,θt+1)​dz≤log∫z​p(z∣x,θt)p(z∣x,θt)p(z∣x,θt+1)​dz=log∫z​p(z∣x,θt+1)dz=log1=0这里应用了Jensen不等式：logj∑​λj​yj​≥j∑​λj​logyj​,其中λj​≥0，logj∑​λj​=1也可以使用KL散度来证明∫z​p(z∣x,θt)⋅logp(z∣x,θt)p(z∣x,θt+1)​dz≤0：两个概率分布P(x)和Q(x)的KL散度的定义为DKL​(P∣∣Q)=Ex∼P​[logQ(x)P(x)​]，KL散度是恒≥0的。因此∫z​p(z∣x,θt)⋅logp(z∣x,θt)p(z∣x,θt+1)​dz=−KL(p(z∣x,θt)∣∣p(z∣x,θt+1))≤0

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

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

三、EM的算法的导出

1. ELBO+KL散度的方法

$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}}$logp(x∣θ)=logp(x,z∣θ)−logp(z∣x,θ)=logq(z)p(x,z∣θ)​−logq(z)p(z∣x,θ)​q(z)​=0以上引入一个关于z的概率分布q(z)，然后式子两边同时求对q(z)的期望左边=∫z​q(z)⋅logp(x∣θ)dz=logp(x∣θ)∫z​q(z)dz=logp(x∣θ)右边=ELBO(evidencelowerbound)∫z​q(z)logq(z)p(x,z∣θ)​dz​​KL(q(z)∣∣p(z∣x,θ))−∫z​q(z)logq(z)p(z∣x,θ)​dz​​

因此我们得出$log\; p(x|\theta )=ELBO+KL(q||p)$logp(x∣θ)=ELBO+KL(q∣∣p)，由于KL散度恒$\geq0$≥0，因此$log\; p(x|\theta )\geq ELBO$logp(x∣θ)≥ELBO，则$ELBO$ELBO就是似然函数$log\; p(x|\theta )$logp(x∣θ)的下界。

使得$log\; p(x|\theta )=ELBO$logp(x∣θ)=ELBO时，就必须有$KL(q||p)=0$KL(q∣∣p)=0，也就是$q(z)=p(z|x,\theta )$q(z)=p(z∣x,θ)时。

在每次迭代中我们取$q(z)=p(z|x,\theta ^{t})$q(z)=p(z∣x,θt)，就可以保证$log\; p(x|\theta ^{t})$logp(x∣θt)与$ELBO$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$logp(x∣θ)=ELBO∫z​p(z∣x,θt)logp(z∣x,θt)p(x,z∣θ)​dz​​KL(p(z∣x,θt)∣∣p(z∣x,θ))−∫z​p(z∣x,θt)logp(z∣x,θt)p(z∣x,θ)​dz​​当θ=θt时，logp(x∣θt)取ELBO，即logp(x∣θt)=ELBO∫z​p(z∣x,θt)logp(z∣x,θt)p(x,z∣θt)​dz​​=0−∫z​p(z∣x,θt)logp(z∣x,θt)p(z∣x,θt)​dz​​=ELBO

也就是说$log\; p(x|\theta )$logp(x∣θ)与$ELBO$ELBO都是关于$\theta$θ的函数，且满足$log\; p(x|\theta )\geq ELBO$logp(x∣θ)≥ELBO，也就是说$log\; p(x|\theta )$logp(x∣θ)的图像总是在$ELBO$ELBO的图像的上面。对于$q(z)$q(z)，我们取$q(z)=p(z|x,\theta ^{t})$q(z)=p(z∣x,θt)，这也就保证了只有在$\theta =\theta ^t$θ=θt时$log\; p(x|\theta )$logp(x∣θ)与$ELBO$ELBO才会相等，因此使$ELBO$ELBO取极大值的$\theta ^{t+1}$θt+1一定能使得$log\; p(x|\theta ^{t+1})\geq log\; p(x|\theta ^{t})$logp(x∣θt+1)≥logp(x∣θ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 )]}}$θt+1=θargmax​ELBO=θargmax​∫z​p(z∣x,θt)logp(z∣x,θt)p(x,z∣θ)​dz=θargmax​∫z​p(z∣x,θt)logp(x,z∣θ)dz−与θ无关θargmax​∫z​p(z∣x,θt)p(z∣x,θt)dz​​=θargmax​∫z​p(z∣x,θt)logp(x,z∣θ)dz=θargmax​Ez∣x,θt​[logp(x,z∣θ)]

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

2. ELBO+Jensen不等式的方法

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

$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]$t∈[0,1]c=ta+(1−t)bϕ=tf(a)+(1−t)f(b)凹函数恒有f(c)≥ϕ也就是f(ta+(1−t)b)≥tf(a)+(1−t)f(b)当t=21​时有f(2a​+2b​)≥2f(a)​+2f(b)​可以理解为先求期望再求函数恒≥先求函数再求期望即f(E)≥E[f]

接下来应用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)}]}}$logp(x∣θ)=log∫z​p(x,z∣θ)dz=log∫z​q(z)p(x,z∣θ)​⋅q(z)dz=logEq(z)​[q(z)p(x,z∣θ)​]≥ELBOEq(z)​[logq(z)p(x,z∣θ)​]​​

这里应用了Jensen不等式得到了上面出现过的$ELBO$ELBO，这里的$f(x)$f(x)函数也就是$log$log函数，显然这是一个凹函数。当$log\frac{P(x,z|\theta )}{q(z)}$logq(z)P(x,z∣θ)​这个函数是一个常数时会取得等号：

$\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 )}}$q(z)p(x,z∣θ)​=C⇒q(z)=Cp(x,z∣θ)​⇒∫z​q(z)dz=∫z​C1​p(x,z∣θ)dz⇒1=C1​∫z​p(x,z∣θ)dz⇒C=p(x∣θ)将C代入q(z)=Cp(x,z∣θ)​得q(z)=p(x∣θ)p(x,z∣θ)​=p(z∣x,θ)

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

$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种方法一样了。$logp(x∣θ)≥ELBOEq(z)​[logq(z)p(x,z∣θ)​]​​当q(z)=p(z∣x∣θ)时取等号，因此在迭代更新过程中取q(z)=p(z∣x,θt)。接下来的推导过程就和第1种方法一样了。

四、广义EM

上面介绍的EM算法属于狭义的EM算法，它是广义EM的一个特例。在上面介绍的EM算法的E步中我们假定$q(z)=p(z|x,\theta ^{t})$q(z)=p(z∣x,θt)，但是如果这个后验$p(z|x,\theta ^{t})$p(z∣x,θ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)$logp(x∣θ)=∫z​q(z)logq(z)p(x,z∣θ)​dz−∫z​q(z)logq(z)p(z∣x,θ)​dz=ELBO+KL(q∣∣p)

当我们对于固定的$\theta$θ，我们希望$KL(q||p)$KL(q∣∣p)越小越好，这样就能使得$ELBO$ELBO更大：

$固定\theta ，\hat{q}=\underset{q}{argmin}\; KL(q||p)=\underset{q}{argmax}\; ELBO$固定θ，q^​=qargmin​KL(q∣∣p)=qargmax​ELBO

$ELBO$ELBO是关于$q$q和$\theta$θ的函数，写作$L(q,\theta)$L(q,θ)。以下是广义EM算法的基本思路：

$E\; step：q^{t+1}=\underset{q}{argmax}\; L(q,\theta^{t})\\ M\; step：\theta^{t+1}=\underset{q}{argmax}\; L(q^{t+1},\theta)$Estep：qt+1=qargmax​L(q,θt)Mstep：θt+1=qargmax​L(qt+1,θ)

再次观察一下$ELBO$ELBO：

$ELBO=L(q,\theta )=E_{q}[log\; p(x,z)-log\; q]\\ =E_{q}[log\; p(x,z)]\underset{H[q]}{\underbrace{-E_{q}[log\; q]}}$ELBO=L(q,θ)=Eq​[logp(x,z)−logq]=Eq​[logp(x,z)]H[q]−Eq​[logq]​​

因此，我们看到，⼴义 EM 相当于在原来的式⼦中加⼊熵这⼀项。

五、EM的变种

EM 算法类似于坐标上升法，固定部分坐标，优化其他坐标，再⼀遍⼀遍的迭代。如果在 EM 框架中，⽆法求解$z$z后验概率，那么需要采⽤⼀些变种的 EM 来估算这个后验：
①基于平均场的变分推断，VBEM/VEM
②基于蒙特卡洛的EM，MCEM