混合模型与EM算法

混合模型与EM算法

一、前言

在用机器学习方法对数据进行建模时,我们一方面总是希望所用的模型尽量简单,这样对模型相关的各类问题能够直接求出解析解。但简单的模型表达能力有限,无法对复杂的数据特性进行建模。通过将简单模型升级改造为复杂模型可以大大提高模型的表达能力,但往往要面临维度灾难问题,很有可能当你将模型升级为可以表达期望的数据关系时,因为维度灾难问题模型已经变得不可求解了。因此,我们希望找到一种方法,既能保持简单模型的简单求解性,又能使模型具有丰富的表达能力。混合模型即是这样一种方法,它以简单模型作为基础模块,通过合理的引入隐藏变量,将基础模块组合在一起形成具有丰富表达能力的“组合式”模型——混合模型(Mixture Models)。这是进行数据建模非常常用的一种方法,《隐马尔科夫模型》一文即详细介绍了此类模型的一个实例,它展示了如何利用混合模型的思想在控制模型复杂度的同时表达变量之间长距离的相关性。除了隐马尔科夫模型之外,我们平常听到的很多模型其本质都是混合模型,例如K均值算法,混合高斯模型,以及线性动态系统(卡尔曼滤波、粒子滤波)等等。因为混合模型是一种引入隐变量的建模方法,因此也叫隐变量模型,根据引入的隐变量是离散的还是连续的,隐变量模型又分离散隐变量模型连续隐变量模型,前面说的线性动态系统就是连续隐变量的例子。这篇文章主要探讨离散隐变量模型的思想和本质,以及离散隐变量模型的最大似然问题的求解方法——EM算法的算法原理。

本文先从传统视角介绍K均值聚类算法和混合高斯模型的最大似然问题。随后通过引入隐藏变量,解释K均值聚类算法和混合高斯模型为何是离散隐变量模型的实例,并借助求解隐变量模型的最大似然问题的有力工具——EM算法重新推导K均值聚类算法和混合高斯模型的最大似然解,从而揭示EM算法和传统方法的关联。最后,从更一般化的角度证明EM算法为何能够收敛到最大似然问题的极大值。增加K均值聚类相关的章节是为了展示离散隐变量模型的广泛应用,以及用更多的实例解释EM算法,从而加深读者对EM算法的理解。如果您对K均值聚类不感兴趣,完全可以跳过相关章节(第二、七节)。虽然您也可以跳过混合高斯模型的相关章节(第三、四、六节),直接阅读有关EM算法的理论介绍(第五、八节),但我还是建议您不要这样做,因为混合高斯模型是引入EM算法非常合适的例子,有具体的例子做辅助,会更容易理解EM算法理论。如果您已经对混合高斯模型相关章节的内容非常熟悉了,此时便可直接参考EM算法理论介绍的章节。

二、K均值聚类

作为混合模型最简单的一个实例,K均值聚类在很多领域都有广泛的应用。本节先从传统的机器学习视角介绍K均值算法,然后在深入介绍混合模型和EM算法之后,再解释K均值算法如何是混合模型的一个实例,以及如何用EM算法的理论求解K均值聚类问题。

K均值聚类是要解决这样一类问题:给定 N 个数据点 \left \{ \bm{x}_1, \bm{x}_2, \cdots, \bm{x}_N\right \} ,这些点的维数不重要,可以是一维的也可以是多维的,我们的目标是要将这些点分成 K 组,使得同一组内的点之间更“相近”或“相似”,不同组之间的点“距离”更远,一般假设 K 是事先给定的值。为了用更形式化的语言描述这个问题,我们需要引入一些记号。为了表达 K 个分组,我们引入 K 个与数据点维数相同的向量 \bm{\mu}_k ,叫做分组的“原型(prototype)”,我们可以暂且将其理解为分组的聚类中心。为了表达将每个数据点分给哪个分组,对每一个数据点 \bm{x}_n ,我们引入一个对应的二值指示变量 \bm{\gamma}_n ,该变量为 K 维的,每个元素的取值只能是0或1,且同时只能有一个元素取1,其它都取0。如果 \gamma_{n,k} = 1 ,即表示将 \bm{x}_n 分配给第 k 个分组。

如果用欧式距离表示数据点之间的“距离”或“不相似性”的话,那么K均值聚类问题的目标就是,找到 K 个分组的聚类中心 \bm{\mu}_k ,并将每个数据点分配到其中一个分组,使得所有数据点到其所属的聚类中心的距离之和最小,目标函数为

J = \sum_{n=1}^{N} { \sum_{k=1}^{K} {\gamma_{n,k} ||\bm{x}_n - \bm{\mu}_k || ^ 2}} (1)

该问题可用迭代方法求解:首先,假定 \bm{\mu}_k 已知且固定,让 J\gamma_{n,k} 求最小。因为 J\gamma_{n,k} 的线性函数,且不同 n 之间相互独立,因此固定 n 为某一个值,则有由 ||\bm{x}_n - \bm{\mu}_k||^2 给定的 K 个距离值,且每个距离值有一个二值变量 \gamma_{n,k} 作为系数。因为这 K 个二值变量系数同时只能有一个取1,因此,为了使得最终的加权距离和最小,只能让距离最小者的系数取1。也就是说,对每一个数据点 \bm{x}_n ,我们将其分配给与之距离最小的那个分组:

