个人学习笔记（十一）EM算法及其推广

       当模型含有隐变量时，无法使用极大似然估计法或贝叶斯估计法估计模型参数，这时适合使用EM算法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。EM算法的每次迭代由两步组成：E步，求期望(expectation)；M步，求极大(maximization)。所以这一算法称为期望极大算法(expectation maximization algorithm)，简称EM算法。

1、EM算法的引入

       EM算法通过迭代求$L(\theta )=logP(Y|\theta )$的极大似然估计。在介绍EM算法之前，需要先定义一个$Q$函数($Q$ function)，这是EM算法的核心。
       完全数据的对数似然函数$logP(Y,Z|\theta )$关于条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为$Q$函数，即
$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]$$       下面是EM算法的计算过程：
       首先选择参数的初值${\theta }^{(0)}$，进入下面的迭代。记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步中，计算$Q$函数
$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$       在M步中，求使$Q$函数极大化的$\theta$，作为第$i+1$次迭代的参数估计值${\theta }^{(i+1)}$
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$       重复上面的步骤直到收敛。
       需要注意的是，EM算法对初值是敏感的，每次迭代实际在求$Q$函数（E步）及其极大（M步）。
       那么，为什么上面的算法能够求解模型参数呢？该如何理解它呢？我们可以用近似求解对数似然函数的极大化问题来导出EM算法。
       首先，我们面对一个含有隐变量的概率模型，目标是极大化不完全数据$Y$关于参数$\theta$的对数似然函数，即极大化
$$L(\theta )=logP(Y|\theta )$$       由于概率模型含有隐变量$Z$，因此似然函数可写成
$$P(Y|\theta )=\sum _{Z}P(Y,Z|\theta )=\sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$       则对数似然函数便为
$$L(\theta )=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$       为了极大化上面的对数似然函数，可以考虑迭代方法，保证每一代所得的$\theta$能逐步逼近$L(\theta )$的最大值，即在第$i$代时，需要求解新的$\theta$使$L(\theta )>L({\theta }^{(i)})$。为此，考虑两者的差
$$L(\theta )-L({\theta }^{(i)})=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )-logP(Y|{\theta }^{(i)})$$       将$L(\theta )$变换一下，令
$$L(\theta )=log\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}$$       用Jensen不等式将$log$移进${\sum }_Z$内的系数右侧，公式为
$$log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log{y}_j,其中{\lambda }_j\ge 0,\sum _j{\lambda }_j=1$$       得
$$L(\theta )\ge \sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}$$       由于$L({\theta }^{(i)})$的公式与$Z$无关，可写作
$$L({\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|{\theta }^{(i)})$$       两式相减，得
$$L(\theta )-L({\theta }^{(i)})\ge \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)})}$$       为表述方便，引入$B(\theta ,{\theta }^{(i)})$的表达式
$$B(\theta ,{\theta }^{(i)})=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)})}$$       因此有
$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$       同时，当$\theta ={\theta }^{(i)}$时，有
$$B({\theta }^{(i)},{\theta }^{(i)})=L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|{\theta }^{(i)})}{P(Y,Z|{\theta }^{(i)})}=L({\theta }^{(i)})$$       因此，任何可以使$B(\theta ,{\theta }^{(i)})$增大的$\theta$，都能使$L(\theta )$增大，为此我们可以选择令$B(\theta ,{\theta }^{(i)})$最大的$\theta$作为${\theta }^{(i+1)}$，即
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})$$       将$B(\theta ,{\theta }^{(i)})$中与自变量$\theta$无关的项去除后，可得
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|Z,\theta )P(Z|\theta )$$       合并$P(Y|Z,\theta )P(Z|\theta )$得
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$       这正是$Q$函数$Q(\theta ,{\theta }^{(i)})$的表达式。
       因此，在每一次迭代中极大化$Q$函数，等价于极大化$B(\theta ,{\theta }^{(i)})$，也就是求取对数似然函数$L(\theta )$下界的极大值。这样，EM算法在迭代中便保证了对数似然函数$L(\theta )$的增加，但无法保证最后找到全局最优值。

2、EM算法的收敛性

       关于EM算法的收敛性，有两个定理。
       定理1：$P(Y|{\theta }^{(i)})$是单调递增的，即
