EM算法及其推广

EM算法是一种迭代算法。用于含有隐变量(hidden variable)的概率模型参数的极大似然估计，或极大后验概率估计。

EM算法每次迭代由两步组成:

• E步，求期望
• M步，求极大

因此，EM算法也称为期望极大算法。

本章首先叙述EM算法，然后讨论算法的收敛性，接着介绍EM算法中的高斯混合模型的学习，最后叙述EM算法的推广——GEM算法。

12.1 EM算法的引入

概率模型有时既含有观测变量，又含有隐变量或潜在变量。

如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法或贝叶斯估计法估计模型参数。

如果概率模型含有隐变量，则不能简单地使用以上估计方法。

EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

在这里，我们只讨论极大似然估计法。

12.1.1 EM算法

一般地，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z连在一起称为完全数据，观测数据Y又称为不完全数据。

假设给定观测数据Y，其概率分布是$P(Y|\theta)$，其中$\theta$是需要估计的模型参数，那么不完全数据Y的似然函数是$P(Y|\theta)$，对数似然函数$L(\theta)=\log P(Y|\theta)$；假设Y和Z的联合概率分布是$P(Y,Z|\theta)$，那么完全数据的对数似然函数是$\log P(Y,Z|\theta)$。

EM算法通过迭代求$L(\theta)=\log P(Y|\theta)$的极大似然估计。每次迭代包含两步：E步，求期望；M步，求极大化。

算法12.1(EM算法)

输入：观测变量数据Y，隐变量数据Z，联合分布$P(Y,Z|\theta)$，条件分布$P(Z|Y,\theta)$；

输出：模型参数$\theta$。

(1)选择参数的初值$\theta^{(0)}$，开始迭代；

(2)E步：记$\theta^{(i)}$为第i次迭代参数$\theta$的估计值，在第i+1次迭代的E步，计算：

$$\begin{aligned} Q(\theta,\theta^{(i)}) &= E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}] \\ &= \sum_Z \log P(Y,Z|\theta)P(Z|Y,\theta^{(i)}) \ (12.1)\end{aligned}$$

这里，$P(Z|Y,\theta^{(i)}$是在给定观测数据Y和当前的参数估计$\theta^{(i)}$下隐变量数据Z的条件概率分布；

(3)M步：求使$Q(\theta,\theta^{(i)})$极大化的$\theta$，确定第i+1次迭代的参数的估计值$\theta^{(i+1)}$

$$\theta^{(i+1)} = \arg \max_{\theta} Q(\theta,\theta^{(i)}) \ (12.2)$$

(4)重复第(2)步和第(3)步，直到收敛。

式(12.1)的函数$Q(\theta,\theta^{(i)})$是EM算法的核心，称为Q函数。

定义12.1(Q函数) 完全数据的对数似然函数$\log P(Y,Z|\theta)$关于在给定观测数据Y和当前参数$\theta^{(i)}$下对未观测数据Z的条件概率分布$P(Z|Y,\theta^{(i)})$的期望称为Q函数，即

$$Q(\theta,\theta^{(i)}) = E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}] \ (12.3)$$

下面关于EM算法作几点说明：

步骤(1) 参数的初值可以任意选择，但需注意EM算法对初值是敏感的。

步骤(2) E步求$Q(\theta,\theta^{(i)})$。Q函数式中Z是未观测数据，Y是观测数据。(注意：$Q(\theta,\theta^{(i)})$的第1个变量表示要极大化的参数，第2个变量表示参数的当前估计值)。每次迭代实际在求Q函数及其极大。

步骤(3) M步求$Q(\theta,\theta^{(i)})$的极大化，得到$\theta^{(i+1)}$，完成一次迭代$\theta^{(i)} \rightarrow \theta^{(i+1)}$，后面将证明每次迭代使似然函数增大或达到局部极值。

步骤(4) 给出停止迭代的条件，一般是对较小的正数$\varepsilon_1$，$\varepsilon_2$，若满足

$$\|\theta^{(i+1)}-\theta^{(i)}\| < \varepsilon_1$$

或者

$$\|Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})\|<\varepsilon_2$$

则停止迭代。

12.1.2 EM算法的导出

为什么EM算法能够近似实现对观测数据的极大似然估计呢？

下面通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法，由此可以清楚地看出EM算法的作用。

对于一个含有隐变量的概率模型，目标是极大化观测数据(不完全数据)Y关于参数$\theta$的对数似然函数，即极大化

