机器学习（十）——期望值最大算法(EM算法)

10.期望值最大算法(EM算法)

1.Jensen不等式

设 fff 为一个函数，其定义域（domain）为整个实数域（set of real numbers）。这里要回忆一下，如果函数 fff 的二阶导数 f′′(x)≥0f''(x) \ge 0f′′(x)≥0 （其中的 x∈Rx \in Rx∈R），则函数 fff 为一个凸函数（convex function）。如果输入的为向量变量，那么这个函数就泛化了，这时候该函数的海森矩阵（hessian） HHH 就是一个半正定矩阵（positive semi-definite H≥0H \ge 0H≥0）。如果对于所有的 xxx ，都有二阶导数 f′′(x)>0f''(x) > 0f′′(x)>0，那么我们称这个函数 fff 是严格凸函数（对应向量值作为变量的情况，对应的条件就是海森矩阵必须为正定，写作 H>0H > 0H>0）。这样就可以用如下方式来表述 Jensen 不等式：

定理（Theorem）： 设 fff 是一个凸函数，且设 XXX 是一个随机变量（random variable）。然后则有：
E[f(X)]≥f(EX).E[f(X)] \ge f(EX). E[f(X)]≥f(EX).
Jensen 不等式也适用于凹函数（concave）fff，但不等式的方向要反过来，也就是对于凹函数，E[f(X)]≤f(EX)E[f(X)] \le f(EX)E[f(X)]≤f(EX)。

2.期望最大算法(EM算法)

假如我们有一个估计问题（estimation problem），其中由训练样本集 {x(1),...,x(m)}\{x^{(1)}, ..., x^{(m)}\}{x(1),...,x(m)} 包含了 mmm 个独立样本。我们用模型 p(x,z)p(x, z)p(x,z) 对数据进行建模，拟合其参数（parameters），其中的似然函数（likelihood）如下所示：
l(θ)=∑i=1mlog⁡p(x;θ)=∑i=1mlog⁡∑zp(x,z;θ)\begin{aligned} l(\theta) &= \sum_{i=1}^m\log p(x;\theta) \\ &= \sum_{i=1}^m\log\sum_z p(x,z;\theta) \end{aligned} l(θ)=i=1∑mlogp(x;θ)=i=1∑mlogz∑p(x,z;θ)
对于每个 iii，设 QiQ_iQi 是某个对 zzz 的分布（∑zQi(z)=1,Qi(z)≥0\sum_z Q_i(z) = 1, Q_i(z)\ge 0∑zQi(z)=1,Qi(z)≥0）。则有下列各式1^11：
∑ilog⁡p(x(i);θ)=∑ilog⁡∑z(i)p(x(i),z(i);θ)(1)=∑ilog⁡∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))(2)≥∑i∑z(i)Qi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))(3)\begin{aligned} \sum_i\log p(x^{(i)};\theta) &= \sum_i\log\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)&(1) \\ &= \sum_i\log\sum_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} &(2)\\ &\ge \sum_i\sum_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}&(3) \end{aligned} i∑logp(x(i);θ)=i∑logz(i)∑p(x(i),z(i);θ)=i∑logz(i)∑Qi(z(i))Qi(z(i))p(x(i),z(i);θ)≥i∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)(1)(2)(3)

1 如果 zzz 是连续的，那么 QiQ_iQi 就是一个密度函数（density），上面讨论中提到的对 zzz 的求和（summations）就要用对 zzz 的积分（integral）来替代。

上面推导（derivation）的最后一步使用了 Jensen 不等式（Jensen’s inequality）。其中的 f(x)=logxf(x) = log xf(x)=logx 是一个凹函数（concave function），因为其二阶导数 f′′(x)=−1/x2<0f''(x) = -1/x^2 < 0f′′(x)=−1/x2<0 在整个定义域（domain） x∈R+x\in R^+x∈R+ 上都成立。

由(2)式到(3)式证明：

