EM算法原理和应用

      EM算法是带隐变量概率模型的推断算法。今天我们介绍 EM 算法的原理和应用。我们先介绍推导出 EM 算法的一般方法,再介绍另一种 EM 算法推导方法,最后将 EM 算法应用于高斯混合模型。

1.EM 算法的推导

      我们让 (\pmb{y}) 表示可见变量,(\pmb{x}) 表示隐变量,(\theta) 表示模型参数。EM算法的目标是得到一个参数(\theta),使得概率模型 (p(\pmb{y}|\theta)) 似然值极大。
\begin{eqnarray}
\theta &=& argmax_{\theta} p(\pmb{y}|\theta) \nonumber \
&=& argmax_{\theta}lnp(\pmb{y}|\theta) \nonumber \
&=& argmax_{\theta} ln\sum_{\pmb{x}}p(\pmb{x},\pmb{y}|\theta) \nonumber
\end{eqnarray}

在通常情况下, 直接优化 (lnp(\pmb{y}|\theta)) 很困难。我们希望找到 (lnp(\pmb{y}|\theta)) 的一个下界函数,然后极大这个下界。EM 算法处于第 i 步时,参数为 (\theta^i)。
\begin{eqnarray}
lnp(\pmb{y}|\theta) &=& ln\sum_{\pmb{x}}p(\pmb{x},\pmb{y}|\theta) \nonumber \
&=& ln\sum_{\pmb{x}} \frac{p(\pmb{x},\pmb{y}|\theta)}{p(\pmb{x}|\pmb{y},\theta^i)}p(\pmb{x}|\pmb{y},\theta^i) \nonumber \
&\ge& \sum_{\pmb{x}} p(\pmb{x}|\pmb{y},\theta^i) lnp(\pmb{x},\pmb{y}|\theta)+H(p(\pmb{x}|\pmb{y},\theta^i)) \quad ln是凹函数 \nonumber \
&\ge& E_{\pmb{x}}[lnp(\pmb{x},\pmb{y}|\theta)|\pmb{y},\theta^i] \nonumber
\label{eq:emsimple}
\end{eqnarray}
定义 Q 函数等于这个下界函数。
\begin{eqnarray}
Q(\theta,\theta^{i}) = E_{\pmb{x}}[lnp(\pmb{x},\pmb{y}|\theta)|\pmb{y},\theta^i] \nonumber
\end{eqnarray}

      有了 Q 函数,EM 算法交替如下 E 步和 M 步。

      E步,计算出 Q 函数;
      M步,求出\(\theta^{i+1}\)以极大化Q函数,即\(\theta^{i+1}=argmax_{\theta}Q(\theta,\theta^{i})\)。

2.EM 算法的另一种推导

2.1. 交替优化 F 函数

      假设隐变量的概率分布为 (\bar{p}(\pmb{x})), 我们定义 F 函数
\begin{eqnarray}
F(\theta,\bar{p}) = E_{\bar{p}}[lnp(\pmb{x},\pmb{y}|\theta)]+H(\bar{p}) \nonumber \
\end{eqnarray}
其中(p(\pmb{x},\pmb{y}|\theta))是(\theta)的连续函数。F函数是 (lnp(\pmb{y}|\theta)) 的下界函数。我们用交替优化算法 (Alternative Method) 优化 F 函数。交替优化方法适用优化两组参数,先固定一组优化另一组,然后固定刚刚优化的那组优化另一组,交替进行下去。故而 EM 算法交替执行如下 E 步和 M 步。

      E步,固定参数 \(\theta^i\), 求隐变量分布以极大化 F 函数,用拉格朗日法求得隐变量分布如下所示 (求解细节见文末附录)。
\begin{eqnarray}
\bar{p}_{\theta^i}(\pmb{x})=p(\pmb{x}|\theta^i,\pmb{y}) \nonumber \\
\end{eqnarray}
      M步,固定隐变量分布,极大化 F 函数求得 \(\theta\)