$$\begin{aligned} L(\theta) &= \log P(Y|\theta) = \log \sum_Z P(Y,Z|\theta)\\&= \log \left(\sum_Z P(Y|Z,\theta)P(Z|\theta) \right) \end{aligned} \ (12.4)$$

式(12.4)中有未观测数据并有包含和(或积分)的对数。

事实上，EM算法是通过迭代逐步近似极大化$L(\theta)$的。假设在第i次迭代后$\theta$的估计值是$\theta^{(i)}$。我们希望新估计值$\theta$能使$L(\theta)$增加，即$L(\theta)>L(\theta^{(i)})$，并逐步达到极大值。为此，考虑两者的差:

$$L(\theta)-L(\theta^{(i)}) = \log \left(\sum_Z P(Y|Z,\theta)P(Z|\theta)\right)-\log P(Y|\theta^{(i)})$$

利用Jensen不等式得到其下界:

Jensen不等式：$\log \sum_j \lambda_j y_j \ge \sum_j \lambda_j \log y_j$，其中$\lambda_j \ge 0$，$\sum_j \lambda_j = 1$。

如下所示：使用jensen不等式得到第二个式子，其实相当于计算$P(Z|Y,\theta^{(i)})$与$P(Y|Z,\theta)P(Z|\theta)$的相对熵或者说时两个分布KL散度。

相对熵:两个概率密度函数p(x)和q(x)之间的相对熵定义为:

$$D_{KL}(p\|q) = \sum_{x \in \mathcal{X}} p(x) \log \frac{p(x)}{q(x)}$$

$$\begin{aligned}L(\theta)-L(\theta^{(i)}) &= \log \left(\sum_Z P(Z|Y,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}\right)-\log P(Y|\theta^{(i)}) \\ &\ge \sum_Z P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-\log P(Y|\theta^{(i)}) \\ &= \sum_Z P(Z|Y,\theta^{(i)}) \log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\end{aligned}$$

令

$$B(\theta,\theta^{(i)}) \hat{=} L(\theta^{(i)})+\sum_Z P(Z|Y,\theta^{(i)}) \log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})} \ (12.5)$$

则

$$L(\theta) \ge B(\theta,\theta^{(i)}) \ (12.6)$$

即函数$B(\theta,\theta^{(i)})$是$L(\theta)$的一个下界，而且由式(12.5)可知，

$$L(\theta^{(i)}) = B(\theta^{(i)},\theta^{(i)}) \ (12.7)$$

因此，任何可以使$B(\theta,\theta^{(i)})$增大的$\theta$，也可以使$L(\theta)$增大，为了使$L(\theta)$有尽可能大的增长，选择$\theta^{(i+1)}$使$B(\theta,\theta^{(i)})$达到极大，即

$$\theta^{(i+1)} = \arg \max_{\theta} B(\theta,\theta^{(i)}) \ (12.8)$$

现在求$\theta^{(i+1)}$的表达式，省去对$\theta$的极大化而言是常数的项，由式(12.8)、式(12.5)及式(12.2)，有

