一般形式的EM算法

问题定义

先给出一般EM算法的问题定义。 给定一个训练集\(X=\{x^{\{1\}},\ldots,x^{\{m\}}\}\)。根据模型的假设,我们希望能够通过这些数据来估计参数\(\theta\),但是每个我们观察到的\(x^{i}\) 还受着一个我们观察不到的隐含变量\(z^{i}\)(也是\(\theta\)生成的)控制,我们记\(Z=\{z^{\{1\}},\ldots,z^{\{m\}}\}\), 整个模型\(p(X,Z|\theta)\)是个关于变量\(X,Z\)的联合分布。

我们通过求最大似然\(L(\theta|X)\)来估计\(\theta\)的值。

\[ \theta=\arg\max_{\theta} L(\theta|X)\]

其中

\[ L(\theta|X)=\ln p(X|\theta)=\ln \sum_{Z}{p(X,Z|\theta)}\]

这里要注意\(p(X|\theta)\)是一个边缘分布,需要通过聚合\(Z\)得到。由于右边的式子在对数函数中又存在加和,所以我们直接对似然求导是没办法得到\(\theta\)的。

这里EM算法可以比较好地解决这个问题。对于EM的一般算法,有从期望的角度解释,用Jensen不等式证明正确性的方法(见参考2)。这里我讲的是另一种从相对熵角度解释,来自PRML的方法。这两种方法各有特定,对比如下:

  • 使用期望的角度来讲,可以地了解EM算法是如何绕过未知\(Z\)下难以计算的问题,即通过最大化它的期望。但是对E和M步如何逼近极大似然的过程,需要对Jensen不等式和单调逼近\(\theta\)最优值过程的理解。

  • 使用相对熵角度解释,可以看清似然函数和它的下界函的迭代增长过程。但是对\(Q(Z)\)和KL散度的引入,以及为什么要最大化下界函数这一点讲的不是很明白(个人看法)。2013年9月25日更新:这里应该是设一个分布Q(z)来近似P(Z),KL散度描述了这两个分布之间的差异,这里只是没有提到下界函数刚好是似然关于后验概率的期望,要求的是下界函数也就是期望的最大值。

推导过程

可以看到,求联合分布\(p(X,Z|\theta)\)的极值比较简单。我们根据贝叶斯公式,对\(\ln p(X|\theta)\)变形:

\[\begin{split}\ln p(X|\theta) =& \ln \frac{p(X,Z|\theta)}{p(Z|X,\theta)} \\ =& \ln p(X,Z|\theta)-\ln p(Z|X,\theta) \\ \end{split}\] 这里我们假设存在一个随机变量\(Z\)的概率分布\(q(Z)\),并有技巧地插入式子: \[\begin{split}\ln p(X|\theta)=& \ln p(X,Z|\theta)-\ln p(Z|X,\theta) \\=&\ln p(X,Z|\theta)-\ln q(Z)-[\ln p(Z|X,\theta)-\ln q(Z)]\\=& \ln \frac{ p(X,Z|\theta)}{q(Z)}-\ln \frac{ p(Z|X,\theta)}{q(Z)} \\\end{split}\]

等式两边同时乘以\(q(Z)\),并对\(Z\)求和

\[\sum_{z}q(Z) \ln p(X|\theta)=\sum_{z}[q(Z)(\ln \frac{ p(X,Z|\theta)}{q(Z)}-\ln \frac{ p(Z|X,\theta)}{q(Z)})]\]

因为\(Z\)\(p(X|\theta)\)独立且\(\sum_{z}q(Z)=1\),所以

\[ \ln p(X|\theta)=\sum_{z}q(Z)\ln \frac{ p(X,Z|\theta)}{q(Z)}-\sum_{z}q(Z)\ln \frac{ p(Z|X,\theta)}{q(Z)}\]

我们将其写作

\[\ln p(X|\theta)=\mathcal{L}(q,\theta)+KL(q||p)\]

其中:

\[\mathcal{L}(q,\theta)=\sum_{z}q(Z)\ln \frac{ p(X,Z|\theta)}{q(Z)}\]

\[ KL(q||p)=\sum_{z}q(Z)\ln \frac{ q(Z)}{p(Z|X,\theta)}\]