令
g(z(i))=p(x(i),z(i);θ)Qi(z(i))g(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} g(z(i))=Qi(z(i))p(x(i),z(i);θ)
因为：
z(i)∼Qi(z(i))z^{(i)}\sim Q_i(z^{(i)}) z(i)∼Qi(z(i))
所以：
g(z(i))∼Qi(z(i))g(z^{(i)})\sim Q_i(z^{(i)}) g(z(i))∼Qi(z(i))
于是：
E(g(z))=∑z(i)g(z(i))P(g(z(i)))=∑z(i)p(x(i),z(i);θ)Qi(z(i))Qi(z(i))\begin{aligned} E\left(g(z)\right)&=\sum_{z^{(i)}}g(z^{(i)})P\left(g(z^{(i)})\right)\\ &=\sum_{z^{(i)}}\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}Q_i(z^{(i)}) \end{aligned} E(g(z))=z(i)∑g(z(i))P(g(z(i)))=z(i)∑Qi(z(i))p(x(i),z(i);θ)Qi(z(i))
由Jenson不等式：
log(E(g(z)))≥E[log(g(z))]log(E\left(g(z)\right))\ge E[log(g(z))] log(E(g(z)))≥E[log(g(z))]
不妨再令
h(z)=log(g(z))=log⁡p(x(i),z(i);θ)Qi(z(i))h(z)=log(g(z))=\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} h(z)=log(g(z))=logQi(z(i))p(x(i),z(i);θ)
所以同理可以得到
E[log(g(z))]=E(h(z))=∑z(i)h(z(i))P(h(z(i)))=∑z(i)h(z(i))Qi(z(i))E[log(g(z))]=E(h(z))=\sum_{z^{(i)}}h(z^{(i)})P\left(h(z^{(i)})\right)=\sum_{z^{(i)}}h(z^{(i)})Q_i(z^{(i)}) E[log(g(z))]=E(h(z))=z(i)∑h(z(i))P(h(z(i)))=z(i)∑h(z(i))Qi(z(i))
所以可以得到：
log⁡∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))≥∑z(i)Qi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))\log\sum_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \ge \sum_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} logz(i)∑Qi(z(i))Qi(z(i))p(x(i),z(i);θ)≥z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)
得证。

为了使Jenson不等式求得等号，即使得：
f(E(x))=E(f(x))f(E(x)) = E(f(x)) f(E(x))=E(f(x))
则必须变量xxx为一个常量，则f(x)f(x)f(x)也将变为一个常量，所以有：
E(x)=xE(f(x))=f(x)E(x)=x\\ E(f(x))=f(x) E(x)=xE(f(x))=f(x)
所以可以推出：
f(E(x))=f(x)=E(f(x))f(E(x))=f(x)=E(f(x)) f(E(x))=f(x)=E(f(x))
也就是需要：
p(x(i),z(i);θ)Qi(z(i))=c\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c Qi(z(i))p(x(i),z(i);θ)=c
其中常数 ccc 不依赖 z(i)z^{(i)}z(i)。要实现这一条件，只需满足：
Qi(z(i))∝p(x(i),z(i);θ)Q_i(z^{(i)})\propto p(x^{(i)},z^{(i)};\theta) Qi(z(i))∝p(x(i),z(i);θ)
所以：
Qi(z(i))=p(x(i),z(i);θ)cQ_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{c} Qi(z(i))=cp(x(i),z(i);θ)
又因为z(i)z^{(i)}z(i)是一个分布,所以：
∑zQi(z(i))=1=∑zp(x(i),z(i);θ)c\begin{aligned} \sum_z Q_i(z^{(i)}) &= 1\\ &=\frac{\sum_z p(x^{(i)},z^{(i)};\theta)}{c} \end{aligned} z∑Qi(z(i))=1=c∑zp(x(i),z(i);θ)
因此我们可以得出：
c=∑zp(x(i),z(i);θ)=p(x(i);θ)c=\sum_z p(x^{(i)},z^{(i)};\theta)=p(x^{(i)};\theta) c=z∑p(x(i),z(i);θ)=p(x(i);θ)
所以：
Qi(z(i))=p(x(i),z(i);θ)p(x(i);θ)=p(z(i)∣x(i);θ)\begin{aligned} Q_i(z^{(i)}) &= \frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \\ &= p(z^{(i)}|x^{(i)};\theta) \end{aligned} Qi(z(i))=p(x(i);θ)p(x(i),z(i);θ)=p(z(i)∣x(i);θ)
因此，在给定 x(i)x^{(i)}x(i) 和参数 θ\thetaθ 的设置下，我们可以简单地把 QiQ_iQi 设置为 z(i)z^{(i)}z(i) 的后验分布（posterior distribution）。