\gamma_{n,k} = \{ \begin{align} & 1, \quad if \ k = \arg \min_{j} || \bm{x}_n - \bm{\mu}_j||^2 \\ &0, \quad otherwise. \end{align} (2)

然后,假定 \gamma_{n,k} 已知且固定,让 J\bm{\mu}_k 求最小。因为 J\bm{\mu}_k 的二次函数,因此对其求导并令导数为0即可:

\frac{\partial{J}}{\partial{\bm{\mu}_k}} = \sum_{n}{2 \gamma_{n,k} (\bm{x}_n - \bm{\mu}_k)} = 0 (3)

容易求解得到
\bm{\mu}_k = \frac{\sum_{n}{\gamma_{n,k}\bm{x}_n}}{\sum_{n}{\gamma_{n,k}}} (4)

该表达式分母表示分配给第 k 个分组的数据点个数,因此整个表达式表示将聚类中心设置为分配给该分组的所有数据点的平均值。正因为这个原因,该算法被叫做K均值聚类算法。

以上两个步骤交互迭代进行,因为每一步都会减小目标值,因此算法最终收敛到极小值。因为以上介绍中我们用欧式距离表达两点间距离,因此该极小值一定为全局最小值。但如果使用更复杂的距离表达式,该算法有可能收敛到局部极小值。有研究者专门研究了K均值算法的收敛问题(MacQueen 1967),有兴趣的读者可以参考相关研究成果。

如果我们用更一般的形式 V(\bm{x}, \bm{x^{'}}) 表示两个数据点之间的距离,则目标函数将变为

\tilde{J} = \sum_{n=1}^{N} {\sum_{k=1}^{K}{\gamma_{n,k} V(\bm{x}_n, \bm{\mu}_k)}} (5)

此即K中心点算法(K-medoids algorithm)。第一步与K均值算法一样,也是将每个数据点分配到与之距离最近的分组。而第二步稍微复杂一点,尤其当距离函数过于复杂时往往无法直接求解。此时一般限制聚类中心必须为该分组中的某一个数据点,因此求聚类中心的问题就简化为选择该分组中的某一个点,使得当该点为聚类中心时类内距离和最小。该简化版本只要求已知 \bm{x}\bm{x^{'}} 时,式 V(\bm{x}, \bm{x^{'}}) 可以求值即可。

三、混合高斯模型:表示

本节介绍混合模型另一个应用非常广泛的实例——混合高斯模型,它比K均值算法稍微复杂一点,但能带领我们更接近EM算法的神秘面貌。混合高斯模型一般用来对数据进行概率建模,跟上一节一样,我们还是先从传统的机器学习视角,建立最大似然问题,通过纯数学的公式推导,得到该最大似然问题的迭代求解方法。随后逐步引入EM算法,在深入介绍EM算法之后,将解释该迭代算法如何被装进EM算法的框架,即如何从EM算法的角度推导出相同的结论。

高斯分布因其简单、常见(大部分自然规律都与高斯分布有关)的特性,在机器学习领域有着非常广泛的应用。然而单独的高斯分布由于太简单,表达能力有限,面对复杂的实际问题往往无能为力。如果我们取多个高斯分布,并以合适的系数加权平均,就能得到表达能力更丰富一类模型——高斯混合模型:

p(\bm{x}) = \sum_{k = 1}^{K} {\pi_k N(\bm{x} | \bm{\mu}_k, \Sigma_k)} (6)

其中 K 为高斯分布的个数, \bm{\mu}_k\Sigma_k 为第 k 个高斯分布的均值向量和方差, \pi_k 为第 k 个高斯分布的加权系数,也叫混合系数。同样的,我们不关心 \bm{x} 的维数,可以是一维的,也可以是多维的,姑且用 D 表示 \bm{x} 的维数。因为概率密度函数要满足非负和归一化的条件,因此混合系数满足

\begin{align} 0 \le \pi_k \le 1 \\ \sum_{k = 1}^{K}{\pi_k} = 1\end{align} (7)

接下来,我们通过引入隐藏变量,构造混合高斯模型的另一种表示方法。我们可以这样想象一个符合混合高斯模型分布的数据点的产生过程:首先以一定的概率从 K 个高斯分布中选择一个,其中第 k 个高斯分布被选中的概率为 \pi_k ,然后依据被选中的高斯分布生成一个数据点。重复以上过程生成多个数据点,这些数据点将符合(6)式的混合高斯分布。为了表示选中哪一个高斯分布,我们引入一个二值指示变量 \bm{z} ,该变量为 K 维的,每个元素的取值只能是0或1,且同时只能有一个元素取1,其它都取0。如果 z_{k} = 1 ,即表示选择第 k 个高斯分布。由于第 k 个高斯分布被选中的概率为 \pi_k ,即

p(z_k=1) = \pi_k (8)

因此,随机变量 \bm{z} 的分布可表示为

p(\bm{z}) = \prod_{k = 1} ^ {K} {\pi_k^{z_k}} (9)

所以,当已知选中第 k 个高斯分布时, x 的条件分布为

p(\bm{x} | z_k = 1) = N(\bm{x} | \bm{\mu}_k, \Sigma_k) (10)

或表示为

p(\bm{x} | \bm{z}) = \prod_{k = 1} ^ {K} {N(\bm{x} | \bm{\mu}_k, \Sigma_k) ^ {z_k}} (11)

从而, \bm{x} 的边缘分布为

p(\bm{x}) = \sum_{\bm{z}} {p(\bm{x} | \bm{z}) p(\bm{z})} = \sum_{\bm{z}} {\prod_{k = 1} ^ {K} { \left[ \pi_k N(\bm{x} | \bm{\mu}_k, \Sigma_k)\right] ^ {z_k}}} (12)

上式等号右边又是求和、又是连乘积,看似很复杂,其实非常简单。因为 \bm{z} 总共只有 K 种取值,对每一个取值,求和符号里面的连乘积的结果只剩下一项,因为 \bm{z} 的每一个取值( K 维二值向量)只有一个元素值为1。例如, 当 \bm{z} 取第 k 种取值( z_k = 1 )时,连乘积只剩下 \pi_k N(\bm{x} | \bm{\mu}_k, \Sigma_k) 一项,因此(12)式简化为

p(\bm{x}) = \sum_{\bm{z}} {p(\bm{x} | \bm{z}) p(\bm{z})} = \sum_{k = 1} ^ {K} {\pi_k N(\bm{x} | \bm{\mu}_k, \Sigma_k)} (13)

可见,用隐藏变量重新表达混合高斯模型后,最终得到的 \bm{x} 的边缘分布与(6)式一模一样,这是必要的。因为无论你用什么方式重新表达混合高斯模型,最后 \bm{x} 的边缘分布一定要“回归”到(6)式,否则你的表达是毫无意义的,或者说是错误的。

如果将(6)式看做混合高斯模型的“判别式”表示方法,则(8)~(13)式可以看做混合高斯模型的“产生式”表示方法,因为这是从数据产生的角度对模型进行的表示。

从目前来看,我们用隐藏变量的方法重新表示混合高斯模型,似乎没什么收益。事实上,从后面的推导将会看出,这种表示方法能够大大简化最大似然问题的求解过程。因为如果直接用(6)式建模,似然函数将会是高斯函数和的连乘积形式,求解此种形式的最优化问题将会是非常困难的。但是引入隐藏变量后,我们就能够使用联合分布 p(\bm{x}, \bm{z}) 简化问题,因为联合分布直接是高斯函数的连乘积形式,这种形式便于求解析解,我们从接下来的推导中将会看到。在此之前,我们先看看变量 \bm{z} 的条件概率的形式。如果用 \gamma(z_k) 表示条件概率 p(z_k = 1 | \bm{x}) ,利用贝叶斯定理,

\begin{align} \gamma(z_k) = p(z_k = 1 | \bm{x}) &= \frac{p(\bm{x} | z_k = 1) p(z_k = 1)}{\sum_{j = 1} ^{K} {p(\bm{x} | z_j = 1) p(z_j = 1)}} \\ &= \frac{\pi_k N(\bm{x} | \bm{\mu}_k, \Sigma_k)}{\sum_{j = 1} ^{K} {\pi_j N(\bm{x} | \bm{\mu}_j, \Sigma_j)}} \end{align} (14)

上式是在已知观测变量 \bm{x} 的条件下 z_k = 1 的后验概率,而其先验概率已经知道为 \pi_k

接下来我们就可以形式化地表示混合高斯模型的最大似然问题了。已知 N 个数据点 \left \{ \bm{x}_1, \bm{x}_2, \cdots, \bm{x}_N\right \} ,我们希望用混合高斯模型对其进行概率密度建模。对每一个数据点 \bm{x}_n ,都有一个与之对应的二值指示变量 \bm{z}_n 。为了记号的紧凑,我们用一个 N \times D 的矩阵 X 表示所有的数据点,其中矩阵的第 n 行为 \bm{x}_n^T 。类似地,我们用一个 N \times K 的矩阵 Z 表示所有的二值指示变量,其中矩阵的第 n 行为 \bm{z}_n^T 。如果我们假设给定的数据集独立同分布,那么根据(6)式,对数似然函数为

p(X | \bm{\pi}, \bm{\mu}, \Sigma) = \sum_{n = 1} ^ {N} {\ln { \{ \sum_{k = 1} ^ {K} {\pi_k N (\bm{x} | \bm{\mu}_k, \Sigma_k) \}}}} (15)

在求解该最大似然问题之前,我们需要首先讨论一下应用最大似然方法求解混合高斯模型参数时将要面临的两个问题。第一个是奇点问题。简单起见,假设某一个高斯分布的协方差矩阵为对角阵,即 \Sigma_k = \sigma^2 I ,其中 I 为单位阵。需要说明的是,不仅仅是对角阵,奇点问题对任意形式的协方差矩阵都是存在的。假设这个高斯分布的均值向量正好等于某一个数据点,即 \bm{\mu}_k = \bm{x}_n ,则该数据点在这个高斯分布下的似然值为

N(\bm{x}_n | \bm{\mu}_k, \sigma^2 I) = \frac{1}{(2 \pi) ^{D / 2}} \frac{1}{ { \sigma}} (16)

\sigma \to 0 时,上式将趋近于无穷大,因此似然函数也将趋近于无穷大。因此混合高斯模型的最大似然问题不是一个良好定义(well-posed)的问题,因为一旦其中某一个高斯分布“坍缩”到某个数据点上,似然值将变为无穷大,此时求得的参数将是不可信的。这可以看做是所有最大似然法求解面临的过拟合问题。这个问题并不是无法解决的,以后我们将会看到,通过采用贝叶斯方法可以解决最大似然求解的过拟合问题,不过暂时我们不打算这么做。因此,在应用最大似然法求解混合高斯模型的参数时,需要注意避免奇点问题带来的病态解,可以在迭代求解的过程中,不断检测是否有高斯分布坍缩到某个数据点,倘若如此,则将其均值重新初始化为随机值,并将其协方差重置为某个较大的值,然后继续迭代。

第二个问题是唯一性问题。这个问题很好理解,假定我们已经求得了混合高斯模型的一组解,然后不改变这组解的值,仅仅将 K 个参数重新排列组合,则可以得到另一组等价的解。这样的组合共有 K! 种。以 K = 3 的一维模型为例,假设我们已经求得一组解: \bm{\pi} = [0.5, 0.3, 0.2] ^ T\mu_1 = 1\mu_2 = 2\mu_3 = 3\Sigma_1 = 0.1 ^ 2\Sigma_2 = 0.2 ^ 2\Sigma_3 = 0.3 ^ 2 。如果交换前两组参数的顺序,可以得到另一组等价的解: \bm{\pi} = [0.3, 0.5, 0.2] ^ T\mu_1 = 2\mu_2 = 1\mu_3 = 3\Sigma_1 = 0.2 ^ 2\Sigma_2 = 0.1 ^ 2\Sigma_3 = 0.3 ^ 2 。这些不同组合的等价解共有 3! = 6 种。唯一性问题在需要对求解的参数进行解释时需要认真对待,但如果仅仅是为了找到数据集的概率密度模型,那么最终得到哪一组解并不重要,因为任意一组解都是等价的。

四、混合高斯模型:求解

求解的目的是要对(15)式求最大值,并求得取最大值时所有参数的取值。为了简化求导过程,一般将原始问题转化为对数似然的最大化问题,对(15)式两边取对数

\ln {p(X | \bm{\pi}, \bm{\mu}, \Sigma)} = \sum _ {n = 1} ^ {N} {\ln \{ \sum _ {k = 1} ^ {K} {\pi_k N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)} \}} (17)

我们不妨先看看(17)式取极大值时参数应该满足什么条件。

先让(17)式对 \bm{\mu}_k 求导并令导数为0:

\frac{\partial P(X | \bm{\pi}, \bm{\mu}, \Sigma)}{\partial {\bm{\mu}_k}} = -\sum_{n = 1} ^ {N} {\frac{\pi_k N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)}{\sum_{j = 1} ^ {K} {\pi_j N(\bm{x}_n | \bm{\mu}_j, \Sigma_j)}} \Sigma_k^{-1} (\bm{x}_n - \bm{\mu}_k)} = 0 (18)

上式最外层求和符号里面的第一项分式恰好就是由(14)式定义的后验概率,仅仅多一个下标 n ,可以简记为 \gamma(z_{nk}) ,带入上式容易求解得到

\bm{\mu}_k = \frac{\sum_n{\gamma(z_{nk})\bm{x}_n}}{\sum_n{\gamma(z_{nk})}} = \frac{1}{N_k}\sum_n{\gamma(z_{nk})\bm{x}_n} (19)

其中,定义了

N_k = \sum_n{\gamma(z_{nk})} (20)

由于 \gamma(z_{nk}) 可以理解为数据点 \bm{x}_n 来自第 k 个高斯分布的可能性,或者第 k 个高斯分布对生成该数据点承担的“责任(responsibility)”,因此 N_k 可以理解为所有数据点中来自第 k 个高斯分布的数据点的“有效个数”。从而(19)式拥有良好的解释性,它表示,第 k 个高斯分布的均值 \bm{\mu}_k 等于所有数据点的加权平均,其中加权系数为该高斯分布应当为产生该数据点承担的责任。

注意到,

\begin{align} & \sum_k{\gamma(z_{nk})} = 1 \\ & \sum_n{\sum_{k} \gamma(z_{nk})} = N \end{align} (21)

上式的结论随后将会用到。

然后让(17)式对 \Sigma_k 求导并令导数为0,

\frac{\partial P(X | \bm{\pi}, \bm{\mu}, \Sigma)}{\partial {\Sigma_k}} = - \frac{1}{2} \sum _ {n = 1} ^ {N} {\gamma(z_{nk}) \left[ |\Sigma_k| ^ {-1} \frac{\partial{|\Sigma_k|}}{\partial{\Sigma_k}} + \frac{\partial{(\bm{x}_n - \bm{\mu}_k) ^ T \Sigma_k^{-1} (\bm{x}_n - \bm{\mu}_k)}}{\partial{\Sigma_k}}\right]} = 0 (22)

要求解该式,需要一些矩阵求导相关的知识。由于矩阵求导不是本文的重点,此处直接给出相关结论,对详细推导过程感兴趣的读者可以参考相关资料。

\frac{\partial{|A|}}{\partial{A}} = |A| (A ^ {-1})T (23)

\frac{\partial{A ^ {-1}}}{\partial{x}} = -A^{-1} \frac{\partial{A}}{\partial{x}} A^{-1} (24)

将(23)、(24)两式带入(22)式将其进一步简化为

- \frac{1}{2} \sum_{n = 1} ^ {N} {\gamma(z_{nk}) \Sigma_k^{-1} \left[ I - (\bm{x}_n - \bm{\mu}_k) (\bm{x}_n - \bm{\mu}_k)^T \Sigma_k^{-1} \right]} = 0 (25)

求解得到

\Sigma_k = \frac{1}{N_k} \sum_n {\gamma(z_{nk}) (\bm{x}_n - \bm{\mu}_k) (\bm{x}_n - \bm{\mu}_k)^T} (26)

与(19)式类似,上式可解释为,第 k 个高斯分布的协方差矩阵 \Sigma_k 等于所有数据点以 \bm{\mu}_k 为中心的协方差矩阵的加权平均,其中加权系数为该高斯分布应当为产生该数据点承担的责任。

最后,让(17)式对 \pi_k 求极大值。因为 \pi_k 满足(7)式的归一化约束,因此,在原目标函数的基础上增加一个拉格朗日乘子,

\ln{p(X | \bm{\pi}, \bm{\mu}, \Sigma) } + \lambda (\sum_k{\pi_k} - 1) (27)

让上式对 \pi_k 求导并令导数为0:

\sum_{n = 1} ^ {N} {\frac{N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)}{\sum_{j = 1} ^ {K} {\pi_j N(\bm{x}_n | \bm{\mu}_j, \Sigma_k)}}} + \lambda = 0 (28)

