EM算法的深入理解

在现实应用中，我们往往会遇到“不完整”的训练样本，如某些样本的属性值未知(未观测)，而这些未知的属性值统称为隐变量。由于EM算法在这些未观测属性的情形下仍能对模型参数进行估计，所以其应用及其广泛，比如高斯混合模型，隐马尔可夫模型的非监督学习。本文将从不同角度理解并解释EM算法，最后通过手写数字图像的非监督聚类对EM算法进行进一步解释。

1. EM算法的理解

假设$X$为观测变量，$Z$表示未观测变量，$\theta$为模型参数。在概率模型中，我们常常根据最大似然估计方法进行参数的估计，即极大化观测数据关于参数的对数似然函数：

$$log \,P(X|\theta)=log \, \sum_Z P(X,Z|\theta)$$
注意到上式包含隐变量以及和的对数，而无法有效求解。而EM算法这一利器正是避开了上式的求解，而是通过Jensen不等式( $log\sum_j \lambda_jy_j \geq \sum_j \lambda_j log y_j$，其中 $\lambda_j \geq 0$，且 $\sum_j \lambda_j=1$)找到其下界(lower bound)，通过不断求解下界的极大化来逼近求解对数似然函数极大化。为了使用Jensen不等式，我们必须把对数似然函数构造成类似的式子如下：
$$log \, \sum_Z P(X,Z|\theta)=log \, \sum_Z q(Z)\frac{P(X,Z|\theta)}{q(Z)}$$
其中上式中 $q(Z)$为概率分布，很显然满足Jensen不等式的条件，因此我们进一步得到
$$\begin{array}{c} log \, \sum_Z P(X,Z|\theta) &= log \, \sum_Z q(Z)\frac{P(X,Z|\theta)}{q(Z)} \\ &\geq \sum_Z q(Z) log \frac{P(X,Z|\theta)}{q(Z)} \end{array}$$
那么对数似然函数的下界是 $L(q,\theta)=\sum_Z q(Z) log \frac{P(X,Z|\theta)}{q(Z)}$。可以看出下界是关于 $q(Z)$和 $\theta$的函数，因此EM算法分为两步：

• 固定$\theta$，得到$q(Z)$的分布。如果进一步求隐变量$Z$的期望，则对应EM算法的E-step。这里值得注意的是，$q(Z)$的分布该如何确定。现在我们回顾一下在Jensen不等式当中，只有当$y$为常数时，等式成立，也即当$\frac{P(X,Z|\theta)}{q(Z)}=C$。通过这个条件，我们能轻易得到$q(Z)=p(Z|X,\theta)$。

• 固定$q(Z)$，优化$\theta$，对应EM算法的M-step。这两个步骤不断重复，直至收敛到局部最优解。

一般来讲，EM算法主要用于含有隐变量的概率模型的学习，针对不完全数据$X$的最大对数似然函数找到局部最优解。下面从另一个角度解释EM算法：在现实中，我们很难得到完全数据$\{X,Z\}$，但是我们能根据已有知识能得到隐变量$Z$的后验概率分布$P(Z|X,\theta)$。虽然我们不能建立完全数据的对数似然函数，但是可以考虑隐变量在其后验分布下的期望值，这样有了隐变量的期望值，则可以构建完全数据的对数似然函数。进一步，参数的优化目标为

$$\begin{array}{c} \mathop{max}\limits_{\theta}\mathop{E}\nolimits_{p(Z|X,\theta^{old})} [log \, p(X,Z|\theta)] \\ =\mathop{max}\limits_{\theta}\sum_Z p(Z|X,\theta^{old})\, log \, p(X,Z|\theta) \end{array}$$
我们首先利用当前参数 $\theta^{old}$来得到隐变量的后验分布，在根据完全数据的对数似然函数在后验分布下的期望下更新 $\theta$。

如果我们对对数似然函数的Jensen不等式进行进一步的分析，我们会发现如下等式成立：

$$log p(X|\theta)=L(q,\theta)+KL(q||p)$$
其中
$$\begin{array}{c} L(q,\theta) = \sum_Z q(Z) log \frac{p(X,Z|\theta)}{q(Z)} \\ KL(q||p)=-\sum_Z q(Z) log \frac{p(Z|X,\theta)}{q(Z)} \end{array}$$
$L(q,\theta)$为下界，包含 $X$与$Z$的联合概率分布；而 $KL(q||p)$包含 $Z$的条件分布。Kullback-Leibler divergene 具有非负性，$KL(q||p)\geq 0$，且只有当 $q(Z)=p(Z|X,\theta)$时，等号成立。则此时下界 $L(q,\theta)$等于不完全数据的对数似然函数。

2. 优化最大后验概率

