期望最大化算法

1. 历史背景

1955 年，Cedric Smith 提出一种用于估计等位基因频率的基因计数方法。H.O. Hartley 于 1958 年、Hartley and Hocking 于 1977 年分别提出类似的方法，都是解决不完整数据中的最大似然求解问题，它们均为 EM 算法的起源。

1971 年，Rolf Sundberg 在他的几篇论文中探讨了指数族变量不完整数据的最大似然理论，并给出了非常详细的处理。Sundberg 将后来称为 Sundberg 公式的公式归功于其老师 Per Martin-Löf 和 Anders Martin-Löf 以前的手稿。

1977 年，Arthur Dempster、Nan Laird 和 Donald Rubin Dempster-Laird-Rubin 论文概括了 EM 方法，并对更广泛的问题进行了收敛性分析。 Dempster-Laird-Rubin 的论文将 EM 方法确立为统计分析的重要工具。

1983 年，由于 Dempster-Laird-Rubin 算法的收敛性分析存在缺陷，C. F. Jeff Wu 发表了正确的收敛性分析。Wu 的证明确立了 EM 方法在指数族之外的收敛性，正如 Dempster-Laird-Rubin 所声称的那样。

2. 基础预备

推导一些有关的不等式，部分是高中知识。

定义1　如果 $f$ 是定义在区间 $I=[a,b]$ 上的实值函数，且 $\forall x_1,x_2\in I,\lambda\in I$，满足

$$f(\lambda x_1+(1-\lambda)x_2)\le \lambda f(x_1)+(1-\lambda)f(x_2)$$

则称 $f$ 为凸函数，如果不等式是严格的，则 $f$ 是严格凸的。

实际上 $\lambda x_1+(1-\lambda)x_2=\lambda(x_1-x_2)+x_2$，表示的是区间 $[x_1,x_2]$ 中的一个点。不等式的含义就是区间 $[x_1,x_2]$ 的函数值永远不会高于点 $(x_1,f(x_1))$ 到 $(x_2,f(x_2))$ 的直线。$x^2$ 是凸函数，与高数中定义相反。

定义2　如果 $f$ 是凹（严格凹）函数，则 $-f$ 是凸（严格凸）函数。

定理1　如果 $f(x)$ 在 $[a,b]$ 上是两次可微的，并且 $f’’(x)\ge 0$，则 $f(x)$ 在 $[a,b]$ 是凸的。

证明：假设 $x\le y\in[a,b]$，并且 $\lambda\in[0,1]$，让 $z=\lambda y+(1-\lambda)x$，根据定义 1 ，如果满足

$$f(z)\le \lambda f(y)+(1-\lambda)f(x)$$

则可证明 $f$ 是凸函数，$f(z)$ 可以表示为

$$f(z)=\lambda f(z)+(1-\lambda)f(z)$$

联立两个式子可以得到

$$\lambda(f(y)-f(z))\ge(1-\lambda)(f(z)-f(x))$$

根据中值定理，存在 $z\le t\le y$

$$f(y)-f(z)=f'(t)(y-z)$$

同理，存在 $x\le s\le z$

$$f(z)-f(x)=f'(s)(z-x)$$

由于 $f’’(x)\ge 0$，并且 $t\le s$，可知 $f’(t)\le f’(s)$

根据 $z=\lambda y+(1-\lambda)x$，可以得到

$$(1-\lambda)(z-x)=\lambda(y-z)$$

最终，联立以上各式

$$\begin{aligned} \lambda(f(y)-f(z))&=\lambda f'(t)(y-z) \\ &\ge f'(s)\lambda(y-z)\\ &=f'(s)(1-\lambda)(z-x)\\ &=(1-\lambda)(f(z)-f(x)) \end{aligned}$$

推论1　$-\ln(x)$ 是在区间 $(0,\infty)$ 是严格凸函数。