$$\begin{aligned}\theta^{(i+1)}&=\arg \max_{\theta} \left(L(\theta^{(i)}+\sum_Z P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\right) \\ &= \arg \max_{\theta} \left(\sum_Z P(Z|Y,\theta^{(i)})\log(P(Y|Z,\theta)P(Z|\theta))\right)\\ &= \arg \max_{\theta} \left(\sum_Z P(Z|Y,\theta^{(i)}) \log P(Y,Z|\theta)\right) \\ &= \arg \max_{\theta} Q(\theta,\theta^{(i)})\end{aligned} \ (12.9)$$

式(12.9)等价于EM算法的一次迭代，即求Q函数及其极大化。EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

下图给出EM算法的直观解释。图中上方曲线为$L(\theta)$，下方曲线为$B(\theta,\theta^{(i)})$。

根据式(12.6)，$B(\theta,\theta^{(i)})$为对数似然函数$L(\theta)$的下界。

根据式(12.7)，两个函数在点$\theta=\theta^{(i)}$处相等。

根据式(12.8)和式(12.9)，EM算法找到下一个点$\theta^{(i+1)}$使函数$B(\theta,\theta^{(i)})$极大化，也使函数$Q(\theta,\theta^{(i)})$极大化。

由于$L(\theta) \ge B(\theta,\theta^{(i)})$，函数$B(\theta,\theta^{(i)})$的增加，保证对数似然函数$L(\theta)$在每次迭代中也是增加的。

EM算法在点$\theta^{(i+1)}$重新计算Q函数值，进行下一次迭代。在这个过程中，对数似然函数$L(\theta)$不断增大。从图可以推断出EM算法不能保证找到全局最优值。

12.1.3 EM算法在非监督学习中的应用

当训练数据只有输入没有相应的输出标签，如$\{(x_1,\bullet),(x_2,\bullet),\cdots,(x_N,\bullet)\}$，从这样的数据学习模型称为非监督学习问题。

EM算法可以用于生成模型的非监督学习。生成模型由联合概率分布P(X,Y)表示，可以认为非监督学习训练数据是联合概率分布产生的数据。X为观测数据，Y为未观测数据。

12.2 EM算法的收敛性

EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。对于EM算法，需要回答两个问题：EM算法得到的估计序列是否收敛？如果收敛，是否收敛到全局最大值或局部极大值？

下面给出EM算法收敛性相关的两个定理：

定理12.1 设$P(Y|\theta)$为观测数据的似然函数，$\theta^{(i)}(i=1,2,\cdots)$为EM算法得到的参数估计序列，$P(Y|\theta^{(i)})(i=1,2,\cdots)$为对应的似然函数序列，则$P(Y|\theta^{(i)})$是单调递增的，即

$$P(Y|\theta^{(i+1)}) \ge P(Y|\theta^{(i)}) \ (12.10)$$

证明 由于

$$P(Y|\theta) = \frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}$$

取对数有

$$\log P(Y|\theta) = \log P(Y,Z|\theta)-\log P(Z|Y,\theta)$$

由式(12.3)

$$Q(\theta,\theta^{(i)}) = \sum_Z \log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})$$

令

$$H(\theta,\theta^{(i)}) = \sum_Z \log P(Z|Y,\theta)P(Z|Y,\theta^{(i)}) \ (12.11)$$

于是对数似然函数可以写成

$$\log P(Y|\theta) = Q(\theta,\theta^{(i)})-H(\theta,\theta^{(i)}) \ (12.12)$$

在式(12.12)中分别取$\theta$为$\theta^{(i)}$和$\theta^{(i+1)}$并相减，有

$$\begin{aligned}\log &P(Y|\theta^{(i+1)})-\log P(Y|\theta^{(i)}) \\&=[Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})]-[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)})] \end{aligned} \ (12.13)$$

为证式(12.10)，只需证式(12.13)右端是非负的。式(12.13)右端的第1项，由于$\theta^{(i+1)}$使$Q(\theta,\theta^{(i)})$达到极大，所以有：

$$Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)}) \ge 0 \ (12.14)$$

其第2项，由式(12.11)可知:

$$\begin{aligned}H(\theta^{(i+1)}&,\theta^{(i)})-H(\theta^{(i)},\theta^{(i)}) \\ &= \sum_Z \left(\log \frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}\right)P(Z|Y,\theta^{(i)}) \\ &\le \log \left(\sum_Z \frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})} P(Z|Y,\theta^{(i)})\right) \\ &= \log \left(\sum_Z P(Z|Y,\theta^{(i+1)})\right) = 0\end{aligned} \ (12.15)$$

注意:这里的不等号由Jensen不等式得到。

由式(12.14)和式(12.15)可知式(12.13)右端是非负的。

定理12.2 设$L(\theta)=\log P(Y|\theta)$为观测数据的对数似然函数，$\theta^{(i)}(i=1,2,\cdots)$为EM算法得到的参数估计序列，$L(\theta^{(i)})(i=1,2,\cdots)$为对应的对数似然函数序列。

(1)如果$P(Y|\theta)$有上界，则$L(\theta^{(i)})=\log P(Y|\theta^{(i)})$收敛到某一值$L^*$