接下来，对 QiQ_iQi 的选择，等式 (3)(3)(3) 就给出了似然函数对数（log likelihood）的一个下限，而似然函数（likelihood）正是我们要试图求最大值（maximize）的。这就是 EEE 步骤。接下来在算法的 MMM 步骤中，就最大化等式 (3)(3)(3) 当中的方程，然后得到新的参数 θ\thetaθ。重复这两个步骤，就是完整的 EMEMEM 算法，如下所示：

重复下列过程直到收敛（convergence）: {

（EEE 步骤）对每个 iii，设
Qi(z(i)):=p(z(i)∣x(i);θ)Q_i(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta) Qi(z(i)):=p(z(i)∣x(i);θ)
（MMM 步骤）设
θ:=argmax⁡θ∑i∑z(i)Qi(z(i))logp(x(i),z(i);θ)Qi(z(i))\theta := arg\max_\theta\sum_i\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} θ:=argθmaxi∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)
}

下面我们要来证明这个算法是递增的(收敛的)：

设 θ(t)\theta^{(t)}θ(t) 和 θ(t+1)\theta^{(t+1)}θ(t+1) 是上面 EMEMEM 迭代过程中的某两个参数（parameters）

证明： l(θ(t))≤l(θ(t+1))l(\theta^{(t)})\le l(\theta^{(t+1)})l(θ(t))≤l(θ(t+1))

由条件可知：
l(θ(t))=∑i∑z(i)Qi(t)(z(i))logp(x(i),z(i);θ(t))Qi(t)(z(i))l(\theta^{(t)})=\sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} l(θ(t))=i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t))
参数 θ(t+1)\theta^{(t+1)}θ(t+1) 可以通过对上面等式中等号右侧进行最大化而得到。所以有：
l(θ(t+1))≥∑i∑z(i)Qi(t)(z(i))logp(x(i),z(i);θ(t+1))Qi(t)(z(i))\begin{aligned} l(\theta^{(t+1)}) &\ge \sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \end{aligned} l(θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t+1))
上面的第一个不等式推自：
l(θ)≥∑i∑z(i)Qi(z(i))logp(x(i),z(i);θ)Qi(z(i))l(\theta)\ge \sum_i\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} l(θ)≥i∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)
上面这个不等式对于任意值的 QiQ_iQi 和 θ\thetaθ 都成立，尤其当 Qi=Qi(t),θ=θ(t+1)Q_i = Q_i^{(t)}, \theta = \theta^{(t+1)}Qi=Qi(t),θ=θ(t+1)。由于：
θ(t+1)=argmax⁡θ∑i∑z(i)Qi(z(i))logp(x(i),z(i);θ)Qi(z(i))\theta^{(t+1)} =arg\max_\theta \sum_i\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} θ(t+1)=argθmaxi∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)
所以
l(θ(t+1))≥∑i∑z(i)Qi(t)(z(i))logp(x(i),z(i);θ(t+1))Qi(t)(z(i))≥∑i∑z(i)Qi(t)(z(i))logp(x(i),z(i);θ(t))Qi(t)(z(i))=l(θ(t))\begin{aligned} l(\theta^{(t+1)}) &\ge \sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \\ &\ge \sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \\ &= l(\theta^{(t)}) \end{aligned} l(θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t))=l(θ(t))

