EM 算法

  本内容主要介绍 EM 算法 以及其公式的详细推导过程。

  EM 算法 全称 Expectation Maximization Algorithm，译作 期望最大算法 或最大期望化算法，它是一种迭代算法，用于含有 隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。1977 年由 Dempster 等人总结提出。

1.1 Jensen 不等式

  在开始介绍 EM 算法前，我们先介绍和了解一下凹凸函数和 Jensen 不等式。如果您已经了解它们，可以直接跳过。

1.1.1 凸函数和凹函数

  凸函数 是一个定义在某个向量空间的子集 $C$（区间）上的实值函数 $f$，如果在其定义域 $C$ 上的任意两点 $x_1$，$x_2$，$0<t<1$，有

$$tf(x_1)+(1-t)f(x_2)\ge f(t{x}_1+(1-t)x_2)\text{\hspace{1em}(1.1)}$$

也就是说凸函数任意两点的割线位于函数图形上面，这也是 Jensen 不等式的两点形式。如果总有 $tf(x_1)+(1-t)f(x_2)>f(t{x}_1+(1-t)x_2)$，则称函数 $f$ 为严格凸的。

  凹函数 满足
$$tf(x_1)+(1-t)f(x_2)\le f(t{x}_1+(1-t)x_2)\text{\hspace{1em}(1.2)}$$
同理，如果总有 $tf(x_1)+(1-t)f(x_2)<f(t{x}_1+(1-t)x_2)$，则称函数 $f$ 为严格凹的。

  从导数的角度理解，若在其定义域 $C$ 上，函数 $f$ 二次可微，如果 $f^{{}^{\prime \prime }}(x)\ge 0$，那么 $f$ 为凸函数；反之，如果 $f^{{}^{\prime \prime }}(x)\le 0$，那么 $f$ 为凹函数。

图 1.1 凸函数

1.1.2 Jensen 不等式

1.1.2.1 Jensen 不等式

  对于任意点集 { x i } \{x_i\} {xi​}，若 ${\lambda }_i\ge 0$ 且 ${\sum }_i{\lambda }_i=1$，函数 $f(x)$ 满足：

$$\begin{array}{r}f(\sum _{i=1}^m{\lambda }_i{x}_i)\le \sum _{i=1}^m{\lambda }_{i}f(x_i),如果f(x)为凸函数\\ f(\sum _{i=1}^m{\lambda }_i{x}_i)\ge \sum _{i=1}^m{\lambda }_{i}f(x_i),如果f(x)为凹函数\end{array}\text{\hspace{1em}(1.3)}$$

上式被称为 Jensen 不等式。

  在概率论中，如果把 ${\lambda }_i$ 看成取值为 $x_i$ 的离散变量 $X$ 的概率分布，那么 Jensen 不等式可以写成
$$\begin{array}{r}f(E[X])\le E[f(X)],如果f(x)为凸函数\\ f(E[X])\ge E[f(X)],如果f(x)为凹函数\end{array}\text{\hspace{1em}(1.4)}$$
其中，$E[\cdot ]$ 表示期望。此外，如果 $f(x)$ 为严格凸函数或者凹函数，当且仅当 $X=E[X]$时，$E[f(X)]=f(E[X])$ 成立（即当 $X$ 为一个常数时）。

  对于连续变量，Jensen 不等式给出了积分的函数值和函数的积分值间的关系：

$$\begin{array}{r}f(\int xp(x)dx)\le \int f(x)p(x)dx,如果f(x)为凸函数\\ f(\int xp(x)dx)\ge \int f(x)p(x)dx,如果f(x)为凹函数\end{array}\text{\hspace{1em}(1.5)}$$

1.1.2.2 Jensen 不等式的证明过程

  可以使用数学归纳法证明 Jensen 不等式成立，下面我们来证明当 $f(x)$ 为凸函数的情况，为凹函数时可以采用同样的方法进行证明。

  当 $i=1,2$ 时，由凸函数的定义式（1.1），可知其成立。

  假设当 $i=m$ 时，其成立；当 $i=m+1$ 时，得

$$f(\sum _{i=1}^{m+1}{\lambda }_i{x}_i)=f({\lambda }_{m+1}{x}_{m+1}+\sum _{i=1}^m{\lambda }_i{x}_i)\text{\hspace{1em}(1.6)}$$

我们令 ${\eta }_i=\frac{{\lambda }_i}{1-{\lambda }_{m+1}}$，得

