漫谈 Clustering (番外篇): Expectation Maximization

本文是“漫谈 Clustering 系列”中的第 5 篇，参见本系列的其他文章。

Expectation Maximization (EM) 是一种以迭代的方式来解决一类特殊最大似然 (Maximum Likelihood) 问题的方法，这类问题通常是无法直接求得最优解，但是如果引入隐含变量，在已知隐含变量的值的情况下，就可以转化为简单的情况，直接求得最大似然解。

我们会看到，上一次说到的 Gaussian Mixture Model 的迭代求解方法可以算是 EM 算法最典型的应用，而最开始说的 K-means 其实也可以看作是 Gaussian Mixture Model 的一个变种（固定所有的 $\Sigma_k = \epsilon\mathbf{I}$ ，并令 $\epsilon \rightarrow 0$ 即可）。然而 EM 实际上是一种通用的算法（或者说是框架），可以用来解决很多类似的问题，我们最后将以一个中文分词的例子来说明这一点。

为了避免问题变得太抽象，我们还是先从上一次的 Gaussian Mixture Model 说起吧。回顾一下我们之前要解决的问题：求以下 Log-likelihood function 的最大值：

\[ \displaystyle \sum_{i=1}^N \log \left\{\sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)\right\} \]

但是由于在 $\log$ 函数里面又有加和，没法直接求。考虑一下 GMM 生成 sample 的过程：先选一个 Component ，然后再从这个 Component 所对应的那个普通的 Gaussian Distribution 里进行 sample 。我们可以很自然地引入一个隐含变量 $z$ ，这是一个 $K$ 维向量，如果第 $k$ 个 Component 被选中了，那么我们讲其第 $k$ 个元素置为 1 ，其余的全为 0 。那么，再来看看，如果除了之前的 sample $x_i$ 的值之外，我们同时还知道了每个 $x_i$ 所对应的隐含变量 $z_i$ 的值，情况又会变成怎么样呢？

因为我们同时观察到了 $x$ 和 $z$ ，所以我们现在要最大化的似然函数也变为 $\prod_{i=1}^N p(x_i, z_i)$ 。注意到 $p(x, z)$ 可以表示为：

\[ \displaystyle \begin{aligned} p(x,z) & = p(z)p(x|z) \\ & = p(z)\prod_{k=1}^K \mathcal{N}(x|\mu_k, \Sigma_k)^{z_k} \end{aligned} \]

而 $z$ 的概率，当 $z$ 的第 $k$ 个元素为 1 的时候，亦即第 $k$ 个 Component 被选中的时候，这个概率为 $\pi_k$ ，统一地写出来就是：

\[ \displaystyle p(z) = \prod_{k=1}^K \pi_k^{z_k} \]

带入上面个式子，我们得到 $p(x, z)$ 的概率是一个大的乘积式（对比之前 $p(x)$ 是一个和式）。再替换到最开始的那个 Log-likelihood function 中，得到新的同时关于 sample $x$ 和隐含变量 $z$ 的 Log-likelihood:

\[ \displaystyle \sum_{i=1}^N\sum_{k=1}^K z_{ik}\left\{\log\pi_k + \log\mathcal{N}(x_i|\mu_k, \Sigma_k)\right\} \]

情况瞬间逆转，现在 $\log$ 和求和符号换了个位置，直接作用在普通的高斯分布上了，一下子就变成了可以直接求解的问题。不过，事情之所以能发展得如此顺利，完全依赖于一个我们伪造的假设：隐含变量的值已知。然而实际上我们并不知道这个值。问题的结症在这里了，如果我们有了这个值，所有的问题都迎刃而解了。回想一下，在类似的地方，我们是如何处理这样的情况的呢？一个很类似的地方就是（比如，在数据挖掘中）处理缺失数据的情况，一般有几种办法：

• 用取值范围类的随机值代替。
• 用平均值代替。
• 填 0 或者其他特殊值。

这里我们采取一种类似于平均值的办法：取期望。因为这里我们至少有 sample $x$ 的值，因此我们可以把这个信息利用起来，针对后验概率 $p(z|x)$ 来取期望。前面说过，$z$ 的每一个元素只有 0 和 1 两种取值，因此按照期望的公式写出来就是：

\[ \displaystyle \begin{aligned} \mathbb{E}[z_k] &= 0\times p(z_k=0|x) + 1\times p(z_k=1|x) \\ & = p(z_k=1|x) \\ & = \frac{p(z_k=1)p(x|z_k=1)}{p(x)} \\ & = \frac{\pi_k\mathcal{N}(x|\mu_k,\Sigma_k)}{\sum_{j=1}^K\pi_j\mathcal{N}(x|\mu_j, \Sigma_j)} \end{aligned} \]

中间用贝叶斯定理变了一下形，最后得到的式子正是我们在漫谈 GMM 中定义的 $\gamma(i, k)$ 。因此，对于上面那个可以直接求解的 Log-likelihood function 中的 $z_{ik}$ ，我们只要用其期望 $\mathbb{E}(z_ik)$ 亦即 $\gamma(i, k)$ 代替即可。