3.高斯混合

EEE 步骤很简单。还是按照上面的算法推导过程，只需要计算：
wj(i)=Qi(z(i)=j)=P(z(i)=j∣x(i);ϕ,μ,Σ)w_j^{(i)}=Q_i(z^{(i)}=j)=P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma) wj(i)=Qi(z(i)=j)=P(z(i)=j∣x(i);ϕ,μ,Σ)
这里面的 “Qi(z(i)=j)”“Q_i(z^{(i)} = j)”“Qi(z(i)=j)” 表示的是在分布 QiQ_iQi上 z(i)z^{(i)}z(i) 取值 jjj 的概率。

接下来在 MMM 步骤，就要最大化关于参数 ϕ,μ,Σ\phi,\mu,\Sigmaϕ,μ,Σ的值：
∑i=1m∑z(i)Qi(z(i))logp(x(i),z(i);ϕ,μ,Σ)Qi(z(i))=∑i=1m∑j=1kQi(z(i)=j)logp(x(i)∣z(i)=j;μ,Σ)p(z(i)=j;ϕ)Qi(z(i)=j)=∑i=1m∑j=1kwj(i)log1(2π)n/2∣Σj∣1/2exp(−12(x(i)−μj)TΣj−1(x(i)−μj))⋅ϕjwj(i)\begin{aligned} \sum_{i=1}^m&\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_i(z^{(i)})}\\ &= \sum_{i=1}^m\sum_{j=1}^kQ_i(z^{(i)}=j)log\frac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_i(z^{(i)}=j)} \\ &= \sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}log\frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp(-\frac 12(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\cdot\phi_j}{w_j^{(i)}} \end{aligned} i=1∑mz(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);ϕ,μ,Σ)=i=1∑mj=1∑kQi(z(i)=j)logQi(z(i)=j)p(x(i)∣z(i)=j;μ,Σ)p(z(i)=j;ϕ)=i=1∑mj=1∑kwj(i)logwj(i)(2π)n/2∣Σj∣1/21exp(−21(x(i)−μj)TΣj−1(x(i)−μj))⋅ϕj
先关于 μl\mu_lμl 来进行最大化。如果去关于 μl\mu_lμl 的（偏）导数（derivative），得到：
∇μl∑i=1m∑j=1kwj(i)log1(2π)n/2∣Σj∣1/2exp(−12(x(i)−μj)TΣj−1(x(i)−μj))⋅ϕjwj(i)=−∇μl∑i=1m∑j=1kwj(i)12(x(i)−μj)TΣj−1(x(i)−μj)=12∑i=1mwl(i)∇μl(2μlTΣl−1x(i)−μlTΣl−1μl)=∑i=1mwl(i)(Σl−1x(i)−Σl−1μl)\begin{aligned} \nabla_{\mu_l}&\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}log\frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp(-\frac 12(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j))\cdot\phi_j}{w_j^{(i)}} \\ &= -\nabla_{\mu_l}\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}\frac 12(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j) \\ &= \frac 12\sum_{i=1}^m w_l^{(i)}\nabla_{\mu_l}(2\mu_l^T\Sigma_l^{-1}x^{(i)}-\mu_l^T\Sigma_l^{-1}\mu_l) \\ &= \sum_{i=1}^m w_l^{(i)}(\Sigma_l^{-1}x^{(i)}-\Sigma_l^{-1}\mu_l) \end{aligned} ∇μli=1∑mj=1∑kwj(i)logwj(i)(2π)n/2∣Σj∣1/21exp(−21(x(i)−μj)TΣj−1(x(i)−μj))⋅ϕj=−∇μli=1∑mj=1∑kwj(i)21(x(i)−μj)TΣj−1(x(i)−μj)=21i=1∑mwl(i)∇μl(2μlTΣl−1x(i)−μlTΣl−1μl)=i=1∑mwl(i)(Σl−1x(i)−Σl−1μl)
设上式为零，然后解出 μl\mu_lμl 就产生了更新规则（update rule）：
μl:=∑i=1mwl(i)x(i)∑i=1mwl(i)\mu_l := \frac{\sum_{i=1}^m w_l^{(i)}x^{(i)}}{\sum_{i=1}^m w_l^{(i)}} μl:=∑i=1mwl(i)∑i=1mwl(i)x(i)
推导在 MMM 步骤中参数 ϕj\phi_jϕj 的更新规则。把仅关于参数 ϕj\phi_jϕj 的表达式结合起来，就能发现只需要最大化下面的表达式：
∑i=1m∑j=1kwj(i)logϕj\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}log\phi_j i=1∑mj=1∑kwj(i)logϕj
然而，还有一个附加的约束，即 ϕj\phi_jϕj 的和为111，因为其表示的是概率 ϕj=p(z(i)=j;ϕ)\phi_j = p(z^{(i)} = j;\phi)ϕj=p(z(i)=j;ϕ)。为了保证这个约束条件成立，即 ∑j=1kϕj=1\sum^k_{j=1}\phi_j = 1∑j=1kϕj=1，我们构建一个拉格朗日函数（Lagrangian）：
L(ϕ)=∑i=1m∑j=1kwj(i)logϕj+β(∑j=1kϕj−1)\mathcal L(\phi)=\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}log\phi_j+\beta(\sum^k_{j=1}\phi_j - 1) L(ϕ)=i=1∑mj=1∑kwj(i)logϕj+β(j=1∑kϕj−1)
其中的 β\betaβ 是拉格朗日乘数（Lagrange multiplier）2^22 。求导，然后得到：
∂∂ϕjL(ϕ)=∑i=1mwj(i)ϕj+1\frac{\partial}{\partial{\phi_j}}\mathcal L(\phi)=\sum_{i=1}^m\frac{w_j^{(i)}}{\phi_j}+1 ∂ϕj∂L(ϕ)=i=1∑mϕjwj(i)+1