$$f(\sum _{i=1}^{m+1}{\lambda }_i{x}_i)=f({\lambda }_{m+1}{x}_{m+1}+(1-{\lambda }_{m+1})\sum _{i=1}^m{\eta }_i{x}_i)\text{\hspace{1em}(1.7)}$$

由凸函数的定义得

$$f(\sum _{i=1}^{m+1}{\lambda }_i{x}_i)\le {\lambda }_{m+1}f(x_{m+1})+(1-{\lambda }_{m+1})f(\sum _{i=1}^m{\eta }_i{x}_i)\text{\hspace{1em}(1.8)}$$

  因为 ${\sum }_{i=1}^{m+1}{\lambda }_i=1$，可得 ${\sum }_{i=1}^m{\lambda }_i=1-{\lambda }_{m+1}$，所以

$$\sum _{i=1}^m{\eta }_i=\frac{{\sum }_{i=1}^m{\lambda }_i}{1-{\lambda }_{m+1}}=1\text{\hspace{1em}(1.9)}$$

由式（1.3）和（1.9）得

$$f(\sum _{i=1}^m{\eta }_i{x}_i)\le \sum _{i=1}^m{\eta }_{i}f(x_i)\text{\hspace{1em}(1.10)}$$

由式（1.8）和（1.10）得

$$\begin{array}{rl}f(\sum _{i=1}^{m+1}{\lambda }_i{x}_i)& \le {\lambda }_{m+1}f(x_{m+1})+(1-{\lambda }_{m+1})\sum _{i=1}^m{\eta }_{i}f(x_i)\\ & ={\lambda }_{m+1}f(x_{m+1})+\sum _{i=1}^m{\lambda }_{i}f(x_i)\\ & =\sum _{i=1}^{m+1}{\lambda }_{i}f(x_i)\end{array}\text{\hspace{1em}(1.11)}$$

因此当 $i=m+1$ 时，Jensen 不等式也成立。

  综上，Jensen 不等式成立。

1.2 EM 算法

  阅读下面的内容，需要您对似然函数和极大似然估计有一定了解，如果您还不了解或者想温习一下，可以参考 这里。

  概率模型有时既含有观测变量（observable variable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。EM 算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

  下面我们从 Chuong B Do & Serafim Batzoglou 在论文《What is the expectation maximization algorithm?》中的抛硬币的例子开始。论文中的截图如下：

图 1.2 抛硬币实验

  我们将论文中的实验简单整理后如下：

---

场景：假设有两枚硬币（硬币 $A$ 和硬币 $B$），硬币 $A$ 正面朝上的概率为 ${\theta }_A$，硬币 $B$ 正面朝上的概率为 ${\theta }_B$。

实验：从两枚硬币中随机选择 1 枚硬币（选中的概率是相等的），使用选中的硬币抛十次。重复 5 次这样的操作，最终将得到 50 次的抛硬币结果。

目的：根据实验结果，推断出 ${\theta }_A$ 和 ${\theta }_B$ 的值。

---

  如果我们知道每次抛的硬币是硬币 $A$ 或者硬币 $B$，那我们可以直接使用极大似然估计求解 ${\theta }_A$ 和 ${\theta }_B$ 的值。但是当前我们不知道每次抛的到底是哪个硬币，即存在隐变量，这个时候我们需要使用 EM 算法求解 ${\theta }_A$ 和 ${\theta }_B$ 的值。

EM 算法流程：

（1）参数初始化：假设 ${\hat{\theta }}_A^{(0)}=0.60$，${\hat{\theta }}_B^{(0)}=0.50$。

（2）E 步：计算后验概率，即每一次抛的硬币分别为 $A$ 和 $B$ 的概率。

  比如第一次抛硬币的结果为 HTTTHHTHTH，即 5 正 5 反，该硬币分别为 $A$ 和 $B$ 的概率为
$$P(A)=\frac{{\theta }_A^5(1-{\theta }_A{)}^5}{{\theta }_A^5(1-{\theta }_A{)}^5+{\theta }_B^5(1-{\theta }_B{)}^5}=\frac{0.6^5\ast 0.4^5}{0.6^5\ast 0.4^5+0.5^5\ast 0.5^5}\approx 0.45$$

$$P(B)=1-P(A)\approx 0.55$$

 使用上面的计算方法，可以求出 5 次实验中，硬币分别为 $A$ 和 $B$ 的概率如表 1.1 所示。