(2)在函数$Q(\theta,\theta^{'})$与$L(\theta)$满足一定条件下，由EM算法得到的参数估计序列$\theta^{(i)}$的收敛值$\theta^{*}$是$L(\theta)$的稳定点。

证明

(1)由$L(\theta)=\log P(Y|\theta^{(i)})$的单调性及$P(Y|\theta)$的有界性立即得到。

(2)证明从略，参考文献:A View of the Em Algorithm that Justifies Incremental, Sparse, and other Variants

EM算法的收敛性包含关于对数似然函数序列$L(\theta^{(i)})$的收敛性和关于参数估计序列$\theta^{(i)}$的收敛性两层意思。

在应用中，初值的选择变得非常重要，常用的方法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

12.3 EM算法在高斯混合模型学习中的应用

在许多情况下，EM算法是学习高斯混合模型的有效方法。

12.3.1 高斯混合模型

定义12.2(高斯混合模型) 高斯混合模型是指具有如下形式的概率分布模型:

$$P(y|\theta) = \sum_{k=1}^K \alpha_k \phi(y|\theta_k) \ (12.16)$$

其中，$\alpha_k$是系数，$\alpha_k \ge 0$，$\sum_{k=1}^K \alpha_k = 1$；$\phi(y|\theta_k)$是高斯分布密度，$\theta_k = (\mu_k, \sigma_k^2)$，

$$\phi(y|\theta_k) = \frac{1}{\sqrt{2\pi}\sigma_k} \exp \left(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\right) \ (12.17)$$

称为第k个分模型。

一般混合模型可以由任意概率分布密度代替式(12.17)中的高斯分布密度，我们只介绍最常见的高斯混合模型。

12.3.2 高斯混合模型参数估计的EM算法

假设观测数据$y_1,y_2,\cdots,y_N$由高斯混合模型生成，

$$P(y|\theta) = \sum_{k=1}^K \alpha_k \phi(y|\theta_k) \ (12.18)$$

其中，$\theta=(\alpha_1,\alpha_2,\cdots,\alpha_K;\theta_1,\theta_2,\cdots,\theta_K)$。我们用EM算法估计高斯混合模型的参数$\theta$。

1. 明确隐变量，写出完全数据的对数似然函数

观测数据$y_j,j=1,2,\cdots,N$，是这样产生的:

1. 首先依概率$\alpha_k$选择第k个高斯分布分模型$\phi(y|\theta_k)$；

2. 然后依第k个分模型的概率分布$\phi(y|\theta_k)$生成观测数据$y_j$。这时，观测数据$y_j,j=1,2,\cdots,N$，是已知的；

   反映观测数据$y_j$来自第k个分模型的数据是未知的，k=1,2,…,K，以隐变量$\gamma_{jk}$表示，其定义如下:

   $$\gamma_{jk}=\left\{\begin{aligned}&1,\mathrm{j-th\ sample \ come \ from \ k-th \ model} \\ &0, \mathrm{others}\end{aligned}\right. \ (12.19)$$ $$j=1,2,\cdots,N;k=1,2,\cdots,K$$

   $\gamma_{jk}$是0-1随机变量。

有了观测数据$y_j$及未观测数据$\gamma_{jk}$，那么完全数据是

$$(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}) , j=1,2,\cdots,N$$

于是，可以写出完全数据的似然函数:

$$\begin{aligned}P(y,\gamma|\theta) &= \prod_{j=1}^N P(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}|\theta) \\ &= \prod_{k=1}^K \prod_{j=1}^N [\alpha_k \phi(y_j|\theta_k)]^{\gamma_{jk}} \\ &= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N [\phi(y_j|\theta_k)]^{\gamma_{jk}} \\ &=\prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N \left[\frac{1}{\sqrt{2\pi}\sigma_k} \exp \left(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\right)\right]^{\gamma_{jk}}\end{aligned}$$

式中，$n_k = \sum_{j=1}^N \gamma_{jk}$，$\sum_{k=1}^K n_k = N$。

那么，完全数据的对数似然函数为:

$$\log P(y,\gamma|\theta) = \sum_{k=1}^K \left\{n_k \log \alpha_k + \sum_{j=1}^N \gamma_{jk} \left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2\right]\right\}$$

2.EM算法的E步:确定Q函数

$$\begin{aligned}Q(\theta,\theta^{(i)}) &= E[\log P(y,\gamma|\theta)|y,\theta^{(i)}] \\ &= E \left\{\sum_{k=1}^K \left\{n_k \log \alpha_k + \sum_{j=1}^N \gamma_{jk} \left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2\right]\right\}\right\}\\&=\sum_{k=1}^K\left\{\sum_{j=1}^N(E\gamma_{jk})\log \alpha_k+\sum_{j=1}^N(E\gamma_{jk})\left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2\right]\right\}\end{aligned} \ (12.20)$$

这里需要计算$E(\gamma_{jk}|y,\theta)$，记为$\hat{\gamma}_{jk}$。