证明：$f(x)=-\ln(x)$，由于 $\displaystyle f’’(x)=\frac{1}{x^2}>0$，因此 $-\ln(x)$ 是在区间 $(0,\infty)$ 是严格凸函数。

将上述凸性的概念扩展到 $n$ 个点，得到的就是 Jensen 不等式。

定理2（Jensen 不等式）　如果 $f$ 是区间 $I$ 上的凸函数，如果 $x_1,x_2,\dots,x_n\in I$ ， $\lambda_1,\lambda_2,\dots,\lambda_n \ge 0$ 且 $\displaystyle\sum_{i=1}^n \lambda_i = 1$，则有

$$f\left(\sum_{i=1}^n \lambda_ix_i\right)\le \sum_{i=1}^n \lambda_if(x_i)$$

证明：$n=2$ 对应于凸性的定义，为了证明这对所有自然数都是正确的，我们通归纳来证明，假设定理对某些 $n$ 成立，则

$$\begin{aligned} f\left(\sum_{i=1}^{n+1} \lambda_ix_i\right) &=f\left(\lambda_{n+1}x_{n+1}+\sum_{i=1}^n \lambda_ix_i\right)\\ &=f\left(\lambda_{n+1}x_{n+1}+(1-\lambda_{n+1})\frac{1}{1-\lambda_{n+1}}\sum_{i=1}^n \lambda_ix_i\right)\\ &\le \lambda_{n+1}f(x_{n+1})+(1-\lambda_{n+1})f\left(\frac{1}{1-\lambda_{n+1}}\sum_{i=1}^n \lambda_ix_i\right) \\ &= \lambda_{n+1}f(x_{n+1})+(1-\lambda_{n+1})f\left(\sum_{i=1}^n \frac{\lambda_i}{1-\lambda_{n+1}}x_i\right) \\ &\le \lambda_{n+1}f(x_{n+1})+(1-\lambda_{n+1})\sum_{i=1}^n \frac{\lambda_i}{1-\lambda_{n+1}}f(x_i) \\ &= \lambda_{n+1}f(x_{n+1})+\sum_{i=1}^n \lambda_if(x_i)\\ &=\sum_{i=1}^{n+1} \lambda_if(x_i) \end{aligned}$$

由于 $\ln$ 是凹函数，因此通过 Jensen 不等式可以推导出

$$\ln\sum_{i=1}^n \lambda_ix_i \ge \sum_{i=1}^n \lambda_i\ln(x_i)$$

这个公式可以得到一个和的对数的下限，可用于推导 EM 算法。

Jensen不等式可以证明算术平均值大于等于几何平均值，即

推论2　

$$\frac{1}{n}\sum_{i=1}^n x_i\ge \sqrt[n]{x_1x_2\dots x_n}$$

证明： 如果 $x_1,x_2,\dots,x_n\ge0$，由于 $\ln$ 是凹函数，可以得到

$$\begin{aligned} \ln\left(\frac{1}{n}\sum_{i=1}^{n}x_i\right) &\ge \sum_{i=1}^{n}\frac{1}{n} \ln(x_i) \\ &= \frac{1}{n}\ln(x_1x_2\dots x_n) \\ &=\ln(x_1x_2\dots x_n)^{\frac{1}{n}} \end{aligned}$$

因此，有

$$\frac{1}{n}\sum_{i=1}^n x_i\ge \sqrt[n]{x_1x_2\dots x_n}$$

3. 算法推导

在概率模型当中，已知样本观察数据，想要得到模型参数，一般使用极大似然估计。

极大似然估计需要对似然函数的所有未知值、参数求导，并同时求解得到的的方程。而在某些统计模型当中，除了已知的观察数据之外，还有未知的潜在变量。话句话说，也就是模型中缺失了部分的观察变量，而缺失的观察变量和已知的之间存在某种概率关系。如果用极大似然算法，它们之间的导数是一组互锁的方程，其中参数的解需要潜在变量的值，反之亦然，但是将一组方程带入另一组方程会产生一个不可解的方程组。