\begin{eqnarray}
\theta &=& argmax_{\theta} F(\theta,\bar{p}) \nonumber \\
&=& argmax_{\theta}\{E_{\pmb{x}}[lnp(\pmb{x},\pmb{y}|\theta)|\theta^i,\pmb{y}]+H(p(\pmb{x}|\theta^i,\pmb{y})\} \nonumber \\
&=& argmax_{\theta}Q(\theta,\theta^{i}) \quad H(p(\pmb{x}|\theta^i,\pmb{y}) 是定值 \nonumber
\end{eqnarray}

2.2. 利用 F 函数推导 EM 算法性质

      EM 算法两种推导方法的结果都是极大化 (Q(\theta,\theta^{i}) )。F 函数 让 EM 算法的单调性和一致性比较直观。

      直观地推导 EM 算法的单调性。 EM算法得到的参数序列为 (\theta^0,\theta^1,…,\theta^k,..) ,则有 (lnp(y|\theta^i) \le lnp(y|\theta^{i+1}), \forall i)。EM 执行第 i 次迭代的 E 步之后,这是F函数为 (F(\theta^i,\bar{p}{\theta^i}))
\begin{eqnarray}
F(\theta^i,\bar{p}
{\theta^i}) &=& E_{\pmb{x}}[lnp(\pmb{x},\pmb{y}|\theta)|\theta^i,\pmb{y}]+H(p(\pmb{x}|\theta^i,\pmb{y})) \quad \bar{p}{\theta^i} = p(\pmb{x}|\pmb{y},\theta^i) \nonumber \
&=& E
{\pmb{x}}[ln\frac{p(\pmb{x},\pmb{y}|\theta^i)}{p(\pmb{x}|\theta^i,\pmb{y})}|\theta^i,\pmb{y}] \nonumber \
&=& E_{\pmb{x}}[lnp(Y|\theta^{i})|\theta^i,\pmb{y}]= lnp(\pmb{y}|\theta^i)
\label{eq:fproperty2}
\end{eqnarray}
这说明我们每次执行完E步,F函数等于目标函数(lnp(y|\theta^i))。又由于(\theta^0,\theta^1,…,\theta^k,..)是EM算法得到的参数序列,F函数随着i增大而增大,则容易证明(lnp(y|\theta^i))随着i增大而增大。如下图所示,黑线是F函数,红线是(lnp(\pmb{y}|\theta))。F函数是(lnp(\pmb{y}|\theta))的下界函数,每次完成M和E步骤,F都上升,并和(lnp(\pmb{y}|\theta))交汇一次。

em_figure

      直观地推导 EM 算法的一致性。F 函数在(\theta^{})和(\bar{p}^{})达到局部最大值,则 (lnp(y|\theta)) 也在 (\theta^{}) 局部最大值。(lnp(y|\theta) = F(\bar{p}{\theta}, \theta)) 对于任意 (\theta) 成立。特别地,对于使 (F(\bar{p}{\theta}, \theta)) 极大的参数 (\theta^),有
\begin{eqnarray}
lnp(y|\theta^) = F(\bar{p}_{\theta^}, \theta^)
\end{eqnarray}
我们需证明(lnp(y|\theta^
))也是极大值。假设(lnp(y|\theta^))不是极大值,则存在接近(在任意邻域内)(\theta^)的参数(\theta^{}),使得(F(\bar{p}_{\theta^*}, \theta^*) = lnp(y|\theta^*) < lnp(y|\theta^{}) =F(\bar{p}_{\theta^{}}, \theta^{}))。因为(\bar{p_{\theta}})是(\theta)的连续函数,则(\bar{p}_{\theta^{**}})接近(\bar{p}_{\theta^{*}})。这与F函数在(\theta^{*})和(\bar{p}^{*})达到极大值矛盾。

3. EM 算法求解高斯混合模型

      高斯混合模型 (Gauss Mixed Model, GMM) 是概率生成模型,其按多项式分布从 K 个类别中选择一个类别,然后按照该类别对应的高斯分布生成样例。让 (\pmb{\pi}) 表示多项分布的参数,(\pi_{i}) 表示第 i 个项目被选择的概率;让(\pmb{u}) 和 (\pmb{\sigma}) 分别表示高斯分布的期望和标准差, 其中(u_{i}) 和 (\sigma_{i}) 分别表示第 i 个高斯分别的期望和标准差。下图是 k = 2, (\pi_0 = 1/3), (\pi_1 = 2/3), (u_0 = -2), (u_1=0),(\sigma_0 = 0.571), (\sigma_1 = 0.418) 时高斯混合模型的概率密度近似图。

gmm visualization

      
用 EM 算法求解高斯混合模型,可见变量 (\pmb{y}) 是生成的 n 个样例,隐变量 (\pmb{x}) 是被选中的类别 (其中 (x_i) 表示生成第 i 个样例时选中的类别), 参数 (\theta) 包括(\pmb{\pi}),(\pmb{u}) 和 (\pmb{\sigma})。

3.1. E 步推导

      EM 算法的 E 步是计算 (Q(\theta,\theta^i))。让(Z = \sum_{\pmb{x}} \prod_{s=1}^{n} \pi_{x_s}^i N(u_{x_s}^i,\sigma_{x_s}^i)),我们有
\begin{eqnarray}
&&Q(\theta,\theta^i) \nonumber \
&=& \sum_{\pmb{x}} lnp(\pmb{y}|\pmb{x},\theta) p(\pmb{x}|\pmb{y},\theta^{i}) \nonumber \
&=& \sum_{\pmb{x}} {\sum_{t=1}^{n} ln\pi_{x_t}N(y_t;u_{x_t},\sigma_{x_t})}{\frac{\prod_{s=1}^{n} \pi_{x_s}^i N(y_s;u_{x_s}^i,\sigma_{x_s}^i)}{Z}} \nonumber \
&=& \sum_{t=1}^{n} \sum_{\pmb{x}} { ln\pi_{x_t}N(y_t;u_{x_t},\sigma_{x_t})\frac{\prod_{s=1}^{n} \pi_{x_s}^i N(y_s;u_{x_s}^i,\sigma_{x_s}^i)}{Z}} \nonumber \
&=& \sum_{t=1}^{n} \sum_{k=1}^{K} { ln \pi_{k}N(y_t;u_{k},\sigma_{k}) \sum_{\pmb{x}:x_t=k}\frac{\prod_{s=1}^{n} \pi_{x_s}^i N(y_s;u_{x_s}^i,\sigma_{x_s}^i)}{Z}}
\nonumber \
&=& \sum_{t=1}^{n} \sum_{k=1}^{K} { ln \pi_{k}N(y_t;u_{k},\sigma_{k}) \frac{ \pi_{k}^i N(y_t;u_{k}^i,\sigma_{k}^i)}{\sum_{k=1}^{K}\pi_{k}^i N(y_t;u_{k}^i,\sigma_{k}^i)}}\nonumber \
&=& \sum_{t=1}^{n} \sum_{k=1}^{K} { ln \pi_{k}N(y_t;u_{k},\sigma_{k}) w_{t,k}} \nonumber
\end{eqnarray}
E 步计算 (w_{t,k}),便得 (Q(\theta,\theta^i)),代码 (代码地址:GitHub) 如下所示。

def E_step(y, gmm):
    n = len(y);
    K = len(gmm.pi);
    w = [[0.0 for j in xrange(K)] for i in xrange(n)];
    
    for i in xrange(n):
        sum1 = 0.0;
        for j in xrange(K):
            sum1 += gmm.pi[j] * normal(gmm.u[j], gmm.d[j], y[i]);
        for j in xrange(K):
            w[i][j]  = gmm.pi[j] * normal(gmm.u[j],gmm.d[j], y[i]);
            w[i][j] /= sum1;
    return w; 

3.2. M 步推导

      EM 算法的 M 步要极大化 (Q(\theta,\theta^i))。在给定 (w_{t,k})的情况下,我们用拉格朗日方法在 (\sum_{k=1}^{K}\pi_{k} = 1) 约束下极大化 (Q(\theta,\theta^i)) ,求得更新公式。
\begin{eqnarray}
\pi_{k} &=& \frac{N_k}{N} \nonumber \
u_{k} &=& \frac{1}{N_k} \sum_{t=1}^{n} w_{t,k} y_t \nonumber \
\sigma_{k} &=& \frac{1}{N_k} \sum_{t=1}^{n} w_{t,k} (y_t – u_k)^2 \nonumber
\end{eqnarray}

其中 (N = \sum_{k=1}^{K}\sum_{t=1}^{n} w_{t,k}) 和
(N_k = \sum_{t=1}^{n} w_{t,k} )。 M 步的代码 (代码地址:GitHub)。如下所示。

def M_step(y, gmm, w):
    n = len(y);
    K = len(gmm.pi);
    N_k = [0.0 for i in xrange(K)];
    N   = 0.0;
    for j in xrange(K):
        for i in xrange(n):
            N_k[j] += w[i][j];
        N += N_k[j];

    for j in xrange(K):
        gmm.pi[j] = N_k[j] / N;
        gmm.u[j]  = 0.0;
        gmm.d[j]  = 0.0;
        for i in xrange(n):
            gmm.u[j] += w[i][j] * y[i];
            gmm.d[j] += w[i][j] * math.pow(y[i] - gmm.u[j],2)
        gmm.u[j] /= N_k[j];
        gmm.d[j] /= N_k[j];
        gmm.d[j]  = math.sqrt(gmm.d[j]);

    return gmm;

3.2. 一个小实验

      下面是我们做一个小实验,实验代码已经放到 GitHub 上了。下面是实验结果,纵坐标是高斯混合模型的对数似然值的近似值,横坐标是 EM 算法的迭代次数。可以看到,随着 EM 算法的推进,对数似然值在上升。

em gmm result

      

附录

       证明公式(\bar{p}{\theta^i}(\pmb{x})=p(\pmb{x}|\theta^i,\pmb{y}))。对固定的(\theta),极大化F函数是有约束的优化问题:
\begin{eqnarray}
&max& \quad E
{\bar{p}}[lnp(\pmb{x},\pmb{y}|\theta)]+H(\bar{p}(\pmb{x})) \nonumber \
&sub& \quad \sum_{\pmb{x}} \bar{p}(\pmb{x}) = 1 \nonumber \
\end{eqnarray}
为此,引入拉格朗日系数 (\lambda),拉格朗日函数如下
\begin{eqnarray}
L= E_{\bar{p}}[lnp(\pmb{x},\pmb{y}|\theta)] – E_{\bar{p}}[ln\bar{p}(\pmb{x})]+\lambda(1 – \sum_{\pmb{x}} \bar{p}(\pmb{x}) ) \nonumber \
\end{eqnarray}
求拉格朗日函数对(\bar{p})的偏导数:
\begin{eqnarray}
\frac{\partial L}{\partial \bar{p}(\pmb{x})}= lnp(\pmb{x},\pmb{y}|\theta)-ln\bar{p}(\pmb{x})-1-\lambda = 0 \nonumber \
\end{eqnarray}
得(\bar{p}(\pmb{x}) = \frac{p(\pmb{x},\pmb{y}|\theta)}{\exp(1+\lambda)}), 又(\sum_{\pmb{x}} \bar{p}(\pmb{x}) = 1), 所以我们有
\begin{eqnarray}
\bar{p} (\pmb{x})=p(\pmb{x}|\theta,\pmb{y}) \nonumber \
\end{eqnarray}
即公式成立。

      文章结尾欢迎关注我的公众号,每周日的更新就会有提醒哦~

weixin_saomiao

      

此条目发表在算法荟萃分类目录,贴了标签。将固定链接加入收藏夹。

EM算法原理和应用》有 13 条评论

  1. 大话王说:

    公式2的熵的哪项前面为啥负号呢。。

  2. 很不错,你推得和这个是一样的。http://www.cs.cmu.edu/~tom/10-702/Zoubin-702.pdf

  3. larry说:

    倒数第4号的对“求拉格朗日函数对p¯的偏导数”式子中“-1“怎么来的?

  4. O(∩_∩)O说:

    为什么叫做另一种解释呢?

  5. 匿名说:

    文章倒数第4行开始,p¯(x)=p(x,y|θ) / exp(1+λ), 又∑x p¯(x)=1 是怎么证明得出p¯(x)=p(x|θ,y),这里这几天一直没弄懂,lz能帮忙解答一下吗

O(∩_∩)O进行回复 取消回复

电子邮件地址不会被公开。