$$\begin{aligned}\hat{\gamma}_{jk} &= E(\gamma_{jk}|y,\theta)=P(\gamma_{jk}=1|y,\theta) \\ &= \frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^K P(\gamma_{jk}=1,y_j|\theta)}\\&=\frac{P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^K P(y_j|\gamma_{jk}=1,\theta)P(\gamma_jk=1|\theta)}\\&=\frac{\alpha_k \phi(y_j|\theta_k)}{\sum_{k=1}^K \alpha_k \phi(y_j|\theta_k)}\end{aligned}$$ $$j=1,2,\cdots,N;k=1,2,\cdots,K$$

关于$\hat{\gamma}_{jk}$是在当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据$y_j$的响应度。

将$\hat{\gamma}_{jk}=E\gamma_{jk}$及$n_k = \sum_{j=1}^N E\gamma_{jk}$代入式（12.20）即得

$$Q(\theta,\theta^{(i)}) = \sum_{k=1}^K \left\{n_k \log \alpha_k + \sum_{j=1}^N \hat{\gamma}_{jk} \left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2\right]\right\} \ (12.21)$$

3. 确定EM算法的M步

迭代的M步是求函数$Q(\theta,\theta^{(i)})$对$\theta$的极大值，即求新一轮迭代的模型参数:

$$\theta^{(i+1)} = \arg \max_{\theta} Q(\theta,\theta^{(i)})$$

用$\hat{\mu}_k$，$\hat{\sigma}_k^2$及$\hat{\alpha}_k$，k=1,2,…,K，表示$\theta^{(i+1)}$的各参数。求$\hat{\mu}_k$，$\hat{\sigma}_k^2$只需将式(12.21)分别对$\mu_k$，$\sigma_k^2$求偏导数并令其为0，即可得到；求$\hat{\alpha}_k$是在$\sum_{k=1}^K \alpha_k = 1$条件下求偏导数并令其为0得到的(此处使用拉格朗日乘子法求解)。结果如下：

$$\hat{\mu}_k = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}},k=1,2,\cdots,K \ (12.22)$$ $$\hat{\sigma}_k^2 = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}} , k=1,2,\cdots,K \ (12.23)$$ $$\hat{\alpha}_k = \frac{n_k}{N}=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N} , k=1,2,\cdots,K \ (12.24)$$

重复以上计算，直到对数似然函数值不再有明显的变化为止。

算法12.2(高斯混合模型参数估计的EM算法)

输入：观测数据$y_1,y_2,\cdots,y_N$,高斯混合模型

输出：高斯混合模型参数

(1)取参数的初始值开始迭代

(2)E步：依据当前模型参数，计算分模型k对观测数据$y_j$的响应度

$$\hat{\gamma}_{jk} = \frac{\alpha_k \phi(y_j|\theta_k)}{\sum_{k=1}^K \alpha_k \phi(y_j|\theta_k)} , j=1,2,\cdots,N;k=1,2,\cdots,K$$

(3)M步：计算新一轮迭代的模型参数

$$\hat{\mu}_k = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j}{\sum_{j=1}^N\hat{\gamma}_{jk}} , k=1,2,\cdots,K$$ $$\hat{\sigma}_k^2 = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}} , k=1,2,\cdots,K$$ $$\hat{\alpha}_k = \frac{\sum_{j=1}^N \hat{\gamma}_{jk}}{N},k=1,2,\cdots,K$$

(4)重复第(2)步和第(3)步，直到收敛。

12.4 EM算法的推广

EM算法还可以解释为F函数的极大-极大算法。

12.4.1 F函数的极大-极大算法

首先引进F函数并讨论其性质。

定义12.3(F函数) 假设隐变量数据Z的概率分布为$\tilde{P}(Z)$，定义分布$\tilde{P}$与参数$\theta$的函数$F(\tilde{P},\theta)$如下：

$$F(\tilde{P},\theta) = E_{\tilde{P}}[\log P(Y,Z|\theta)]+H(\tilde{P}) \ (12.25)$$

称为F函数。式中$H(\tilde{P})=-E_{\tilde{P}}\log \tilde{P}(Z)$是分布$\tilde{P}(Z)$的熵。

通常假设$P(Y,Z|\theta)$是$\theta$的连续函数，因而$F(\tilde{P},\theta)$是$\tilde{P}$和$\theta$的连续函数。