等式两边同时乘以 \pi_k

\sum_n{\gamma(z_{nk})} + \pi_k \lambda = 0 (29)

上式两边对 k 求和,并利用(7)式和(21)式的结论可得

\lambda = -N (30)

带入(29)式得

\pi_k = \frac{N_k}{N} (31)

因此上式可以解释为,混合系数 \pi_k 等于第 k 个高斯分布为产生所有数据点承担的责任的平均值。

注意到,(19)、(26)、(31)式并不构成该最大似然问题的解析解,因为 \bm{\mu}_k\Sigma_k\pi_k 的求解依赖 \gamma(z_{nk}) ,而 \gamma(z_{nk}) 的求解又反过来依赖 \bm{\mu}_k\Sigma_k\pi_k ,详见(14)式。不过这很容易让我们想到用迭代的方法求得参数值。第一步,随机选取 \bm{\mu}_k\Sigma_k\pi_k 的初始值,依据(14)式求得 \gamma(z_{nk}) 。第二步,将第一步求得的 \gamma(z_{nk}) 依次带入(19)、(26)、(31)式更新 \bm{\mu}_k\Sigma_k\pi_k 的值。如此往复迭代,直至达到收敛条件(似然值或参数值的变化在某个阈值以下)。因为每一步的迭代都会增加似然值,因此该迭代算法可以保证找到似然函数的一个极大值,但不能保证找到全局最大值。迭代算法的收敛速度很大程度上受初始值选取的影响,为了加快算法的收敛速度,一般先用K均值算法求得 \bm{\mu}_k\Sigma_k\pi_k 的初始值,然后以此初始值开始迭代。