$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$$       下面证明这一定理。由于对数似然函数$logP(Y|\theta )$与$Z$无关，可变换其形态，令
$$logP(Y|\theta )=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|\theta )=\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )}$$       因此，对数似然函数为
$$logP(Y|\theta )=\sum _{Z}P(Z|Y,{\theta }^{(i)})log[P(Y,Z|\theta )-P(Z|Y,\theta )]$$       令
$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$$$H(\theta ,{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )$$       则对数似然函数可写作
$$logP(Y|\theta )=Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})$$       接着便可开始证明下式成立
$$logP(Y|{\theta }^{(i+1)})\ge logP(Y|{\theta }^{(i)})$$       由于EM算法的每次迭代均为求取$Q$函数的极大值，因此必有
$$Q({\theta }^{(i+1)},{\theta }^{(i)})\ge Q({\theta }^{(i)},{\theta }^{(i)})$$       接着看$H$函数
$$H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}$$       根据Jensen不等式
$$H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})\le log\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}$$       因此可得
$$H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})\le log\sum _{Z}P(Z|Y,{\theta }^{(i+2)})=log1=0$$       这样，便证明了
$$logP(Y|{\theta }^{(i+1)})-logP(Y|\theta (i))\ge 0$$       即
$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$$       接下来是第二个定理。定理2：如果$P(Y|\theta )$有上界，则$L({\theta }^{(i)})=logP(Y|{\theta }^{(i)})$收敛到某一值$L^{\ast }$，在一定条件下EM算法得到的${\theta }^{(i)}$收敛值${\theta }^{\ast }$是$L(\theta )$的稳定点。
       这个定理就不证明了。定理只能保证参数估计序列${\theta }^{(i)}$收敛到对数似然函数序列$L({\theta }^{(i)})$的稳定点，不能保证收敛到极大值点。

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

       最后两节只叙述，不证明。
       首先介绍高斯混合模型(Gaussian misture model)，高斯混合模型是指具有如下概率分布的模型
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$       其中，${\alpha }_k$是系数，满足${\alpha }_k\ge 0,{\sum }_{k=1}^K{\alpha }_k=1$；$\varphi (y|{\theta }_k)$是高斯分布密度，${\theta }_k=({\mu }_k,{\sigma }_k^2)$，具体为
$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi {\sigma }_k}}{e}^{-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}}$$       接着介绍高斯混合模型参数估计的EM算法。有$N$个观测数据及$K$个分模型，首先取参数的初始值开始迭代。
       在E步，计算分模型$k$对观测数据$y_j$的响应度
$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)},j=1,2,\cdots ,N;k=1,2,\cdots ,K$$       在M步，计算新一轮的模型参数
$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N\widehat{{\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-{\hat{\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、EM算法的推广

       EM算法还可以解释为$F$函数($F$ function)的极大-极大算法(maximization-maximization algorithm)，基于这个解释有广义期望极大(generalized expectation maximization, GEM)算法等推广。
       首先引进$F$函数。假设隐变量$Z$的概率分布为$\tilde{P}(Z)$，定义$F$函数如下
$$F(\tilde{P},\theta )=E_{\tilde{P}}[logP(Y,Z|\theta )]+H(\tilde{P})$$       式中，$H(\tilde{P})$是分布$\tilde{P}(Z)$的熵，即
$$H(\tilde{P})=-E_{\tilde{P}}log\tilde{P}(Z)$$       引进$F$函数的概念后，便可给出下面四条定理。
       定理1：对于固定的$\theta$，存在唯一的分布${\tilde{P}}_{\theta }$使$F(\tilde{P},\theta )$极大化，这时的${\tilde{P}}_{\theta }$为
$${\tilde{P}}_{\theta }(Z)=P(Z|Y,\theta )$$       引理2：若${\tilde{P}}_{\theta }(Z)=P(Z|Y,\theta )$，则$F$函数为
$$F(\tilde{P},\theta )=logP(Y|\theta )$$       定理3：如果$F(\tilde{P},\theta )$在${\tilde{P}}^{\ast }$和${\theta }^{\ast }$有局部极大值，那么$L(\theta )$也在${\theta }^{\ast }$有局部极大值，这个结论对全局最大值也使用。
       定理4：EM算法的一次迭代，可由$F$函数的极大-极大算法实现。在第$i+1$次迭代的两步为：对固定的${\theta }^{(i)}$，求${\tilde{P}}^{(i+1)}$使$F(\tilde{P},{\theta }^{(i)})$极大化；对固定的${\tilde{P}}^{(i+1)}$，求${\theta }^{(i+1)}$使$F({\tilde{P}}^{(i+1)},\theta )$极大化。
       从上面的定理可知，由EM算法与F函数的极大-极大算法得到的${\theta }^{(i)}$是一致的。这样便有了EM算法的推广，GEM算法。
       在GEM算法1中，先初始化参数${\theta }^{(0)}$，在第$i+1$次迭代先求${\tilde{P}}^{(i+1)}$使$F(\tilde{P},{\theta }^{(i)})$极大化，再求${\theta }^{(i+1)}$使$F({\tilde{P}}^{(i+1)},\theta )$极大化。
       在GME算法2中，考虑到有时候求${\theta }^{(i+1)}$使$Q(\theta ,{\theta }^{(i)})$极大化比较困难，于是在原本EM算法的基础上放低要求，在M步中求取${\theta }^{(i+1)}$使$Q({\theta }^{(i+1)},{\theta }^{(i)})\ge Q({\theta }^{(i)},{\theta }^{(i)})$。
       在GME算法3中，当参数$\theta$的维数为$d$时，将EM算法的M步分解为$d$次条件极大化，每次只改变参数向量的一个分量。