引理12.1 对于固定的$\theta$，存在唯一的分布$\tilde{P}_{\theta}$极大化$F(\tilde{P},\theta)$，这时$\tilde{P}_{\theta}$由下式给出:

$$\tilde{P}_{\theta}(Z) = P(Z|Y,\theta) \ (12.26)$$

并且$\tilde{P}_{\theta}$随$\theta$连续变化。

证明: 对于固定的$\theta$，可以求得使$F(\tilde{P},\theta)$达到极大的分布$\tilde{P}_{\theta}(Z)$。为此，引进拉格朗日乘子$\lambda$，拉格朗日函数为 :

$$L = E_{\tilde{P}} \log P(Y,Z|\theta)-E_{\tilde{P}}\log \tilde{P}(Z)+\lambda \left(1-\sum_{Z}\tilde{P}(Z)\right) \ (12.27)$$

将其对$\tilde{P}$求偏导数，并令偏导数等于0:

$$\begin{aligned}\frac{\partial L}{\partial \tilde{P}(Z)} &= \sum_Z \left(\log P(Y,Z|\theta)-\log \tilde{P}(Z) -1 -\lambda\right) \\&= \log P(Y,Z|\theta)-\log \tilde{P}(Z) -1 -\lambda \\ &= 0\end{aligned}$$

得出

$$\lambda = \log P(Y,Z|\theta)-\log \tilde{P}_{\theta}(Z)-1$$

由此推出$\tilde{P}_{\theta}(Z)$与$P(Y,Z|\theta)$成比例

$$\frac{P(Y,Z|\theta)}{\tilde{P}_{\theta}} = e^{1+\lambda}$$

再从约束条件$\sum_Z \tilde{P}_{\theta}=1$得式(12.26)(注意：为什么可以得到式(12.26)这个结论)

由假设$P(Y,Z|\theta)$是$\theta$的连续函数，得到$\tilde{P}_{\theta}$是$\theta$的连续函数。

引理12.2 若$\tilde{P}_{\theta}(Z)=P(Z|Y,\theta)$，则

$$F(\tilde{P},\theta) = \log P(Y|\theta) \ (12.28)$$

由以上引理，可以得到关于EM算法用F函数的极大-极大算法的解释。

定理12.3 设$L(\theta)=\log P(Y|\theta)$为观测数据的对数似然函数，$\theta^{(i)}$，i=1,2,…,为EM算法得到的参数估计序列，函数$F(\tilde{P},\theta)$由式(12.25)定义。

如果$F(\tilde{P},\theta)$在$\tilde{P}^*$和$\theta^{*}$有局部极大值，那么$L(\theta)$也在$\theta^{*}$有局部极大值。

如果$F(\tilde{P},\theta)$在$\tilde{P}^*$和$\theta^{*}$达到全局最大值，那么$L(\theta)$也在$\theta^{*}$达到全局最大值。

证明：由引理12.1和引理12.2可知，$L(\theta)=\log P(Y|\theta)=F(\tilde{P}_{\theta},\theta)$对任意$\theta$成立。特别地，对于使$F(\tilde{P},\theta)$，有

$$L(\theta^*) = F(\tilde{P}_{\theta^*},\theta^*)=F(\tilde{P}^*,\theta^*) \ (12.29)$$

为了证明$\theta^*$是$L(\theta)$的极大点，需要证明不存在接近$\theta^*$的点$\theta^{**}$，使$L(\theta^{**})>L(\theta^*)$。

假如存在这样的点$\theta^{**}$，那么应有$F(\tilde{P}^{**},\theta^{**})>F(\tilde{P}^*,\theta^*)$，这里$\tilde{P}^{**}=\tilde{P}_{\theta^{**}}$。但因$\tilde{P}_{\theta}$是随$\theta$连续变化的，$\tilde{P}^{**}$应接近$\tilde{P}^*$，这与$\tilde{P}^*$和$\theta^*$是$F(\tilde{P},\theta)$的局部极大点的假设矛盾。

定理12.4 EM算法的一次迭代可由F函数的极大-极大算法实现。

设$\theta^{(i)}$为第i次迭代参数$\theta$的估计，$\tilde{P}^{(i)}$为第i次迭代函数$\tilde{P}$的估计。在第i+1次迭代的两步为:

(1)对固定的$\theta^{(i)}$，求$\tilde{P}^{(i+1)}$使$F(\tilde{P},\theta^{(i)})$极大化。