下图是用包含两个高斯分布的混合高斯模型对old faithful数据集进行概率密度建模的示例图。图中,红色(R)圆圈和蓝色(B)圆圈分别为两个高斯分布在一个标准差处的等高线,数据点用红色和蓝色的混合表示,红色的多少用该数据点属于红色高斯分布的概率表示,蓝色的多少用该数据点属于蓝色高斯分布的概率表示。因此,如果一个点完全属于蓝色,则该点将用蓝色画出,如果完全属于红色,则用红色画出,如果属于红色和蓝色的概率恰好相等,则用紫色画出。如图(1),初始时,随机设置两个高斯分布的均值和方差,由于此时还未计算数据点属于哪个高斯分布的概率,因此所有数据点用蓝色(G)画出(红色和蓝色通道为0)。如图(2),经过一轮EM迭代后,两个高斯分布的均值和方差均得到了更新,且每个数据点属于两个高斯分布的概率也计算出来了。可以看到,靠近蓝色高斯分布中心的数据点几乎完全为蓝色,靠近红色高斯分布中心的数据点几乎完全为红色,而处在两个高斯分布中间的那些点带着紫色,说明这些点属于哪个高斯分布还不好确定。此时两个高斯分布的方差还比较大,对数据点的概率密度的拟合还比较模糊。如图(3),经过5轮EM迭代后,两个高斯分布的均值和方差开始逐渐收敛。本例设置的收敛条件为数据集的对数似然值的变化量小于1e-4,经过11轮EM迭代后算法收敛,最终收敛结果如图(4)所示。图(5)给出了数据集的对数似然随迭代次数的变化趋势,可见对数似然值在开始几轮迭代中迅速增加,随后到达平台期逐步上升最终达到收敛条件。本例中参数的初始值是随机选取的,因此收敛速度较慢。如果先用K均值算法找到一组“差不多”的参数,并以此作为初始参数开始迭代,可以大大加快收敛速度。文末给出了本例所使用的数据集和matlab代码。