到这里为止，看起来是一个非常完美的方法，不过仔细一看，发现还有一个漏洞：之前的 Log-likelihood function 可以直接求最大值是建立在 $z_{ik}$ 是已知情况下，现在虽然我们用 $\gamma(i, k)$ 来代替了 $z_{ik}$，但是实际上 $\gamma(i, k)$ 是一个反过来以非常复杂的关系依赖所要求参数的一个式子，而不是一个“已知的数值”。解决这个问题的办法就是迭代。

到此为止，我们就勾勒出了整个 EM 算法的结构，同时也把上一次所讲的求解 GMM 的方法又推导了一遍。EM 名字也来源于此：

• Expectation 一步对应于求关于后验概率的期望亦即 $\gamma(i, k)$；
• Maximization 一步则对应于接下来的正常的最大似然的方法估计相关的参数亦即 $\pi_k$、$\mu_k$ 和 $\Sigma_k$ 。

如果你还没有厌烦这些公式的话，让我们不妨再稍微花一点时间，从偏理论一点的角度来简略地证明一下 EM 这个迭代的过程每一步都会对结果有所改进，除非已经达到了一个（至少是局部的）最优解。

现在我们的讨论将不局限于 GMM ，并使用一些稍微紧凑一点的符号。用 $\mathbf{X}$ 表示所有的 sample ，同时用 $\mathbf{Z}$ 表示所有对应的隐含变量。我们的问题是要通过最大似然的方法估计出 $p(\mathbf{X}|\theta)$ 中的参数 $\theta$ 。在这里我们假设这个问题很困难，但是要优化 $p(\mathbf{X}, \mathbf{Z}|\theta)$ 却很容易。这就是 EM 算法能够解决的那一类问题。

现在我们引入一个关于隐含变量的分布 $q(\mathbf{Z})$。注意到对于任意的 $q(\mathbf{Z})$ ，我们都可以对 Log-likelihood Function 作如下分解：

\[ \displaystyle \begin{aligned} \log p(\mathbf{X}|\theta) &= \sum_{\mathbf{Z}}q(\mathbf{Z})\log\left\{\frac{p(\mathbf{X}, \mathbf{Z}|\theta)}{q(\mathbf{Z})}\right\} + \left( -\sum_{\mathbf{Z}}q(\mathbf{Z})\log\left\{\frac{p(\mathbf{Z}|\mathbf{X}, \theta)}{q(\mathbf{Z})}\right\}\right) \\ &\triangleq \mathcal{L}(q, \theta) + \text{KL}(q\|p) \end{aligned} \]

其中 $\text{KL}(q\|p)$ 是分布 $q(\mathbf{Z})$ 和 $p(\mathbf{Z}|\mathbf{X}, \theta)$ 之间的 Kullback-Leibler divergence 。由于 Kullback-Leibler divergence 是非负的，并且只有当两个分布完全相同的时候才会取到 0 。因此我们可以得到关系 $\mathcal{L}(q, \theta) \leqslant \log p(\mathbf{X}|\theta)$ ，亦即 $\mathcal{L}(q, \theta)$ 是 $\log p(\mathbf{X}|\theta)$ 的一个下界。

现在考虑 EM 的迭代过程，记上一次迭代得出的参数为 $\theta^{\text{old}}$，现在我们要选取 $q(\mathbf{Z})$ 以令 $\mathcal{L}(q, \theta^{\text{old}})$ 最大，由于 $\log p(\mathbf{X}|\theta^{\text{old}})$ 并不依赖于 $q(\mathbf{Z})$ ，因此 $\mathcal{L}(q, \theta^{\text{old}})$ 的上限（在 $\theta^{\text{old}}$ 固定的时候）是一个定值，它取到这个最大值的条件就是 Kullback-Leibler divergence 为零，亦即 $q(\mathbf{Z})$ 等于后验概率 $p(\mathbf{Z}|\mathbf{X}, \theta^{\text{old}})$ 。把它带入到 $\mathcal{L}(q, \theta^{\text{old}})$ 的表达式中可以得到

\[ \displaystyle \begin{aligned} \mathcal{L}(q,\theta) &= \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\theta^{\text{old}}) \log p(\mathbf{X}, \mathbf{Z}|\theta) - \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X}, \theta^{\text{old}}) \log p(\mathbf{Z}|\mathbf{X}, \theta^{\text{old}}) \\ &= \mathcal{Q}(\theta, \theta^{\text{old}}) + \text{const} \end{aligned} \]

其中 const 是常量，而 $\mathcal{Q}(\theta, \theta^{\text{old}})$ 则正是我们之前所得到的“同时包含了 sample 和隐含变量的 Log-likelihood function 关于后验概率的期望”，因此这个对应到 EM 中的“E”一步。