(2)对固定的$\tilde{P}^{(i+1)}$，求$\theta^{(i+1)}$使$F(\tilde{P}^{(i+1)},\theta)$极大化

证明

(1)由引理12.1，对于固定的$\theta^{(i)}$，

$$\tilde{P}^{(i+1)}(Z) = \tilde{P}_{\theta^{(i)}}(Z) = P(Z|Y,\theta^{(i)})$$

使$F(\tilde{P},\theta^{(i)})$极大化。此时，

$$\begin{aligned}F(\tilde{P}^{(i+1),\theta}) &= E_{\tilde{P}^{(i+1)}}[\log P(Y,Z|\theta)]+H(\tilde{P}^{(i+1)})\\&=\sum_Z \log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})+H(\tilde{P}^{(i+1)})\end{aligned}$$

由$Q(\theta,\theta^{(i)})$的定义式(12.3)有

$$F(\tilde{P}^{(i+1)},\theta) = Q(\theta,\theta^{(i)})+H(\tilde{P}^{(i+1)})$$

(2)固定$\tilde{P}^{(i+1)}$，求$\theta^{(i+1)}$使$F(\tilde{P}^{(i+1)},\theta)$极大化。得到：

$$\theta^{(i+1)} = \arg \max_{\theta} F(\tilde{P}^{(i+1)},\theta) = \arg \max_{\theta} Q(\theta,\theta^{(i)})$$

通过以上两步完成了EM算法的一次迭代。

由此可知，由EM算法与F函数的极大-极大算法得到的参数估计序列$\theta^{(i)}$，i=1,2,…,是一致的。

12.4.2 GEM算法

算法12.3(GEM算法1)

输入：观测数据,F函数

输出：模型参数

(1)初始化参数$\theta^{(0)}$，开始迭代

(2)第i+1次迭代

第1步:记$\theta^{(i)}$为参数$\theta$的估计值，计算

$$\begin{aligned}Q(\theta,\theta^{(i)})&= E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]\\&= \sum_Z P(Z|Y,\theta^{(i)})\log P(Y,Z|\theta)\end{aligned}$$

(3)第2步:求$\theta^{(i+1)}$使：

$$Q(\theta^{(i+1)},\theta^{(i)}) > Q(\theta^{(i)},\theta^{(i)})$$

(4)重复(2)和(3)，直到收敛。

当参数$\theta$的维数为d$(d \ge 2)$时，可采用一种特殊的GEM算法，它将EM算法的M步分解为d次条件极大化，每次只改变参数向量的一个分量，其余分量不改变。

算法9.5(GEM算法3)

输入:观测数据，Q函数

输出:模型参数

(1) 初始化参数$\theta^{(0)}=(\theta_1^{(0)},\theta_2^{(0)},\cdots,\theta_d^{(0)})$，开始迭代

(2)第i+1次迭代，

第1步：记$\theta^{(i)}=(\theta_1^{(i)},\theta_2^{(i)},\cdots,\theta_d^{(i)})$为参数$\theta=(\theta_1,\theta_2,\cdots,\theta_d)$的估计值，计算

$$\begin{aligned}Q(\theta,\theta^{(i)}) &= E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}] \\ &= \sum_ZP(Z|y,\theta^{(i)})\log P(Y,Z|\theta)\end{aligned}$$

(3)第2步:进行d次条件极大化:

首先，在$\theta_2^{(i)},\cdots,\theta_k^{(i)}$保持不变的条件下求使$Q(\theta,\theta^{(i)})$达到极大的$\theta_1^{(i+1)}$

然后，在$\theta_1=\theta_1^{(i+1)}$，$\theta_j=\theta_j^{(i)}$，j=3,4,…,k的条件下求使$Q(\theta,\theta^{(i)})$达到极大的$\theta_2^{(i+1)}$；

如此继续，经过d次条件极大化，得到$\theta^{(i+1)}=(\theta_1^{(i+1)},\theta_2^{(i+1)},\cdots,\theta_d^{(i+1)})$使得

$$Q(\theta^{(i+1)},\theta^{(i)}) > Q(\theta^{(i)},\theta^{(i)})$$

(4)重复(2)和(3)。直到收敛。

12.5 相对熵角度理解EM

可以参考:PRML第9章

更多阅读资料

• EM算法）The EM Algorithm
• 《统计学习方法》第9章 EM/GMM/F-MM/GEM
• EM算法的九层境界：Hinton和Jordan理解的EM算法