图(1) 初始化状态
图(2) 经过第1论EM迭代后,两个高斯分布的均值和方差均得到了更新。
图(3) 经过5轮EM迭代后,模型已经能较好地拟合观测数据。
图(4) 经过11轮EM迭代后达到收敛。
图(5) 数据集的对数似然随迭代次数的变化趋势。

至此,我们还没有引入EM算法求解混合高斯模型的最大似然问题,不过已经很接近了。EM算法也是一种迭代算法,每一步迭代分E步和M步两步。从接下来的介绍将会看到,上述迭代算法的第一步和第二步恰好对应EM算法的E步和M步。

五、EM算法

此节将会脱离具体问题背景,介绍一般化的EM算法。EM算法用于求解带隐藏变量的模型的最大似然问题。EM的全拼为,Expectation Maximization,即期望最大化算法。

为此首先按照上文类似的方式引入几个记号:用 N \times D 的矩阵 X 表示所有的观测变量,其中矩阵的第 n 行为 \bm{x}_n^T ,数据点的维数为 D ;用 N \times K 的矩阵 Z 表示所有的隐藏变量,其中矩阵的第 n 行为 \bm{z}_n^T ;用 \theta 表示所有的模型参数。因此,对数似然函数为

\ln {p(X | \theta)} = \ln {\{\sum_{Z}{p(X, Z | \theta)}\}} (32)

观察上式可以发现,对数符号直接作用在联合分布的和上,而不是直接作用于联合分布。这使得要求得上式最优化问题的解析解几乎是不可能的。

假设给定数据集 X ,对应的 Z 是知道的,那么 XZ 的联合分布的对数似然函数为 \ln {p(X, Z | \theta)} ,一般来说,最大化该似然函数的解析解是比较好求的。我们把 \{X, Z\} 叫做完全数据集,而把单独的 \{ X\} 的叫做不完全数据集。