表 1.1 硬币分别为 A 和 B 的概率

次数	为硬币 $A$ 的概率	为硬币 $B$ 的概率
1	0.45	0.55
2	0.80	0.20
3	0.73	0.27
4	0.35	0.65
5	0.65	0.35

  然后根据上面得出的概率，计算硬币 $A$ 和硬币 $B$ 出现正反面的期望。比如第一次抛硬币为 5 正 5 反，则硬币 $A$ 分别为正面和反面的期望为

$$E(A_H)=P(A)\ast 5=0.45\ast 5\approx 2.2$$

$$E(A_T)=P(A)\ast 5=0.45\ast 5\approx 2.2$$

  使用上面的计算方法，可以求出 5 次实验中，硬币 $A$ 和硬币 $B$ 出现正反面的期望如表 1.2 所示。

表 1.2 硬币 A 和硬币 B 出现正反面的期望

次数	硬币 $A$	硬币 $B$
1	$\approx$ 2.2 H, 2.2 T	$\approx$ 2.8 H, 2.8 T
2	$\approx$ 7.2 H, 0.8 T	$\approx$ 1.8 H, 0.2 T
3	$\approx$ 5.9 H, 1.5 T	$\approx$ 2.1 H, 0.5 T
4	$\approx$ 1.4 H, 2.1 T	$\approx$ 2.6 H, 3.9 T
5	$\approx$ 4.5 H, 1.9 T	$\approx$ 2.5 H, 1.1 T
合计	$\approx$ 21.3 H, 8.6 T	$\approx$ 11.7 H, 8.4 T

（3）M 步：计算新的参数 ${\hat{\theta }}_A$ 和 ${\hat{\theta }}_B$

$${\hat{\theta }}_A^{(1)}\approx \frac{21.3}{21.3+8.6}\approx 0.71$$

$${\hat{\theta }}_B^{(1)}\approx \frac{11.7}{11.7+8.4}\approx 0.58$$

（4）进行迭代：重复 E 步和 M 步，直到收敛。

  经过十次迭代后，得到 ${\hat{\theta }}_A^{(10)}\approx 0.80$，${\hat{\theta }}_B^{(10)}\approx 0.52$。

1.3 EM 算法的推导

  通过上面抛硬币的例子，我们已经大致了解了 EM 算法，下面我们开始详细介绍 EM 算法的推导过程。

  给定 $m$ 个训练样本 { x ( 1 ) , ⋯   , x ( m ) } \{x^{(1)}, \cdots, x^{(m)}\} {x(1),⋯,x(m)}，假设样本间相互独立，我们希望将模型 $p(x,z)$ 的参数与数据进行拟合，其似然函数为：

$$\begin{array}{rl}l(\theta )& =\sum _{i=1}^m\log p(x;\theta )\\ & =\sum _{i=1}^m\log \sum _{z}p(x,z;\theta )\end{array}\text{\hspace{1em}(1.12)}$$

  但是，直接求解参数 $\theta$ 的极大似然估计一般会比较困难，因为上式存在一个隐变量 $z$。通常情况下，如果确定 $z$ 后，求解 $\theta$ 就很容易了。

  针对存在含有隐变量的情况下，EM 算法提供了一种有效的极大似然估计方法。因为无法直接最大化 $l(\theta )$，所以采用此方法：不断地建立 $l(\theta )$ 的下界（E步），然后优化下界（M步）。这句话比较抽象，我们继续往下看。

  对每一个样例 $i$，让 $Q_i$ 表示表示该样例隐变量 $z$ 的某种分布（存在 ${\sum }_z{Q}_i(z)=1$，$Q_i(z)\ge 0$）需要注意：如果 $Q_i$ 是连续性的，则 $Q_i$ 表示概率密度函数，需要将求和符号换成积分符号。

  对式（1.12）进行变换得：

$$\begin{array}{rl}\sum _i\log p(x^{(i)};\theta )& =\sum _i\log \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )\\ & =\sum _i\log \sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\\ & \ge \sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\end{array}\text{\hspace{1em}(1.13)}$$

对分子和分母同时乘以 $Q_i(z^{(i)})$，所以式（1.13）中的第二个等号成立。根据 Jensen 不等式，式（1.13）中的不等式成立。在这里，$f(x)=\log (x)$，由于 $\log (x)$ 的二阶导数为 $-\frac{1}{x^2}<0$，所以其为凹函数。我们可以把 $Q_i(z^{(i)})$ 看做概率分布 $p$，那么 ${\sum }_{z^{(i)}}{Q}_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$ 可以看做是 $\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$ 的期望。根据 Jensen 不等式可得