2 这里我们不用在意约束条件 ϕj≥0\phi_j \ge 0ϕj≥0，因为很快就能发现，这里推导得到的解会自然满足这个条件的。

设导数为零，然后解方程，就得到了：
ϕj=∑i=1mwj(i)−β\phi_j=\frac{\sum_{i=1}^m w_j^{(i)}}{-\beta} ϕj=−β∑i=1mwj(i)
也就是说，ϕj∝∑i=1mwj(i)\phi_j\propto \sum_{i=1}^m w_j^{(i)}ϕj∝∑i=1mwj(i)。结合约束条件（constraint）Σjϕj=1\Sigma_j \phi_j = 1Σjϕj=1，可以很容易地发现 −β=∑i=1m∑j=1kwj(i)=∑i=1m1=m-\beta = \sum_{i=1}^m\sum_{j=1}^kw_j^{(i)} = \sum_{i=1}^m 1 =m−β=∑i=1m∑j=1kwj(i)=∑i=1m1=m. （这里用到了条件 wj(i)=Qi(z(i)=j)w_j^{(i)} =Q_i(z^{(i)} = j)wj(i)=Qi(z(i)=j)，而且因为所有概率之和等于111，即∑jwj(i)=1\sum_j w_j^{(i)}=1∑jwj(i)=1）。这样我们就得到了在 MMM 步骤中对参数 ϕj\phi_jϕj 进行更新的规则了：
ϕj:=1m∑i=1mwj(i)\phi_j := \frac 1m \sum_{i=1}^m w_j^{(i)} ϕj:=m1i=1∑mwj(i)
接下来对 MMM 步骤中对 Σj\Sigma_jΣj 的更新规则的推导就很容易了。
Σj=∑i=1mwj(i)(x(i)−μj)(x(i)−μj)T∑i=1mwj(i).\Sigma_j=\frac{\sum_{i=1}^m w_j^{(i)}(x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum_{i=1}^m w_j^{(i)}}. Σj=∑i=1mwj(i)∑i=1mwj(i)(x(i)−μj)(x(i)−μj)T.