但实际中,我们只知道 XZ 是不知道的,所有关于 Z 的信息只能从其后验分布 p(Z | X, \theta) 中得知,因此我们无法求得全数据的似然。不过我们可以考虑全数据的似然在后验分布 p(Z | X, \theta) 下的期望,此即期望最大化算法中“期望”二字的来历,对应EM算法的E步。然后,我们将此期望值最大化,就可以求得一组参数值,此即期望最大化算法中“最大化”一词的来历,对应EM算法的M步。从而,EM算法可被描述为包含如下两步的迭代过程。在E步,用参数的当前值 \theta^{old} 求得隐藏变量的后验分布 p(Z | X, \theta^{old}) ;在M步,用此后验分布计算全数据的似然的期望,将此期望记作 Q(\theta, \theta^{old}) ,则

Q(\theta, \theta^{old}) = \sum_{Z} {p(Z | X, \theta^{old}) \ln p(X, Z | \theta)} (33)

然后求得此期望取最大值时的参数

\theta^{new} = \arg\max_{\theta} {Q(\theta, \theta^{old})} (34)

注意到,(33)式中,对数符号直接作用于联合分布 p(X, Z | \theta) 上,因此上述最大化问题一般比较容易求得解析解。

将新求得的 \theta^{new} 作为下一步的 \theta^{old} 重复以上步骤,直至算法收敛。

读到此处,相信的大家心中都有一个疑问,为什么取全数据对数似然的期望并最大化就可以最大化(32)式的似然呢?这个操作看上去似乎毫无依据啊?难道随便蒙了一个恰好就蒙对了?当然不是蒙的,而且这样选择是有明确的理论依据的,问题的答案要在下文证明EM算法时给出。此时,我们先给出EM算法一般化的算法框架。

EM算法: 给定全数据的似然 p(X, Z | \theta) 表达式,目的是求得参数 \theta 使得观测数据的似然 p(X | \theta) 最大。 1. 选择参数的初始值,记为 \theta^{old}
2. E步:计算隐藏变量的后验分布 p(Z | X, \theta^{old}) 3. M步:最大化全数据的对数似然在 p(Z | X, \theta^{old}) 下的期望,从而得到新的参数值 \theta^{new} : \theta^{new} = \arg\max_{\theta} {Q(\theta, \theta^{old})} 其中, Q(\theta, \theta^{old}) Q(\theta, \theta^{old}) = \sum_{Z} {p(Z | X, \theta^{old}) \ln p(X, Z | \theta)} 4. 检查迭代是否收敛(似然值或参数值的变化量是否在给定阈值以下);如果收敛,则停止,算法结束;否则,令 \theta^{old} \leftarrow \theta^{new} ,返回第2步。


、EM视角下的混合高斯模型求解

有了上一节的知识基础,本节我们应用EM算法框架重新求解混合高斯模型的最大似然问题,从而加深对EM算法的理解。可以看到,采用EM算法将会得到与第四节相同的结论。

EM算法要求计算全数据的对数似然。由式(9)和(11)可以得到,全数据的似然函数

\begin{align} p(X, Z | \bm{\pi}, \bm{\mu}, \Sigma) &= \prod_{n = 1} ^ {N} {p(\bm{x}_n, \bm{z}_n | \bm{\pi}, \bm{\mu}, \Sigma)} \\ &= \prod_{n = 1} ^ {N} {p(\bm{x}_n | \bm{z}_n, \bm{\mu}, \Sigma) p(\bm{z}_n | \bm{\pi})} \\ &= \prod_{n = 1} ^ {N} {\prod_{k = 1} ^ {K} {\left[ \pi_k N(\bm{x}_n | \bm{\mu}_k, \Sigma_k) \right] ^ {z_{nk}}}} \end{align} (35)

对两边取对数得

\ln {p(X, Z | \bm{\pi}, \bm{\mu}, \Sigma)} = \sum_{n = 1} ^ {N} {\sum_{k = 1} ^ {K} { z_{nk} \left[ \ln\pi_k + \ln N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)\right]}} (36)

因此,全数据的对数似然在后验分布 p(Z | X, \theta^{old}) 下的期望为

\begin{align} Q(\theta, \theta^{old}) &= \sum_{Z} {p(Z | X, \theta^{old}) \sum_{n = 1} ^ {N} {\sum_{k = 1} ^ {K} { z_{nk} \left[ \ln\pi_k + \ln N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)\right]}}} \\ &= \sum_{n = 1} ^ {N} {\sum_{k = 1} ^ {K} { \sum_{Z} {p(Z | X, \theta^{old})} z_{nk} \left[ \ln\pi_k + \ln N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)\right]}} \end{align} (37)

由于给定 X 时各个隐藏变量之间相互独立,因此

\begin{align} \sum_{Z} {p(Z | X, \theta^{old})} z_{nk} &= \sum_{Z} {\prod_{n = 1} ^ {N} {p(\bm{z}_n | \bm{x}_n, \theta^{old}) z_{nk}}} \\ &= \sum_{\bm{z}_n} {p(\bm{z}_n | \bm{x}_n, \theta^{old}) z_{nk}} \\ &= p(\bm{z}_{nk} = 1 | \bm{x}_n, \theta^{old}) \\ &= \gamma(z_{nk}) \end{align} (38)

从而,

Q(\theta, \theta^{old}) = \sum_{n = 1} ^ {N} {\sum_{k = 1} ^ {K} { \gamma(z_{nk}) \left[ \ln\pi_k + \ln N(\bm{x}_n | \bm{\mu}_k, \Sigma_k)\right]}} (39)

对上式求极大值就比较容易了,可以直接得到解析解。首先,让上式对 \bm{\mu}_k 求导并令导数为0,可以得到与式(19)相同的结论。然后,让上式对 \Sigma_k 求导并令导数为0,可以得到与式(26)相同的结论。最后,考虑 \pi_k 的归一化约束条件,增加相应的拉格朗日乘子后令导数为0,可以得到与式(31)相同的结论。