EM算法不仅能优化求解不完全数据的对数似然函数，还能优化求解在不完全数据下参数的后验概率。一般，根据参数的最大后验概率估计属于频率学派，往往需要事先给定参数的先验分布$p(\theta)$。待优化的目标函数为最大化后验概率的对数：

$$log p(\theta|X)=log p(X|\theta)+log p(\theta)-logp(X)$$
根据前面的Jensen不等式，我们能得到 $log p(X|\theta)$的下界 $L(q,\theta)$，则
$$log p(\theta|X)\geq L(q,\theta)+log p(\theta)-log p(X)$$
其中 $log p(X)$为常数项。上式的优化又可交替的优化 $q$和$\theta$。与标准的最大似然对比，紧紧只增加了参数的先验项。

3. 手写数字图像的建模

下面通过一个手写数字图像的非监督分类来加深对EM算法的理解。假设手写数字图像中笔画所占像素点取值为1，而其它像素点取值为0。那么，整个图像向量(将像素点矩阵取值拉成一列向量)$\, x \,$为一个服从伯努力分布的多元随机变量。假设给定无标签数据集$X=\{x_n\}_{n=1}^N$，其中$x_n \in R^D$，现在我们的任务是对其进行分类。

由于给定的数据集中有0~9等10个类别，因此我们引入混合的伯努力分布对数据集进行模型建立，这里我们取$K=10$个伯努力分布。那么观测数据集$X$的对数似然函数为

$$log \, p(X|u,\pi)=\sum_{n=1}^Nlog \, \{\sum_{k=1}^K\pi_k p(x_n|u_k)\}$$
同样，由于和对数的存在，使得最大似然解没有闭合表达式。因此我们采用EM算法进行求解。针对每一个数据$x_n$，我们引入显示的隐变量$z_n$，且$z_n \in R^K$，其中$z_n$中尤且只有一个元素对应为1，其它为0，即表示对应数据$x_n$所属类别。根据EM算法，我们构造完全数据$\{X,Z\}$的对数形式：

$$\begin{array}{c} log \, p(X,Z|u,\pi)=\sum_n^N log \, p(x_n,z_n|u,\pi) \\ =\sum_n^N log \, p(x_n|z_n,u) \, p(z_n|\pi) \end{array}$$
其中
$$\begin{array}{c} p(x_n|z_n,u)=\prod_k p(x_n|u_k)^{z_{nk}} \\ p(x_n|u_k)=\prod_i^D u_{ki}^{x_{ni}}(1-u_{ki})^{(1-{x_{ni}}) }\\ p(z_n|\pi)=\prod_k \pi_k^{z_{nk}} \end{array}$$
其中，值得注意的是 $x_n$中属性值之间在类别确定的情况下假设相互独立，也即没有考虑像素点之间的邻域相关性。接下来，我们构造完全数据的对数形式的期望
$$\begin{array}{c} E_Z[log \, p(X,Z|u,\pi)]=\sum_n^N\sum_k^KE[z_{nk}] \left\{log \,\pi_k \\ +\sum_i^D[x_{ni}log\, u_{ki}+(1-x_{ni})log\, (1-u_{ki})]\right\} \end{array}$$
最后我们通过EM算法的步骤即可对参数 $\{u_k,\pi_k\}_{k=1}^K$的估计。这里值得指出的是 $E[z_{nk}]$的计算，由于 $z_{nk}$为一个二值(0/1)变量，所以 $E[z_{nk}]=p(z_{nk}=1|x_n,u,\pi)$，那么
$$\begin{array}{c} p(z_{nk}=1|x_n,u,\pi)=\frac{p(x_n|z_{nk}=1,u)\, p(z_{nk}=1|\pi)}{\sum_{z_{nk}} p(x_n|z_{nk}=1,u)\, p(z_{nk}=1|\pi)} \end{array}$$
注意当 $z_{nk}=0$时，表示对应数据 $x_n$不属于 $k$类，则属于其他类别。因此上式转化为
$$\begin{array}{c} p(z_{nk}=1|x_n,u,\pi)=\frac{p(x_n|u_k)\, \pi_k}{\sum_k p(x_n|u_k)\, \pi_k} \end{array}$$
这样有了 $E[z_{nk}]$，我们就得到了完全数据的对数似然函数，只需要对参数 $u_k,\pi_k$分别求偏导数，令其等于0，就能得到其闭合解形式。需要注意的是 $\sum_k\pi_k=1$约束。有了模型的参数，我们就能通过数据对应的隐变量 $z_n$得到其后验概率，进而判别数据 $x_n$的类别。

EM算法使用两个步骤交替计算：第一步期望步，利用当前估计参数来计算对数似然的期望值；第二部最大化步，寻找能使第一步产生的似然期望最大化的参数值。通俗的讲，EM算法首先对缺失数据进行估计，再对整个完全数据进行参数估计。