EM算法详解

1.写在前面的话

本文主要涉及以下知识：

• 概率论
  • 联合概率
  • 条件概率
  • 边际概率
  • 期望
  • 贝叶斯公式
• 极大似然估计
• KL散度的定义及性质

对于本文涉及的上述知识，本文不会用篇幅进行说明。有关极大似然估计和KL散度的内容可以参考

• 一元正态总体参数的极大似然估计_附Python实现
• KL散度(Kullback-Leibler Divergence)介绍及详细公式推导

本文仅作为学习记录用，并非专业的教程，内容大多来自参考书目，同时加以自己的理解而作，内容仅供参考。若要详细了解，请参看参考文献书目。文中如有错误，请联系修改。

2.EM算法

2.1.EM算法推导

极大似然估计是在已观测数据上对模型参数进行估计。但在某些场景中，数据存在未观测变量，称之为潜变量或隐变量，用$Z$表示，而已观测数据又称为不完全数据，用$X$表示。(X,Z)称之为完全数据。由于$Z$并未观测到，这时使用极大似然估计将不可行。为了解决有潜变量的参数估计问题，Dempster et al和McLachlan and Krishnan于1977年提出了EM算法。
先考虑一个概率模型，并将已观测数据联合起来使用$X$表示，将所有隐变量记作$Z$，联合概率$p(X,Z|\theta )$由参数向量$\theta$控制，我们的目标是最大化似然函数
$$p(X|\theta )=\sum _{Z}p(X,Z|\theta )\text{\hspace{1em}(1)}$$

上式的对数似然函数为

$$\ln p(X|\theta )=\ln \sum _{Z}p(X,Z|\theta )\text{\hspace{1em}(2)}$$

对其对数似然函数求导进行优化将非常困难。现在我们假设对完全数据的对数似然函数$\ln p(X,Z|\theta )$进行优化是容易的。接下来我们定义一个在隐变量上的分布$q(Z)$。考虑到$\ln p(X|\theta )$与$q(Z)$无关，则有下式成立
$$\ln p(X|\theta )=\sum _{Z}q(Z)\ln p(X|\theta )\text{\hspace{1em}(3)}$$

又考虑到
$$p(X|\theta )=\frac{p(X,Z|\theta )}{p(Z|X,\theta )}\text{\hspace{1em}(4)}$$

将上式代入(3)式右侧得
$$\begin{array}{rl}\ln p(X|\theta )& =\sum _{Z}q(Z)\ln \frac{p(X,Z|\theta )}{p(Z|X,\theta )}\\ & =\sum _{Z}q(Z)\ln \left\{\frac{p(X,Z|\theta )}{q(Z)}\frac{q(Z)}{p(Z|X,\theta )}\right\}\\ & =\sum _{Z}q(Z)\left\{\ln \frac{p(X,Z|\theta )}{q(Z)}+\ln \frac{q(Z)}{p(Z|X,\theta )}\right\}\\ & =\sum _{Z}q(Z)\left\{\ln \frac{p(X,Z|\theta )}{q(Z)}\right\}+\sum _{Z}q(Z)\ln \left\{\frac{q(Z)}{p(Z|X,\theta )}\right\}\\ & =\mathcal{L}(q,\theta )+D_{KL}(q\parallel p)\end{array}\text{\hspace{1em}(5)}$$

其中
$$\begin{array}{r}\mathcal{L}(q,\theta )=\sum _{Z}q(Z)\left\{\ln \frac{p(X,Z|\theta )}{q(Z)}\right\}\end{array}\text{\hspace{1em}(6)}$$

$$D_{KL}(q\parallel p)=\sum _{Z}q(Z)\ln \left\{\frac{q(Z)}{p(Z|X,\theta )}\right\}\text{\hspace{1em}(7)}$$

显然(5)式右边第二项为$q$到$p$的KL散度，由KL散度的数学性质$D_{KL}(q\parallel p)\ge 0$(当$q(Z)=p(Z|X,\theta )$时等号成立)可得
$$\mathcal{L}(q,\theta )\le \ln p(X|\theta )\text{\hspace{1em}(8)}$$

换句话说，$\mathcal{L}(q,\theta )$给出了对数似然$\ln p(X|\theta )$的下界。由于直接极大化对数似然$\ln p(X|\theta )$是困难的，而$\mathcal{L}(q,\theta )$恰好是对数似然$\ln p(X|\theta )$的下界，故可以极大化$\mathcal{L}(q,\theta )$以使对数似然$\mathcal{L}(q,\theta )$极大化。
现在来详细说明EM算法的步骤。
E步: 假设当前迭代参数为${\theta }^{(i)}$。现固定参数${\theta }^{(i)}$，同时式子$\mathcal{L}(q,\theta )$关于$Z$最大化。根据式子(8)，$\mathcal{L}(q,\theta )$的最大值等于$\ln p(X|\theta )$，而要想$\mathcal{L}(q,\theta )$等于$\ln p(X|\theta )$，只需$D_{KL}(q\parallel p)=0$，意味着$q(Z)=p(Z|X,{\theta }^{(i)})$，将该式代入式(6)有
$$\begin{array}{rl}\mathcal{L}(q,\theta )& =\sum _{Z}p(Z|X,{\theta }^{(i)})\left\{\ln \frac{p(X,Z|\theta )}{p(Z|X,{\theta }^{(i)})}\right\}\\ & =\sum _{Z}p(Z|X,{\theta }^{(i)})\ln p(X,Z|\theta )-\sum _{Z}p(Z|X,{\theta }^{(i)})\ln p(Z|X,{\theta }^{(i)})\end{array}\text{\hspace{1em}(9)}$$