七、EM视角下的K均值聚类算法

前文分别介绍了K均值算法和混合高斯模型,本节将会看到K均值算法其实是混合高斯模型的一个特例,只需要对混合高斯模型增加一个简单的约束即可。因此,也可以将EM算法应用到混合高斯模型的结论应用于K均值算法,可以看到,用EM算法求解K均值聚类问题的结果与第二节介绍的直接最小化聚类误差得到的结果一模一样。

从前文的介绍可以看到,K均值算法和混合高斯模型其实有一些相似之处。例如,都假设模型由 K 个“组件”构成,K均值聚类假设有 K 个聚类,混合高斯模型假设有 K 个高斯分布;都从产生式的角度将数据点分配给各个组件,只不过K均值聚类将每个数据点唯一分配给其中一个组件,具体分配给哪个组件由 \gamma_{nk} 指定((2)式),而混合高斯模型认为每个数据点可以从属于每个高斯分布,而从属的程度由后验分布 \gamma(z_{nk}) 给定。如果说K均值聚类是“硬”分配(hard assignment)的话,那么混合高斯模型则是“软”分配(soft assignment)。

在混合高斯模型中,如果假设所有高斯分布的协方差矩阵相同,且都为 \epsilon{I} ,其中 \epsilon 为已知的常数,则此时

\gamma(z_{nk}) = \frac{\pi_k\exp\{-||\bm{x}_n-\bm{\mu}_k||^2/(2\epsilon)\}}{\sum_{j = 1} ^ {K} {\pi_j \exp\{ -||\bm{x}_n - \bm{\mu}_j|| ^ 2 / (2\epsilon)\}}} (40)

现在让 \epsilon \to 0 ,则上式分子分母中的指数项都将趋近于0,但趋近的速度不一样,其中 ||\bm{x}_n - \bm{\mu}_j||^2 最小的那一项趋近的速度最慢。因此,只有当 k = \arg \min_{j} {||\bm{x}_n - \bm{\mu}_j||^2} 时, \gamma(z_{nk}) = 1k 为其他值时 \gamma(z_{nk}) 均为0。即