$$f\left(E_{z^{(i)}\sim Q_i}\left[\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\right]\right)\ge E_{z^{(i)}\sim Q_i}\left[f\left(\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\right)\right]\text{\hspace{1em}(1.14)}$$

这样就得到式（1.13）中的不等式成立。

  我们可以把上式写成：$L(\theta )\ge J(z,Q)$（$z$ 为隐变量），那么我们可以通过不断的最大化 $J(z,Q)$，来使得 $L(\theta )$ 不断提高，最终达到它的最大值。图 1.2 更形象地描述这个过程：

图 1.3 EM 算法

这里来说说上图的内在含义：首先我们固定 $\theta$，调整 $Q(z)$ 使下界 $J(z,Q)$ 与 $L(\theta )$ 在此点 $\theta$ 处相等（即绿色曲线到蓝色曲线），然后固定 $Q(z)$，调整 $\theta$ 使下界 $J(z,Q)$ 达到最大值（即 ${\theta }^{(t)}$ 到 ${\theta }^{(t+1)}$）；然后再固定 $\theta$，调整 $Q(z)$，……，直到收敛到 $L(\theta )$ 的最大值处 ${\theta }^{\ast }$。

  在上面的迭代过程中，存在以下两个问题：

1. 什么时候下界 $J(z,Q)$ 与 $L(\theta )$ 在此点 $\theta$ 处相等？
2. 为什么一定会收敛？
   
   

1.3.1 什么时候下界 $J(z,Q)$ 与 $L(\theta )$ 在此点 $\theta$ 处相等？

  在前面介绍 Jensen 不等式时提到，当自变量 $X=E(x)$ 时，即为常数的时候，等式成立。换言之，为了使式（1.13）中的不等式取等号，需要满足

$$\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}=c\text{\hspace{1em}(1.15)}$$

其中，$c$ 为常数，不依赖于 $z^{(i)}$。对上面的等式做一下变换，得

$$p(x^{(i)},z^{(i)};\theta )=c{Q}_i(z^{(i)})\text{\hspace{1em}(1.16)}$$

对上面的等式两边对 $z$ 求和，得

$$\begin{array}{rl}\sum _{z}p(x^{(i)},z^{(i)};\theta )& =\sum _{z}c{Q}_i(z^{(i)})\\ & =c\sum Q_i(z^{(i)})\end{array}\text{\hspace{1em}(1.17)}$$

  因为 $\sum Q_i(z^{(i)})=1$（概率之和为 1），得

$$\sum _{z}p(x^{(i)},z^{(i)};\theta )=c\text{\hspace{1em}(1.18)}$$

  由式（1.15）和式（1.18）得

$$\begin{array}{rl}{Q}_i(z^{(i)})& =\frac{p(x^{(i)},z^{(i)};\theta )}{c}\\ & =\frac{p(x^{(i)},z^{(i)};\theta )}{{\sum }_{z}p(x^{(i)},z^{(i)};\theta )}\\ & =\frac{p(x^{(i)},z^{(i)};\theta )}{p(x^{(i)};\theta )}\\ & =p(z^{(i)}|x^{(i)};\theta )\end{array}\text{\hspace{1em}(1.19)}$$

至此，我们推出了在固定参数 $\theta$ 后，$J(z,Q)$ 与 $L(\theta )$ 相等时，$Q_i(z^{(i)})$ 的取值就是其后验概率（在给定 $x^{(i)}$ 和 $\theta$ 后），这样我们同时解决了 $Q_i(z^{(i)})$ 的选值问题。此步就是 E 步，即建立 $L(\theta )$ 的下界。接下来是 M 步，即在给定 $Q_i(z^{(i)})$ 后，调整 $\theta$，从而极大化 $L(\theta )$ 的下界 $J(z,Q)$。不断地重复 E 步和 M 步，直至收敛，这就是 EM 算法。

  EM 算法的完整步骤如下：

1. 参数初始化：随机初始化参数 ${\theta }^{(0)}$。
2. E 步：根据当前参数 ${\theta }^{(t)}$（初始值 ${\theta }^{(0)}$ 或上一次迭代中 M 步求得的 $\theta$ 值）求隐变量的后验概率 $Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};{\theta }^{(t)})$，即式（1.19）。
3. M 步：固定 $Q_i(z^{(i)})$，通过极大化 $J(z,Q)$ 求得新的参数 ${\theta }^{(t)}$。
4. 进行迭代：重复 E 步和 M 步，直到收敛。
   
   