上式右侧第二项中，参数${\theta }^{(i)},X$已知，因此第二项是一个常数。故上式可改写为
$$\begin{array}{rl}\mathcal{L}(q,\theta )& =\sum _{Z}p(Z|X,{\theta }^{(i)})\ln p(X,Z|\theta )+C\\ & =Q(\theta ,{\theta }^{(i)})+C\end{array}\text{\hspace{1em}(10)}$$

M步: 固定分布$q(Z)$，下界$\mathcal{L}(q,\theta )$关于$\theta$进行极大化，得到新的参数${\theta }^{(i+1)}$。由于$C$是常数，并不影响极大化$\mathcal{L}(q,\theta )$的结果，故等价于极大化$Q(\theta ,{\theta }^{(i)})$,即
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(11)}$$

重复以上两个步骤，直至到达指定迭代次数或收敛为止。
上述两个步骤都在对对数似然$\ln p(X|\theta )$的下界$\mathcal{L}(q,\theta )$进行极大化，区别在于，E步固定参数${\theta }^{(i)}$关于分布$q(Z)$极大化$\mathcal{L}(q,\theta )$，而M步固定分布$q(Z)$关于参数$\theta$极大化$\mathcal{L}(q,\theta )$。其中，E步极大化的结果产生了$Q(\theta ,{\theta }^{(i)})$。正是因此，将其称之为期望步。

2.2.图解EM算法

现将上述迭代过程使用图示进行说明。
假设当前迭代参数为${\theta }^{(i)}$，由式(5)知，对数似然是下界$\mathcal{L}(q,{\theta }^{(i)})$与KL散度$D_{KL}(q\parallel p)$的和。如下图所示：

E步固定参数${\theta }^{(i)}$关于分布$q(Z)$极大化$\mathcal{L}(q,{\theta }^{(i)})$，也就是令$q(Z)=p(Z|X,{\theta }^{(i)})$，于是$D_{KL}(q\parallel p)=0$。从而得到下图

M步固定分布$q(Z)$关于参数$\theta$极大化$\mathcal{L}(q,{\theta }^{(i)})$，此时的分布$q(Z)=p(Z|X,{\theta }^{(i)})$。在关于参数$\theta$极大化$\mathcal{L}(q,{\theta }^{(i)})$之后，
$$\begin{gathered} {\theta }^{(i+1)}\leftarrow {\theta }^{(i)} \\ \mathcal{L}(q,{\theta }^{(i+1)})\leftarrow \mathcal{L}(q,{\theta }^{(i)}) \end{gathered}$$

同时，$q(Z)$的后验分布更新为
$$p(Z|X,{\theta }^{(i+1)})\leftarrow p(Z|X,{\theta }^{(i)})$$

由于似然函数极大化后，参数及$q(Z)$的后验分布得到更新，此时，更新后的后验分布与更新前的后验分布$p(Z|X,{\theta }^{(i)})$的KL散度不为零，也即$D_{KL}(q\parallel p)>0$。根据式(5)知，对数似然$\ln p(X|\theta )$的增加量大于$\mathcal{L}(q,{\theta }^{(i)})$的增加量。用图示如下：

在实际应用中，不断重复上述两步，直到达到指定迭代次数或者收敛为止。
上述迭代过程类似于坐标上升法

2.3.EM算法总结

现将EM算法的步骤总结如下：

EM算法：
输入：观测数据$X$，隐变量数据$Z$，初始化参数${\theta }^{(0)}$
输出：模型参数$\theta$
(1)$E$步:
记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计，则在第$i+1$次迭代的$E$步中，计算下式
$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}p(Z|X,{\theta }^{(i)})\ln p(X,Z|\theta )$$

(2)$M$步：
求使得$Q(\theta ,{\theta }^{(i)})$极大化的参数${\theta }^{(i+1)}$，即
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

(3) 重复(1)和(2),直至收敛。

3.总结

本文对EM算法的过程进行推导，并给出了EM算法的一般过程。

参考文献

[1] Goodfellow I, Bengio Y, Courville A. Deep learning[M]. MIT press, 2016.
[2] Bishop C M, Nasrabadi N M. Pattern recognition and machine learning[M]. New York: springer, 2006.
[3] 李航. 统计学习方法[M]. Qing hua da xue chu ban she, 2012.