在接下来的“M”一步中，我们讲固定住分布 $q(\mathbf{Z})$ ，再选取合适的 $\theta$ 以将 $\mathcal{L}(q, \theta)$ 最大化，这次其上界 $\log p(\mathbf{X}|\theta)$ 也依赖于变量 $\theta$ ，并会随着 $\mathcal{L}(q, \theta)$ 的增大而增大（因为我们有前面的不等式成立）。一般情况下 $\log p(\mathbf{X}|\theta)$ 增大的量会比 $\mathcal{L}(q, \theta)$ 要多一些，这时候 Kullback-Leibler divergence 在新的参数 $\theta$ 下又不为零了，因此我们可以进入下一轮迭代，重新回到“E”一步去求新的 $q(\mathbf{Z})$ ；另一方面，如果这里 Kullback-Leibler divergence 在新的参数下还是等于 0 ，那么说明我们已经达到了一个（至少是局部的）最优解，迭代过程可以结束了。

上面的推导中我们可以看到每一次迭代 E 和 M 两个步骤都是在对解进行改进，因此迭代的过程中得到的 likelihood 会逐渐逼近（至少是局部的）最优值。另外，在 M 一步中除了用最大似然之外，也可以引入先验使用 Maximum a Posteriori (MAP) 的方法来做。还有一些很困难的问题，甚至在迭代的过程中 M 一步也不能直接求出最大值，这里我们通过把要求放宽──并不一定要求出最大值，只要能够得到比旧的值更好的结果即可，这样的做法通常称作 Generalized EM (GEM)。

当然，一开始我们就说了，EM 是一个通用的算法，并不只是用来求解 GMM 。下面我们就来举一个用 EM 来做中文分词的例子。这个例子来源于论文 Self-supervised Chinese Word Segmentation 。我在上次 MSTC 搜索引擎系列小课堂讲文本处理的时候提到过。这里为了把注意力集中到 EM 上，我略去一些细节的东西，简单地描述一下基本的模型。

现在我们有一个字符序列 $C = c_1c_2\cdots c_T$，并希望得到一个模型 $\mathcal{M}$（依赖于参数 $\theta$）能用来将其词的序列 $S = s_1s_2\cdots s_M$ 。按照生成模型的方式来考虑，将 $C$ 看成是由 $\mathcal{M}$ 生成的序列的话，我们自然希望找到那个生成 $C$ 的可能性最大的 $\mathcal{M}$ ，亦即通过最大似然的方式来估计参数 $\theta$ 。

然而我们不知道似然函数 $\log p(C|\theta)$ 应该如何去优化，因此我们引入 latent variable $S$ ，如果我们知道 $S$ 的话，似然值很容易求得：

\[ \displaystyle \log p(C, S|\theta) = \sum_{i=1}^M \log p(s_i) \]

其中 $p(s_i)$ 的值是从模型 $\mathcal{M}$ 中已知的。但是现在我们不知道 $S$ 的值，因此我们转而取其关于后验概率的期望：

\[ \displaystyle \sum_S p(S|C, \theta^{\text{old}}) \log p(C, S|\theta) \]

然后将这个期望针对 $\theta$ 最大化即完成了 EM 的一次迭代。具体的做法通常是先把一个初始文本（比如所有的 $C$ 的集合）按照 N-gram 分割（N-gram 在 讲 K-medoids 的那篇文章中介绍过）为 $\{w_i\}_N$，形成一个最初的辞典，而模型 $\mathcal{M}$ 的参数 $\theta$ 实际上就描述了各个 N-gram 的概率 $p(w_i)$ ，初始值可以直接取为频率值。有了辞典之后对于任意的 $C$ ，我们可以根据辞典枚举出所有可能的分割 $S$ ，而每个分割的后验概率 $p(S|C, \theta)$ 就是其中的单词的概率乘积。其他的就是标准的 EM 算法的迭代步骤了。

实际上在实际产品中我们会使用一些带了更多启发式规则的分词方法（比如 MMSEG），而这个方法更大的用处在于从一堆文本中“学习”出一个词典来（也就是 $\mathcal{M}$ ），并且这个过程是全自动的，主要有两个优点：

1. 不需要人参与手工分词、标记等。
2. 能自动从文本中学习，换句话说，你给它一些某个领域的专业文章，它能从中学习到这个领域的专业词汇来。

不管怎样，以上两点看起来都非常诱人。不过理论离实际往往还是有不少差距的。我不知道实际分词系统中有没有在用这样的方法来做训练的。之前我曾经用 Wikipedia （繁体中文）上抓下来的一些数据跑过一下小规模的试验，结果还可以，但是也并不如想像中的完美。因为当时也并没有弄明白 EM 是怎么一回事，加上这个过程本身计算负荷还是非常大的，也就没有深入下去。也许以后有时间可以再尝试一下。

总之，EM 是一种应用广泛的算法，虽然讲了点理论也讲了例子，但是一没有贴代码二没有贴图片，似乎实在是不像我的作风，不过精力和篇幅有限，也只能到此为止了。 😉