1.3.2 EM 算法的收敛性

  我们怎么确保 EM 算法一定会收敛呢？首先，假设 ${\theta }^{(t)}$ 和 ${\theta }^{(t+1)}$ 是 EM 算法第 $t$ 次和 $t+1$ 次迭代后的结果。如果我们证明了 $l({\theta }^{(t)})\le l({\theta }^{(t+1)})$，也就是说对数似然函数单调递增，那么最终就会得到最大值。

证明过程：

  在选定 ${\theta }^{(t)}$ 后，由 E 步得
$$Q_i^{(t)}(z^{(i)}):=p(z^{(i)}|x^{(i)};{\theta }^{(t)})\text{\hspace{1em}(1.20)}$$
  当 $Q_i$ 为后验概率时保证了 Jensen 不等式中的等号成立，即得

$$l({\theta }^{(t)})=\sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{Q_i^{(t)}(z^{(i)})}\text{\hspace{1em}(1.21)}$$

  参数 ${\theta }^{(t+1)}$ 是通过极大化式（1.21）得到。经过一些推导会有一下式子成立

$$\begin{array}{rl}l({\theta }^{(t+1)})& \ge \sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t+1)})}{Q_i^{(t)}(z^{(i)})}\\ & \ge \sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{Q_i^{(t)}(z^{(i)})}\\ & =l({\theta }^{(t)})\end{array}\text{\hspace{1em}(1.22)}$$

  下面来具体解释一下式（1.22）。对下面的式子（即式（1.13））

$$l(\theta )\ge \sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\text{\hspace{1em}(1.23)}$$

中的参数分别取值 $Q_i=Q_i^{(t)}$ 和 $\theta ={\theta }^{(t+1)}$ 得

$$l({\theta }^{(t+1)})\ge \sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t+1)})}{Q_i^{(t)}(z^{(i)})}\text{\hspace{1em}(1.24)}$$

从而式（1.22）中的第一个不等号成立。

  因为参数 ${\theta }^{(t+1)}$ 是通过极大化式（1.21）得到，所以可得

$$\sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t+1)})}{Q_i^{(t)}(z^{(i)})}\ge \sum _i\sum _{z^{(i)}}{Q}_i^{(t)}(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{Q_i^{(t)}(z^{(i)})}\text{\hspace{1em}(1.25)}$$

从而式（1.22）中的第二个不等号成立。

  由式（1.21）可得式（1.22）中的等号成立。

  至此，我们证明了 $l({\theta }^{(t)})\le l({\theta }^{(t+1)})$。这样也就证明了 EM 算使得似然函数 $l(\theta )$ 单调递增，并收敛。在实际应用中，可以使用如下方法判断是否收敛：$l(\theta )$ 不再变化或者变化很小。

  从上面的推导可以看出，EM 算法可以保证收敛到一个稳定点，但是却不能保证收敛到全局的极大值点，因此它是局部最优的算法。所以在应用中，初值的选择变得非常重要，常用的办法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

  从前面的推导中我们知道 $l(\theta )\ge J(Q,\theta )$，EM 算法可以看做是 $J(Q,\theta )$ 的坐标上升法。E 步固定 $\theta$，优化 $Q$；M 步固定 $Q$，优化 $\theta$。

1.4 EM 算法的应用

  K-means 算法是 EM 算法思想的体现，E 步骤为聚类过程，M 步骤为更新类簇中心。高斯混合模型（GMM）也是 EM 算法的一个应用。

  缺点：对初始值敏感，EM 算法需要初始化参数 $\theta$，而参数 $\theta$ 的选择直接影响收敛效率以及能否得到全局最优解。

参考

[1] 李航《统计学习方法》
[2] 周志华《机器学习》
[3] http://cs229.stanford.edu/notes/cs229-notes8.pdf
[4] Chuong B Do & Serafim Batzoglou《What is the expectation maximization algorithm?》
[5] 维基百科：凸函数
[6] 维基百科：凹函数
[7] Jensen不等式初步理解及证明
[8] 如何通俗理解EM算法
[9] 从最大似然到EM算法浅解
[10] EM-最大期望算法
[11] （EM算法）The EM Algorithm