EM 算法就是求解上述问题的一种迭代方法，主要思想是简单地为两组未知数中的一组选择任意值，使用它们来估计第二组，然后使用这些新值来找到对第一组的更好估计，然后在两者之间保持交替，直到结果值都收敛到固定点。

EM 算法是一种通用的算法，无论是离散还是连续，二项式分布或高斯分布，都可以用来求解。

假设 $\mathbf{X}$ 表示观测随机变量的一组结果，在概率模型参数 $\theta$ 下的条件概率为 $\mathcal{P}(\mathbf{X}|\theta)$，找到最大 $\theta$ 使其概率最大的过程就称为最大似然估计。为了估计 $\theta$，通常引入对数似然函数

$$L(\theta)=\ln \mathcal{P}(\mathbf{X}|\theta)$$

如果概率模型中存在一组隐藏变量 $\mathbf{Z}$，它的结果无法通过观测得到，但其结果会对观测变量 $\mathbf{X}$ 产生影响。$\mathbf{X}$ 和 $\mathbf{Z}$ 一起称为完全数据，观测变量 $\mathbf{X}$ 则又称为不完全数据。假定隐藏变量 $\mathbf{Z}$ 的一组取值为 $\mathbf{z}$，则观测结果 $\mathbf{X}$ 是基于 $\mathbf{z}$ 和参数 $\theta$ 的条件之下产生的，

$$\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)=\frac{\mathcal{P}(\mathbf{X}, \mathbf{z}, \theta)}{\mathcal{P}(\mathbf{z}, \theta)}$$

而隐藏变量同样是基于 $\theta$ 产生的，

$$\mathcal{P}(\mathbf{z}|\theta)=\frac{\mathcal{P}(\mathbf{z}, \theta)}{\mathcal{P}(\theta)}$$

将两个条件概率相乘，并将隐藏变量所有可能的取值进行累加（隐藏变量的取值似乎都是可数的），就可以得到 $\mathbf{X}$ 在 $\theta$ 下的条件概率，

$$\sum_{\mathbf{z}}\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta)=\sum_{\mathbf{z}}\frac{\mathcal{P}(\mathbf{X}, \mathbf{z}, \theta)}{\mathcal{P}(\theta)}=\mathcal{P}(\mathbf{X}|\theta)$$

因此将对数似然函数改写为

$$L(\theta) = \ln\mathcal{P}(\mathbf{X}|\theta)=\ln\sum_{\mathbf{z}}\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta)$$

最终是求 $\theta$ 使上式结果最大化，该问题没有解析解，只能通过迭代的方法求解。假设第 $t$ 次迭代后 $\theta$ 的估计值为 $\theta_t$，我们希望新估计值 $\theta$ 能使 $L(\theta)$ 增加，考虑两者差值

$$L(\theta)-L(\theta_t)=\ln\sum_{\mathbf{z}}\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta) - \ln\mathcal{P}(\mathbf{X}|\theta_t)$$

这是一个变量为 $\theta$ 的函数，$\theta_t$ 相当于常数，$\sum$ 在 $\ln$ 里面是无法对求解函数极值的，因此利用上述 Jensen 不等式

$$\ln\sum_{i=1}^n \lambda_ix_i \ge \sum_{i=1}^n \lambda_i\ln(x_i)$$

为了提出 $\lambda_i$，需要凑出 $\displaystyle\sum_{i=1}^n \lambda_i = 1$，可以利用条件概率 $\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)$，我们有 $\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ge 0$ 且 $\displaystyle\sum_\mathbf{z}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)= 1$，因此

$$\begin{aligned} L(\theta)-L(\theta_t)&=\ln\sum_{\mathbf{z}}\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta) - \ln\mathcal{P}(\mathbf{X}|\theta_t)\\ &=\ln\sum_{\mathbf{z}}\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta)\frac{\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)}{\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)} - \ln\mathcal{P}(\mathbf{X}|\theta_t)\\ &\ge\sum_{\mathbf{z}}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ln\frac{\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta)}{\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\mathcal{P}(\mathbf{X}|\theta_t)} \end{aligned}$$

