EM算法的推导、证明和例子

EM

推导 ¹

X = { x 1 , x 2 , ⋯   , x N } \mathcal{X}=\{x_1, x_2, \cdots, x_N\} X={x1​,x2​,⋯,xN​}: 观察数据

z: 潜在变量

似然估计函数
$$\log P(\mathcal{X};\theta )=\log \prod _i^{N}P(x_i;\theta )=\sum _i^N\log P(x_i;\theta )\text{\hspace{1em}(1)}$$
我们的目标是找到最大化似然估计函数的参数$\theta$
$$\underset{\theta }{\mathrm{argmax}}\log P(\mathcal{X};\theta )\text{\hspace{1em}(2)}$$

现有
$$P(x,z;\theta )=P(x;\theta )P(z|x;\theta )\text{\hspace{1em}(3)}$$
因此
$$\log P(x;\theta )=\log \frac{P(x,z;\theta )}{P(z|x;\theta )}=\log P(x,z;\theta )-\log P(z|x;\theta )\text{\hspace{1em}(4)}$$
假设z属于$Q(z;\varphi )$分布
$$\begin{array}{rl}\log P(x;\theta )& =\log P(x,z;\theta )-\log P(z|x;\theta )-\log Q(z;\varphi )-\log Q(z;\varphi )\\ & =(\log P(x,z;\theta )-\log Q(z;\varphi ))-(\log P(z|x;\theta )-\log Q(z;\varphi ))\\ & =\log \frac{P(x,z;\theta )}{Q(z;\varphi )}-\log \frac{\log P(z|x;\theta )}{Q(z;\varphi )}\end{array}\text{\hspace{1em}(5)}$$
公式5左右两边求关于z的期望
$$\text{left}={\int }_{z}Q(z;\varphi )\log P(x;\theta )dz=\log P(x;\theta ){\int }_{z}Q(z;\varphi )dz=\log P(x;\theta )\text{\hspace{1em}(6)}$$
因为${\int }_{z}Q(z;\varphi )dz=1$。

$$\begin{array}{rl}\text{right}& ={\int }_{z}Q(z;\varphi )\log \frac{P(x,z;\theta )}{Q(z;\varphi )}dz-{\int }_{z}Q(z;\varphi )\log \frac{\log P(z|x;\theta )}{Q(z;\varphi )}dz\\ & =\text{ELBO}+\text{KL}(Q(z;\varphi )||P(z|x;\theta ))\end{array}\text{\hspace{1em}(7)}$$
公式7第一项称为ELBO(evidence lower bound)，第二项是KL散度。
因此，我们得到
$$\log P(x;\theta )=\text{ELBO}+\text{KL}(Q(z;\varphi )||P(z|x;\theta ))\text{\hspace{1em}(8)}$$
因为$\text{KL}(\cdot )\ge 0$，因此$\log P(x;\theta )\ge \text{ELBO}$，当且仅当$Q(z;\varphi )=P(z|x;\theta )$时，等号成立。ELBO相当于一个下界，不断地提高ELBO，就能不断提高$\log P(x;\theta )$，达到我们的目的——最大化似然估计函数。

假设我们有${\theta }^{(t)}$，我们想要最大化ELBO，即最小化$\text{KL}(Q(z;\varphi )||P(z|x;\theta ))$:
$${\varphi }^{(t)}=\underset{\varphi }{\mathrm{argmin}}\text{KL}(Q(z;\varphi )||P(z|x;{\theta }^{(t)}))=\underset{\varphi }{\mathrm{argmax}}\text{ELBO}(\varphi ,{\theta }^{(t)})\text{\hspace{1em}(9)}$$
当得到最优的${\varphi }^{(t)}$，有$Q(z;{\varphi }^{(t)})=P(z|x;{\theta }^t)$。实际情况下很难得到最优的${\varphi }^{(t)}$，我们的目的时尽可能最大化ELBO。当我们计算出ELBO后，我们可以反过来求${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmax}}\text{ELBO}({\varphi }^{(t)},\theta )$。

通过不断重复这个两个最大化的过程，我们就可以近似地求出最大化似然估计函数的参数$\theta$。

继续看看$\text{ELBO}({\varphi }^{(t)},\theta )$

