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

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

1.Jensen不等式

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

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

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

假如我们有一个估计问题（estimation problem），其中由训练样本集 {x(1),...,x(m)}\{x^{(1)}, ..., x^{(m)}\}{x(1),...,x(m)} 包含了 $m$ 个独立样本。我们用模型 $p(x,z)$ 对数据进行建模，拟合其参数（parameters），其中的似然函数（likelihood）如下所示：
$$\begin{array}{rl}l(\theta )& =\sum _{i=1}^m\log p(x;\theta )\\ & =\sum _{i=1}^m\log \sum _{z}p(x,z;\theta )\end{array}$$
对于每个 $i$，设 $Q_i$ 是某个对 $z$ 的分布（${\sum }_z{Q}_i(z)=1,Q_i(z)\ge 0$）。则有下列各式${}^1$：
$$\begin{array}{rlr}\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{array}$$

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

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

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

令
$$g(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$
因为：
$$z^{(i)}\sim Q_i(z^{(i)})$$
所以：
$$g(z^{(i)})\sim Q_i(z^{(i)})$$
于是：
$$\begin{array}{rl}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{array}$$
由Jenson不等式：
$$log(E\left(g(z)\right))\ge E[log(g(z))]$$
不妨再令
$$h(z)=log(g(z))=\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(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)})$$
所以可以得到：
$$\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)})}$$
得证。

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

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

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

  （$E$ 步骤）对每个 $i$，设
$$Q_i(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta )$$
  （$M$ 步骤） 设
$$\theta :=arg\underset{\theta }{\max }\sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$
}

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

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

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

由条件可知：
$$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)})}$$
参数 ${\theta }^{(t+1)}$ 可以通过对上面等式中等号右侧进行最大化而得到。所以有：
$$\begin{array}{rl}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{array}$$
上面的第一个不等式推自：
$$l(\theta )\ge \sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$
上面这个不等式对于任意值的 $Q_i$ 和 $\theta$ 都成立，尤其当 $Q_i=Q_i^{(t)},\theta ={\theta }^{(t+1)}$。由于：
$${\theta }^{(t+1)}=arg\underset{\theta }{\max }\sum _i\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$
所以
$$\begin{array}{rl}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{array}$$

3.高斯混合

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

接下来在 $M$ 步骤，就要最大化关于参数 $\varphi ,\mu ,\Sigma$的值：
$$\begin{array}{rl}\sum _{i=1}^m& \sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\varphi ,\mu ,\Sigma )}{Q_i(z^{(i)})}\\ & =\sum _{i=1}^m\sum _{j=1}^k{Q}_i(z^{(i)}=j)log\frac{p(x^{(i)}|z^{(i)}=j;\mu ,\Sigma )p(z^{(i)}=j;\varphi )}{Q_i(z^{(i)}=j)}\\ & =\sum _{i=1}^m\sum _{j=1}^k{w}_j^{(i)}log\frac{\frac{1}{(2\pi {)}^{n/2}|{\Sigma }_j{|}^{1/2}}exp(-\frac{1}{2}(x^{(i)}-{\mu }_j{)}^T{\Sigma }_j^{-1}(x^{(i)}-{\mu }_j))\cdot {\varphi }_j}{w_j^{(i)}}\end{array}$$
先关于 ${\mu }_l$ 来进行最大化。如果去关于 ${\mu }_l$ 的（偏）导数（derivative），得到：
$$\begin{array}{rl}{\nabla }_{{\mu }_l}& \sum _{i=1}^m\sum _{j=1}^k{w}_j^{(i)}log\frac{\frac{1}{(2\pi {)}^{n/2}|{\Sigma }_j{|}^{1/2}}exp(-\frac{1}{2}(x^{(i)}-{\mu }_j{)}^T{\Sigma }_j^{-1}(x^{(i)}-{\mu }_j))\cdot {\varphi }_j}{w_j^{(i)}}\\ & =-{\nabla }_{{\mu }_l}\sum _{i=1}^m\sum _{j=1}^k{w}_j^{(i)}\frac{1}{2}(x^{(i)}-{\mu }_j{)}^T{\Sigma }_j^{-1}(x^{(i)}-{\mu }_j)\\ & =\frac{1}{2}\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{array}$$
设上式为零，然后解出 ${\mu }_l$ 就产生了更新规则（update rule）：
$${\mu }_l:=\frac{{\sum }_{i=1}^m{w}_l^{(i)}{x}^{(i)}}{{\sum }_{i=1}^m{w}_l^{(i)}}$$
推导在 $M$ 步骤中参数 ${\varphi }_j$ 的更新规则。把仅关于参数 ${\varphi }_j$ 的表达式结合起来，就能发现只需要最大化下面的表达式：
$$\sum _{i=1}^m\sum _{j=1}^k{w}_j^{(i)}log{\varphi }_j$$
然而，还有一个附加的约束，即 ${\varphi }_j$ 的和为$1$，因为其表示的是概率 ${\varphi }_j=p(z^{(i)}=j;\varphi )$。为了保证这个约束条件成立，即 ${\sum }_{j=1}^k{\varphi }_j=1$，我们构建一个拉格朗日函数（Lagrangian）：
$$\mathcal{L}(\varphi )=\sum _{i=1}^m\sum _{j=1}^k{w}_j^{(i)}log{\varphi }_j+\beta (\sum _{j=1}^k{\varphi }_j-1)$$
其中的 $\beta$ 是 拉格朗日乘数（Lagrange multiplier）${}^2$ 。求导，然后得到：
$$\frac{\partial }{\partial {\varphi }_j}\mathcal{L}(\varphi )=\sum _{i=1}^m\frac{w_j^{(i)}}{{\varphi }_j}+1$$

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

设导数为零，然后解方程，就得到了：
$${\varphi }_j=\frac{{\sum }_{i=1}^m{w}_j^{(i)}}{-\beta }$$
也就是说，${\varphi }_j\propto {\sum }_{i=1}^m{w}_j^{(i)}$。结合约束条件（constraint）${\Sigma }_j{\varphi }_j=1$，可以很容易地发现 $-\beta ={\sum }_{i=1}^m{\sum }_{j=1}^k{w}_j^{(i)}={\sum }_{i=1}^{m}1=m$. （这里用到了条件 $w_j^{(i)}=Q_i(z^{(i)}=j)$，而且因为所有概率之和等于$1$，即${\sum }_j{w}_j^{(i)}=1$）。这样我们就得到了在 $M$ 步骤中对参数 ${\varphi }_j$ 进行更新的规则了：
$${\varphi }_j:=\frac{1}{m}\sum _{i=1}^m{w}_j^{(i)}$$
接下来对 $M$ 步骤中对 ${\Sigma }_j$ 的更新规则的推导就很容易了。
$${\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)}}.$$