\gamma(z_{nk}) = \{ \begin{align} &1, \quad if \ k = \arg \min _{j} {||\bm{x}_n - \bm{\mu}_j||^2} \\ &0, \quad otherwise \end{align} (41)

此式与(2)式一模一样。可见,此时 \gamma(z_{nk}) \to \gamma_{nk} ,即混合高斯模型的软分配退化成了K均值聚类的硬分配。同时,(19)式的均值更新方程退化成(4)式的聚类中心更新方程。从而可以看到,混合高斯模型退化成了K均值聚类算法。

八、EM算法的证明

前文一直在讲EM算法就是先求隐藏变量的后验分布,然后最大化全数据的对数似然在此后验分布下的期望,这样不断迭代就可以得到EM算法的解。但是相信读者心中还带着疑问,为什么要选择全数据的对数似然?为什么要最大化全数据的对数似然在后验分布下的期望?为什么这样迭代后却会收敛到观测数据的似然 p(X | \theta) 的极大值?

要回答以上问题很简单,不过需要对观测数据的似然做一个非常巧妙的变换:

\begin{align} \ln p(X | \theta) &= \sum_{Z} {q(Z) \ln p(X | \theta)} \\ &= \sum_{Z}{q(Z) \ln \frac{p(X,Z | \theta)}{p(Z|X, \theta)}} \\ &= \sum_Z{q(Z) \ln \frac{p(X, Z | \theta)}{q(Z)}} + \sum_{Z} {q(Z)\ln\frac{q(Z)}{p(Z | X, \theta)}} \\ &= L(q, \theta) + KL(q || p) \end{align} (42)

其中,

L(q, \theta) = \sum_Z{q(Z) \ln \frac{p(X, Z | \theta)}{q(Z)}} (43)

KL(q || p) = \sum_{Z} {q(Z)\ln\frac{q(Z)}{p(Z | X, \theta)}} (44)

因此,将目标函数 \ln p(X | \theta)分解成了 L(q, \theta)KL(q||p) 两项。注意到, L(q, \theta)\theta 的函数,是 q 的泛函,因为 q 本身是一个函数。 KL(q||p) 是分布 qp 之间的KL距离,也是一个泛函,它具有非负的特性,即 KL(q||p) \ge 0 ,当且仅当 q = p 时, KL(q||p) = 0 。此分解的示意图如下:

因此,根据(42)式的分解,可以推导出最大化 \ln p(X | \theta) 的迭代算法,即EM算法。

在E步,假设参数为已知的常数,记为 \theta^{old} ,因此(42)式的左边为一个常数,而等式右边分解成两项,这两项的大小根据 q 的取值不同而不同。因为 KL(q||p) \ge 0 ,所以 L(q, \theta)\ln p(X | \theta) 的下界,我们要让这个下界尽量大。显然,只有当 KL(q||p) = 0q(Z) = p(Z | X, \theta ^ {old}) 时,该下界最大,且正好等于 \ln p(X|\theta^{old}) 。如下图,红色线的高度表示等式左边的值 \ln p(X | \theta^{old}) ,此时它是一个常量。蓝色虚线为分解的位置,虚线下面是 L(q, \theta^{old}) 的值,虚线上面是 KL(q||p) 的值。根据 q 的取值不同蓝色虚线可以上下移动,当 q(Z) = p(Z | X, \theta ^ {old}) 时,蓝色虚线移动到与红色线等高,此时下界 L(q, \theta^{old}) 达到最大。

在接下来的M步,固定 q(Z) 的取值,以 \theta 为参数最大化 L(q, \theta) ,得到新的参数 \theta^{new} 。因为 L(q, \theta)\ln p(X | \theta) 的下界, 下界增大了, \ln p(X | \theta) 也会跟着增大。且此时因为 \theta^{new} \neq \theta^{old} ,从而 q \neq p,所以 KL(q||p) \gt 0 。因此 \ln p(X | \theta) 增大的量比下界增大的量更大。如下图,红色虚线和蓝色虚线分别是经过E步之后目标函数及其下界的值,此时二者相等。在经过M步之后,下界增大到蓝色实线的位置,同时 KL(q||p) 不再等于0,所以目标函数从红色虚线增大到红色实线的位置,增大的量等于两项分解式增大的量的和。

q(Z) = p(Z | X, \theta ^ {old}) 带入(43)式,可以看到我们在M步要最大化的下界变为

\begin{align} L(q, \theta) &= \sum_{Z} {p(Z | X, \theta^{old}) {p(X, Z | \theta)}} - \sum_{Z} {p(Z | X, \theta^{old}) {p(Z | X, \theta^{old})}} \\ &= Q(\theta, \theta^{old}) + const \end{align} (45)

因此,在M步要最大化的实际是全数据的对数似然在后验分布 p(Z | X, \theta ^ {old}) 下的期望,这正是前文在用EM算法求解混合高斯模型时给出的结论。

因为每经过一轮E步和M步,都会增大目标函数的值,因此当迭代收敛时,找到的一定是目标函数的极大值。

EM算法的核心步骤在于求解M步的最大化问题,当 xz 的联合分布属于指数分布族时,对数正好和指数相消,可以较容易求得解析解。但当联合分布是其他复杂的分布时,M步往往不好直接求解,此时需要对标准EM算法进行扩展,其中GEM(Generalized EM)和ECM(Expectation conditional maximization)算法即是其中的例子。在数据量比较大无法一次性载入内存进行计算时,还可以将标准EM算法转化为序列化算法,一次只需要载入一个数据点对参数进行序列化地更新。除此之外,EM算法还可以用来求解参数 \theta 的MAP问题。这些内容由于不在本文的范围之内,就不再赘述了。此文完。

附matlab源码

用包含两个高斯分布的混合高斯模型对old faithful数据集进行概率密度建模的matlab源码如下:

% main.m
clear
clc

load old_faithful
[sof, z_mu, z_sigma] = zscore(old_faithful); % scaled old faithful data

[N, D] = size(sof);
K = 2;
X = sof;
% initialize. You can also initialize parameters using K-means algorithm.
pie = [0.5, 0.5];
mu = [-1.0, 0; 1 1.5];
sigma(:, :, 1) = [0.5, 0; 0, 0.5];
sigma(:, :, 2) = [0.5, 0; 0, 0.5];
index = 1;
% EM iterations
while 1
    % plot before E step (the label of each point is unknown and so ploted
    % in green color.
    figure
    plot(X(:, 1), X(:, 2), 'g.');
    hold on
    h = ezplot(@(x,y)myfunc(x, y, mu(1, :), inv(sigma(:, :, 1))));
    set(h, 'color', 'b', 'lineWidth', 2)
    h = ezplot(@(x,y)myfunc(x, y, mu(2, :), inv(sigma(:, :, 2))));
    set(h, 'color', 'r', 'lineWidth', 2)
    title(sprintf('iteration %d, before E step.', index));
    axis equal
    axis([-2 2 -3 3])
    hold off
    % E step -- calculate responsibilities
    R = zeros(N, K);
    for k = 1 : K
        R(:, k) = mvnpdf(X, mu(k, :), sigma(:, :, k)) * pie(k);
    end
    LL(index) = sum(log(sum(R, 2))); % log likelihood
    R = R ./ repmat(sum(R, 2), 1, K);
    % M step -- update parameters
    Nk = sum(R);
    for k = 1 : K
        mu(k, :) = sum(diag(R(:, k)) * X) / Nk(k);
        MU = repmat(mu(k, :), N, 1);
        sigma(:, :, k) = (X - MU)' * diag(R(:, k)) * (X - MU) / Nk(k);
    end
    pie = Nk / N;
    % plot after M step
    figure
    hold on
    for i = 1 : N
        plot(X(i, 1), X(i, 2), '.', 'color', [R(i, 2) 0 R(i, 1)]);
    end
    h = ezplot(@(x,y)myfunc(x, y, mu(1, :), inv(sigma(:, :, 1))));
    set(h, 'color', 'b', 'lineWidth', 2)
    h = ezplot(@(x,y)myfunc(x, y, mu(2, :), inv(sigma(:, :, 2))));
    set(h, 'color', 'r', 'lineWidth', 2)
    title(sprintf('iteration %d, after M step', index));
    axis equal
    axis([-2 2 -3 3])
    hold off
    mu
    sigma
    pie
    LL(index)
    if index > 1 && LL(index) - LL(index - 1) <= 1e-4
        break;
    end
    index = index + 1;
end


function [ z ] = myfunc( x, y, mu, inv_sigma )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
z = (x - mu(1)) .* inv_sigma(1, 1) .* (x - mu(1)) ...
  + (x - mu(1)) .* inv_sigma(1, 2) .* (y - mu(2)) ...
  + (y - mu(2)) .* inv_sigma(2, 1) .* (x - mu(1)) ...
  + (y - mu(2)) .* inv_sigma(2, 2) .* (y - mu(2)) - 1;

end

old faithful数据集可以从网上公开的地址下载,或者从本人以下网盘地址下载:

链接: pan.baidu.com/s/1nw0kNF 密码: ujnb