令

$$l(\theta|\theta_t)=L(\theta_t)+\sum_{\mathbf{z}}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ln\frac{\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta)}{\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\mathcal{P}(\mathbf{X}|\theta_t)}$$

则有

$$L(\theta)\ge l(\theta|\theta_t)$$

因此我们有一个函数 $l(\theta|\theta_t)$ 来作为对数似然函数 $L(\theta)$ 的下界，并且当 $\theta=\theta_t$ 时，可以推导出

$$l(\theta_t|\theta_t)=L(\theta_t)$$

因此，任何可以使 $l(\theta|\theta_t)$ 增大的 $\theta$，也可以使 $L(\theta)$ 增大，为了使 $L(\theta)$ 尽可能大的增长，选择 $\theta_{t+1}$ 使 $l(\theta|\theta_t)$ 达到极大值，只与 $\theta_t$ 有关的项为常数项，即

$$\begin{aligned} \theta_{t+1}&=\arg\max_{\theta} l(\theta|\theta_t)\\ &=\arg\max_{\theta} \sum_{\mathbf{z}}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ln\mathcal{P}(\mathbf{X}|\mathbf{z}, \theta)\mathcal{P}(\mathbf{z}|\theta)\\ &=\arg\max_{\theta} \sum_{\mathbf{z}}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ln\frac{\mathcal{P}(\mathbf{X}, \mathbf{z}, \theta)}{\mathcal{P}(\mathbf{z}, \theta)}\frac{\mathcal{P}(\mathbf{z}, \theta)}{\mathcal{P}(\theta)}\\ &=\arg\max_{\theta} \sum_{\mathbf{z}}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ln\mathcal{P}(\mathbf{X}, \mathbf{z}|\theta)\\ &=\arg\max_{\theta} E_{\mathbf{z|\mathbf{X},\theta_t}}\{\ln\mathcal{P}(\mathbf{X}, \mathbf{z}|\theta)\}\\ \end{aligned}$$

最后的一项是完全数据的对数似然函数 $\ln\mathcal{P}(\mathbf{X}, \mathbf{z}|\theta)$ 关于给定观测数据 $\mathbf{X}$ 和当前参数 $\theta_t$ 下对未观测数据 $\mathbf{z}$ 的条件概率分布 $\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)$ 的期望。

总结如下

EM 算法

• 输入：观测变量 $\mathbf{X}$
• 输出：模型参数 $\theta$
• 步骤：
  1. 初始化模型参数为 $\theta_0$
  2. （E 步）假设第 $t$ 次迭代参数 $\theta$ 的估计值为 $\theta_t$，在给定观测数据 $\mathbf{X}$ 和当前参数 $\theta_t$，求得完全数据的对数似然函数关于未观测数据 $\mathbf{z}$ 的概率分布的期望 $\displaystyle \sum_{\mathbf{z}}\mathcal{P}(\mathbf{z}|\mathbf{X}, \theta_t)\ln\mathcal{P}(\mathbf{X}, \mathbf{z}|\theta)$，记为 $Q(\theta,\theta_t)$。
  3. （M 步）求使 $l(\theta|\theta_t)$ 最大化的 $\theta$，使 $\theta_{t+1}=\arg\max_{\theta} Q(\theta,\theta_t)$。
  4. 重复 E 步与 M 步。
  5. 如果满足 $|\theta_{t+1}-\theta_t|<\varepsilon_1$ 或 $|Q(\theta_{t+1},\theta_t)-Q(\theta_t,\theta_t)|<\varepsilon_1$，则停止迭代。

参考

1. The Expectation Maximization Algorithm: A short tutorial
2. Expectation–maximization algorithm - Wikipedia
3. Lecture 10: Expectation-Maximization Algorithm
4. EM算法原理总结 - 刘建平

---

Machine Learning

2022-04-08