这里\(KL(q||p)\)(注1)描述了概率分布\(q(Z)\)\(p(Z|X,\theta)\)之间的差别。由于\(KL(q||p)\geq 0\),而且仅在\(q(Z)=p(Z|X,\theta)\)时等于0,所以总有\(\mathcal{L}(q,\theta)\leq \ln p(X|\theta)\)成立,我们把\(\mathcal{L}(q,\theta)\) 看做\(\ln p(X|\theta)\) 的下界函数。

假设参数当前的值为\(\theta^{old}\)

在E-step中,我们已知固定值\(\theta^{old}\),通过寻找\(q(Z)\)来最大化下界\(\mathcal{L}(q,\theta^{old})\) 的值。这个最大化问题比较简单, \(\ln p(X|\theta^{old})\)是一个定值,变量只有\(Z\),最大值在\(KL(q||p)=0\),也就是\(q(Z)= p(Z|X,\theta)\) 时取得(这个同时也解释了为什么要用\(Z\)的后验概率去而不是别的概率求Z 的期望值),此时\(\mathcal{L}(q,\theta) = \ln p(X|\theta^{old})\)(如图1所示)。

图1

在M-step中,我们固定\(q(Z)\),通过寻找新的参数\(\theta^{new}\)来最大化\(\mathcal{L}(q,\theta)\) 的值。通过增大下界函数\(\mathcal{L}(q,\theta)\)的值,来增大似然函数的值。因为概率分布\(q(Z)\)是用\(\theta^{old}\)决定的,当\(\theta\)改变之后,它将不再等于新的后验概率\(p(Z|X,\theta^{new})\)。因此\(KL\)散度也大于0.因此似然函数的增大将来自两个方面:下界函数的增大和KL散度的增大(如图2所示)。

图2

然后我们不断运行E步和M步直到算法收敛。如图3所示,红色曲线(未知)是我们想要最大化的似然函数。在第一个E步,我们初始化参数\(\theta^{old}\),然后通过它求得下界\(\mathcal{L}(q,\theta^{old})\)(蓝色曲线)。可以看到,下界函数在\(\theta^{old}\)时,为似然函数的正切线,所以两条曲线的梯度相同。注意到下界函数是个凸函数(参照国外convex function 定义,国内定义正好相反),所以它有唯一的最大值。在M步中,我们通过\({\theta^{new}}\)来获取下界函数的最大值。然后在接下来的E 步,我们重新构造一个下界函数(如绿色曲线所示),使得它在在\(\theta^{new}\) 时,为似然函数的正切线。

图3

此外,我们不仅可以用EM算法估计最大似然值,还能求最大后验概率(MAP)。在给定先验概率\(p(\theta)\) 的情况下,后验概率 \[ p(\theta|X)=\frac{p(X|\theta)p(\theta)}{p(X)}=\frac{p(\theta,X)}{p(X)}\] 所以我们有 \[\begin{split}\ln p(\theta|X)=&\ln p(\theta,X)-\ln p(X) \\ =&\ln p(X|\theta)+\ln p(\theta)-\ln p(X) \\ =&\mathcal{L}(q,\theta)+KL(q||p) +\ln p(\theta)-\ln p(X) \\\end{split}\]

这里\(p(X)\)为常数。

EM算法步骤

这里的EM算法:

  1. E-step:固定\(\theta^{old}\),那么\(\ln p(\theta^{old})\)为常数,找一个分布\(q(Z)\),使\(\mathcal{L}(q,\theta^{old})\)最大。
  2. M-step:固定\(q(Z)\),寻找参数\(\theta^{new}\),使\(\mathcal{L}(q,\theta^{new}) +\ln p(\theta^{new})\)最大。
  3. 重复1.2直至收敛。

除此之外,如果 M-step中的最大化\(\mathcal{L}(q,\theta)\)\(\theta^{new}\)也很难求的话,我们可以把要求放宽,并不一定要求出最大值,只要能够得到比\(\theta^{old}\)更好的结果就可以,这样的做法叫 Generalized EM (GEM)。

备注

KL散度,在信息论中称为相对熵,用于衡量两个函数或概率分布之间的差异程度。性质如下:

  • \(KL(q||p)\neq KL(p||q)\)>
  • \(KL(q||p)\)越大,差异程度越大。当\(p\)\(q\)完全相同时,\(KL(p||q)=0\)

参考:

  1. 《Pattern_Recognition_and_Machine_Learning》第九章
  2. http://blog.tomtung.com/2011/10/em-algorithm/