$$\begin{array}{rl}\text{ELBO}({\varphi }^{(t)},\theta )& ={\int }_{z}Q(z;{\varphi }^{(t)})\log \frac{P(x,z;\theta )}{Q(z;{\varphi }^{(t)})}dz\\ & ={\int }_{z}Q(z;{\varphi }^{(t)})\log P(x,z;\theta )dz-{\int }_{z}Q(z;{\varphi }^{(t)})\log Q(z;{\varphi }^{(t)}dz\\ & =E_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]-E_{z\backsim Q(z;{\varphi }^{(t)})}[\log Q(z;{\varphi }^{(t)}]\end{array}\text{\hspace{1em}(10)}$$
公式10的最后一行的第一项是关于z的期望，第二项是一个常数（${\varphi }^{(t)}$已知），所以
$$\begin{array}{rl}{\theta }^{(t+1)}& =\underset{\theta }{\mathrm{argmax}}\text{ELBO}({\varphi }^{(t)},\theta )\\ & =\underset{\theta }{\mathrm{argmax}}{E}_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]\end{array}\text{\hspace{1em}(11)}$$
因此EM算法叫做期望最大化算法。

总结，EM算法的迭代过程如下

• E-step: 固定${\theta }^{(t)}$，${\varphi }^{(t)}=\underset{\varphi }{\mathrm{argmax}}{E}_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]$；
• M-step: 固定${\varphi }^{(t)}$，${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmax}}{E}_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]$。

E-step和M-step的顺序可以互换。

EM算法收敛性证明¹

简单的证明：只要${\theta }^{(t)}\to {\theta }^{(t+1)},\log P(x;{\theta }^{(t)})\le \log P(x;{\theta }^{(t+1)})$，就能保证算法收敛。

从公式4出发，两边求关于z的期望
$$\begin{array}{rl}\text{left}& ={\int }_{z}Q(z;{\varphi }^{(t)})\log P(x;\theta )dz\\ & =\log P(x;\theta ){\int }_{z}Q(z;{\varphi }^{(t)})dz\\ & =\log P(x;\theta )\end{array}\text{\hspace{1em}(12)}$$

$$\text{right}={\int }_{z}Q(z;{\varphi }^{(t)})\log P(x,z;\theta )dz-{\int }_{z}Q(z;{\varphi }^{(t)})\log P(z|x;\theta )dz\text{\hspace{1em}(13)}$$

因为${\varphi }^{(t)}$根据公式9求解得到的，假设我们得到的是最优解，则$Q(z;{\varphi }^{(t)})=P(z|x;{\theta }^t)$，代入公式13得
$$\begin{array}{rl}\text{right}& ={\int }_{z}P(z|x;{\theta }^{(t)})\log P(x,z;\theta )dz-{\int }_{z}P(z|x;{\theta }^{(t)})\log P(z|x;\theta )dz\\ & =H_1(\theta ,{\theta }^{(t)})-H_2(\theta ,{\theta }^{(t)})\end{array}\text{\hspace{1em}(14)}$$
我们分别用$H_1$和$H_2$指代公式14的两项。

由公式12和14得
$$\log P(x;\theta )=H_1(\theta ,{\theta }^{(t)})-H_2(\theta ,{\theta }^{(t)})\text{\hspace{1em}(15)}$$

因为${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmax}}{E}_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]$，所以有$H_1({\theta }^{(t+1)},{\theta }^{(t)})\ge H_1({\theta }^{(t)},{\theta }^{(t)})$。接下来，只要证明$-H_2({\theta }^{(t+1)},{\theta }^{(t)})\ge -H_2({\theta }^{(t)},{\theta }^{(t)})$，就能证明$\log P(x;{\theta }^{(t+1)})\ge \log P(x;{\theta }^{(t+1)})$。

现在
$$\begin{array}{rl}& H_2({\theta }^{(t+1)},{\theta }^{(t)})-H_2({\theta }^{(t)},{\theta }^{(t)})\\ =& {\int }_{z}P(z|x;{\theta }^{(t)})\log P(z|x;{\theta }^{(t+1)})dz-{\int }_{z}P(z|x;{\theta }^{(t)})\log P(z|x;{\theta }^{(t)})dz\\ =& {\int }_{z}P(z|x;{\theta }^{(t)})\log \frac{P(z|x;{\theta }^{(t+1)})}{P(z|x;{\theta }^{(t)})}dz\end{array}\text{\hspace{1em}(16)}$$
证明公式16小于等于0：

方法1：公式16是负KL散度$-KL(P(z|x;{\theta }^{(t)})||P(z|x;{\theta }^{(t+1)}))\le 0$。

方法2：公式16等于$E_{z\backsim P(z|x,{\theta }^{(t)})}[\log \frac{P(z|x;{\theta }^{(t+1)})}{P(z|x;{\theta }^{(t)})}]$。根据Jensen不等式$E[\log (x)]\le \log E[x]$，因此
$$\begin{array}{rl}& E_{z\backsim P(z|x,{\theta }^{(t)})}[\log \frac{P(z|x;{\theta }^{(t+1)})}{P(z|x;{\theta }^{(t)})}]\\ \le & \log E_{z\backsim P(z|x,{\theta }^{(t)})}[\frac{P(z|x;{\theta }^{(t+1)})}{P(z|x;{\theta }^{(t)})}]\\ =& \log {\int }_{z}P(z|x;{\theta }^{(t)})\frac{P(z|x;{\theta }^{(t+1)})}{P(z|x;{\theta }^{(t)})}dz\\ =& \log {\int }_{z}P(z|x;{\theta }^{(t+1)})dz\\ =& \log 1=0\end{array}\text{\hspace{1em}(17)}$$

例子 ²

抛硬币，有两个硬币，但是两个硬币的材质不同导致其出现正反面的概率不一样，目前我们只有一组观测数据，要求出每一种硬币投掷时正面向上的概率。总共投了五轮，每轮投掷五次。假设我们不知道每一次投掷用的是哪一种硬币，等于是现在的问题加上了一个隐变量，就是每一次选取的硬币的种类。

(图片来自https://blog.csdn.net/u010834867/article/details/90762296)

设两个硬币分别是AB，P(正|A)=$x_1$，P(反|A)=$1-x_1$，P(正|B)=$x_2$，P(正|B)=$x_2$。

假设第i次实验选择硬币A的概率是$P(z_i=A)=y_i$，选择硬币B的概率是$P(z_i=B)=1-y_i$。

看实验i的数据j，用$x_{ij}$表示，似然估计函数为
$$\log P(x)=\log \prod _i\prod _{j}P(x_{ij})=\sum _i\sum _j\log P(x_{ij})$$
$$P(x_{ij})=\frac{P(x,z)}{P(z|x)}=\frac{P(z)P(x|z)}{P(z|x)}$$

其中$P(z|x)$不好求出来。我们使用EM算法来解$x_1$和$x_2$。

首先，我们求出期望，
$$\begin{array}{rlrl}& E_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]& & \\ =& y_1\log y_1(x_1{x}_1(1-x_1)x_1(1-x_1))& +& (1-y_1)\log (1-y_1)(x_2{x}_2(1-x_2)x_2(1-x_2))\\ +& y_2\log y_2((1-x_1)(1-x_1)x_1{x}_1(1-x_1))& +& (1-y_2)\log (1-y_2)((1-x_2)(1-x_2)x_2{x}_2(1-x_2))\\ +& y_3\log y_3(x_1(1-x_1)(1-x_1)(1-x_1)(1-x_1))& +& (1-y_3)\log (1-y_3)(x_2(1-x_2)(1-x_2)(1-x_2)(1-x_2))\\ +& y_4\log y_4(x_1(1-x_1)(1-x_1)x_1{x}_1)& +& (1-y_4)\log (1-y_4)(x_2(1-x_2)(1-x_2)x_2{x}_2)\\ +& y_5\log y_5((1-x_1)x_1{x}_1(1-x_1)(1-x_1))& +& (1-y_5)\log (1-y_5)((1-x_2)x_2{x}_2(1-x_2)(1-x_2))\end{array}$$

假设$x_1=0.2,x_2=0.7$，代入上式得
$$\begin{array}{rlrl}& E_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]& & \\ =& y_1\log 0.00512y_1& +& (1-y_1)\log 0.03087(1-y_1)\\ +& y_2\log 0.02048y_2& +& (1-y_2)\log 0.01323(1-y_2)\\ +& y_3\log 0.08192y_3& +& (1-y_3)\log 0.00567(1-y_3)\\ +& y_4\log 0.00512y_4& +& (1-y_4)\log 0.03087(1-y_4)\\ +& y_5\log 0.02048y_5& +& (1-y_5)\log 0.01323(1-y_5)\end{array}$$

(图片来自https://blog.csdn.net/u010834867/article/details/90762296)

现在求
$$\max E_{z\backsim Q(z;{\varphi }^{(t)})}[\log P(x,z;\theta )]$$

为了简单运行，我们取$Q(z;{\varphi }^{(t)})$为 y = { 0 , 1 , 1 , 0 , 1 } y=\{0,1,1,0,1\} y={0,1,1,0,1}，即 z = { B , A , A , B , A } z=\{B,A,A,B,A\} z={B,A,A,B,A}。虽然求出来的期望不是最大的，但不影响算法的收敛。因为z的结果已经固定了，可以直接计算${\theta }^{(t+1)}$：
$$x_1=(2+1+2)/15=0.33,x_2=(3+3)/10=0.6$$

接着不断迭代，直到z或者$x_1$和$x_2$收敛。

若有不恰当之处，请指正。

---

1. https://www.bilibili.com/video/av31906558?p=1 ↩︎ ↩︎

2. https://blog.csdn.net/u010834867/article/details/90762296 ↩︎