EM算法详解

作者: oskor | 来源:发表于2019-05-04 23:21 被阅读0次

作为N大机器学习方法的一员，EM算法在各种书籍、博客、网上视频上被描述或者介绍，每次看完总感觉很多地方含糊不清，不能让一个初学者（有一定统计概率基础）接受。最近再B站上，看到徐亦达老师的课程，EM算法这块讲解易于理解和接受，再结合PRML一书的关于混合模型和EM章节内容，对整个EM算法从具体的原理上面有了更深入的理解。
在下文中，更多的是通过公式推导和一些文字说明来梳理EM算法，尽量做到大家一看就明白。

极大似然估计

极大似然估计是概率统计中，估计模型参数的一种方法，如果对似然概念以及极大似然估计的概念不理解，可参考维基百科https://zh.wikipedia.org/wiki/%E6%9C%80%E5%A4%A7%E4%BC%BC%E7%84%B6%E4%BC%B0%E8%AE%A1

这里，我们给出一个一维高斯函数的例子。如图1所示，一维空间里面离散分布的样本点其分布特点可以看成是一种高斯分布。那这些样本的高斯分布函数参数怎么求解呢？可以通过极大似然估计获得。

图1 一维高斯示例

假设我们获取图1中数据样本集合为 $\{x_1,x_2,..x_i,...,x_n\}$ ，其分布函数假设为一维高斯分布，即：
$p(x) =\frac{1}{\sqrt{2\pi}\sigma^2}exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$
那所有数据的联合概率,即似然函数为：
$P(\boldsymbol{x})=\prod_{i=1}^n p(x_i)=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}exp\left(-\frac{(x_i-\mu)^2}{2\sigma^2}\right)$
对等式两边求log，即log-likelihood：
$Ln(P(\boldsymbol{x}))=-nLn(\sqrt{2\pi}\sigma)-\sum_{i=1}^n(x_i-\mu)^2/(2\sigma^2)$
分别对参数求导：
$\frac{\partial{Ln(P(\boldsymbol{x}))}}{\partial\mu}=\frac{\sum_{i=1}^n(x_i-\mu)}{\sigma^2}=0$
$\frac{\partial{Ln(P(\boldsymbol{x}))}}{\partial\sigma^2}=-\frac{n}{2\sigma^2}+\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^4}=0$
可以得到：
$\mu=\frac{1}{n}\sum_{i=1}^n{x_i}$
$\sigma^2=\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2$
通过上述方法，就可以得到基于样本数据和假设分布函数模型情况下，得到样本数据的分布函数。

高斯混合模型

从图2所示中，如果样本数据分布式蓝色点和橙色点的集合，如果依然用高斯分布去进行最大似然估计，得到的分布模型自然是不合适的。很显然，样本数据分布存在两个密集区（蓝色点和橙色点），如果能通过一种方法，确认样本数据里面的蓝色点和橙色点，然后分别对蓝色点集合进行一个高斯估计，对橙色点集进行另外一个高斯估计，两个高斯混合到一起，是不是比单个高斯能更好的匹配样本数据的分布情况。这种情况下的分布函数就是高斯混合模型。

图2 一维混合高斯分布示意

这里我们先直接给出高斯混合模型的分布函数(多维)：
$p(\boldsymbol{x})=\sum_{k=1}^K\pi_k\mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$
在图2例子中，提到如果给每一个数据样本能够先进行分类，即确定属于哪一个集中的簇中，就能比较容易的进行混合模型的求解。这说明了什么呢？我们可以理解，原始样本数据是不完全的（incomplete），引入一个K维的二值随机变量 $\boldsymbol z$ ，这个变量采用"1-of-K"的表示方法，即K维中，只有一维是1，其余所有元素等于0。那么每一个样本数据，即有数据特征 $\boldsymbol x$ ，再匹配一个分配变量 $\boldsymbol z$ ，就可以将图2过程完整描述，即我们认为 $\boldsymbol x$ 和 $\boldsymbol z$ 联合在一起才是完全的（complete）。

数学表示上，我们利用边缘概率分布 $p(\boldsymbol z)$ 和条件概率分布 $p(\boldsymbol x|\boldsymbol z)$ 定义联合概率分布 $p(\boldsymbol x,\boldsymbol z)$ .
$\boldsymbol z$ 的边缘概率分布根据混合系数 $\pi_k$ 进行赋值，即 $p(z_k = 1)=\pi_k,0\leq\pi_k\leq 1,\sum_{k=1}^K\pi_k=1$ ,使得边缘概率分布是合法，那么：
$p(\boldsymbol x) = \prod_{k=1}^K\pi_k^{z_k}$
类似的，在给定 $\boldsymbol z$ 的一个特定的值，即针对某一簇的情况， $\boldsymbol x$ 的条件概率分布是一个高斯分布，即 $p(\boldsymbol x|z_k = 1)=\mathcal N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$
$p(\boldsymbol x|\boldsymbol z) = \prod_{k=1}^K\mathcal N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)^{z_k}$
那么根据贝叶斯公式得到高斯混合模型：
$p(\boldsymbol x) = \sum_{\boldsymbol z}p(\boldsymbol z)p(\boldsymbol x|\boldsymbol z)=\sum_{k=1}^K\pi_k\mathcal N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$
由于我们只能观察样本数据 $\boldsymbol x_1,...,\boldsymbol x_n$ ,无法直接获取每个数据的分配变量 $\boldsymbol z_1,...,\boldsymbol z_n$ ,可理解 $\boldsymbol z_1,...,\boldsymbol z_n$ 为潜在变量（latent）
依据极大似然函数的求解方法，我们可以对高斯混合模型写出数据的对数似然函数：
$Ln(p(\boldsymbol X))=\sum_{i=1}^nLn\left\{ \sum_{k=1}^K\pi_k\mathcal N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\right\}$

由于对数函数内部出现求和，这种情况是没有办法通过求导等于0的方式获取估计参数的解析解的。这种情况下，就需要用到EM算法，通过迭代方式获取估计参数。下面我们先从一般性问题上进行EM算法的理论描述，然后再利用EM算法推导高斯混合模型的计算方法。

EM算法理论证明

EM算法叫做期望最大化方法，首先我们给出EM算法一般性结论或者说步骤，其具体分为两步，即E-step和M-step。

E-step：是求期望的步骤，求谁的期望呢？对联合概率分布的log数据利用潜在变量的条件概率函数求期望，这个期望定义为 $Q(\theta,\theta^{(g)})$ ，我们写连续形式，离散形式类似。
$Q(\theta,\theta^{(g)})=\int_{\boldsymbol z}log(p_{\theta}(\boldsymbol x,\boldsymbol z))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
注：有些书籍或者博客中会将参数写在 $p$ 里面，我个人看着比较不习惯，总会把变量和参数搞混了，这里写在下标上。但如果有了解贝叶斯估计的，那个时候参数也被考虑成变量，其有自己的分布。表达只是一种形式，背后的理论理解掌握就OK啦。
M-step:是求最大的步骤，求谁的最大呢？自然是期望 $Q$ 的最大，从而得到新的估计参数，即
$\theta^{(g+1)}=\mathop{\arg\min}_{\theta}Q(\theta,\theta^{(g)})$

EM算法的步骤，通过高斯混合模型可以直观理解记忆。 $p(\boldsymbol z|\boldsymbol x)$ 是什么意思呢，其含义是在给定数据样本的情况下，潜在变量的概率情况。也就是说在高斯混合模型中，给定样本数据，我们推测样本属于哪一个高斯簇的概率情况。在确定分配情况后，每一个簇都用高斯进行估计，衡量估计好还是不好，用期望形式表示。这样可以帮助理解和记忆Q的定义。那EM算法怎么证明其收敛性呢？

我们要保证：
$p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$
这样每次迭代可以使得似然函数随着参数变大，一直到似然函数不再增大或满足其他终止条件止。

那怎么保证呢？首先，利用贝叶斯公式有：
$p_ {\theta}(\boldsymbol x)=\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{p_\theta(\boldsymbol z|\boldsymbol x)}$

两边同时取log
$log(p_ {\theta}(\boldsymbol x))=log(p_\theta(\boldsymbol x,\boldsymbol z))-log(p_\theta(\boldsymbol z|\boldsymbol x))$
然后两边同时用 $p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)$ 求期望，可以得：
$\int_{\boldsymbol z}log(p_ {\theta}(\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z=\int_{\boldsymbol z}log(p_\theta(\boldsymbol x,\boldsymbol z))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z-\int_{\boldsymbol z}log(p_\theta(\boldsymbol z|\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
等式左边 $p_ {\theta}(\boldsymbol x)$ 和 $\boldsymbol z$ 没有关系， $p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)$ 是概率密度函数，积分是1，所以等式左边依然是 $log(p_ {\theta}(\boldsymbol x))$ .
等式右边第一项就是E-step里面的Q函数，第二项我们记做H函数，则：
$log(p_ {\theta}(\boldsymbol x))=Q(\theta,\theta^{(g)})-H(\theta,\theta^{(g)})$
要保证 $p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$ ，首先 $Q(\theta^{(g+1)},\theta^{(g)})\geq Q(\theta^{(g)},\theta^{(g)})$ ，那么是不是保证 $H(\theta^{(g+1)},\theta^{(g)})\leq H(\theta^{(g)},\theta^{(g)})$ 满足，就一定有 $p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$ ？答案是肯定的。（说明一下，这里面的H和Q函数都是关于 $\theta$ 的函数，每一次迭代 $\theta^{(g)}$ 是已知的，他不是变量）

那下面只要证明 $H(\theta^{(g+1)},\theta^{(g)})\leq H(\theta^{(g)},\theta^{(g)})$ 就可以了。

$H(\theta^{(g)},\theta^{(g)})-H(\theta,\theta^{(g)})$
$=\int_{\boldsymbol z}log(p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z-\int_{\boldsymbol z}log(p_\theta(\boldsymbol z|\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
$=-\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)}\right)p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
$\geq-log\left(\int_{\boldsymbol z}\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)}p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z\right )$
$=-log\left(\int_{\boldsymbol z}p_\theta(\boldsymbol z|\boldsymbol x)d\boldsymbol z\right)=-log(1)=0$

因此可以得到 $H(\theta^{(g+1)},\theta^{(g)})\leq H(\theta^{(g)},\theta^{(g)})$ ，则整个EM算法的收敛性证毕。

注：这里用到了Jessian不等式，如果函数是convex的，则有函数的期望大于期望的函数，即 $E[f(x)]\geq f(E[x])$ .

图3 Jessian不等式示例

EM算法与ELOB和KL散度

上述EM算法的证明，有一些技巧性，而这些技巧性有没有一种是在已知结论基础上硬凑的感觉呢，尤其是用 $p_{\theta^{(g)}}(\boldsymbol z| \boldsymbol x)$ 去求期望，就是为了去构造Q函数。那有没有更好的解释或者更为直观的方式来得到相同的结论呢？答案是有的。

首先，仍然用到：
$log(p_ {\theta}(\boldsymbol x))=log(p_\theta(\boldsymbol x,\boldsymbol z))-log(p_\theta(\boldsymbol z|\boldsymbol x))$
我们稍微变个型，假设一个概率密度函数 $q(\boldsymbol z)$ :
$log(p_ {\theta}(\boldsymbol x))=log\left(\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{q(\boldsymbol z)}\right)-log\left(\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{q(\boldsymbol z)}\right)$
两边同时用 $q(\boldsymbol z)$ 求期望：
$log(p_ {\theta}(\boldsymbol x))=\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{q(\boldsymbol z)}\right)q(\boldsymbol z)d\boldsymbol z-\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{q(\boldsymbol z)}\right)q(\boldsymbol z)d\boldsymbol z$

其中等式右边第一项，我们记做 $\mathcal L(q,\theta)$ ，可以称为ELOB，EvidenceLowerBound，第二项是 $q(\boldsymbol z)$ 和 $p_\theta(\boldsymbol z|\boldsymbol x)$ 的KL散度。

图4 似然函数与ELOB和KL散度

如图4所示，当我固定参数 $\theta$ 时候，ELOB就是关于 $q$ 的泛函（只听过没学过，可理解为函数的函数），那ELOB的上界是什么呢？那首先要理解KL散度，KL散度是衡量两个概率分布之间的差异性，如果两个分布没有差异，则KL散度等于0，如果有差异则大于0，所以KL散度的最小值就是0，那ELOB的上界就是此刻的似然函数。

在EM算法中，当前迭代时，参数 $\theta^{(g)}$ ，则对于E-step，我们需要使得ELOB最大，即KL散度为0，如图5所示，其中 $\theta^{old}$ 为 $\theta^{(g)}$ 。此时， $q(\boldsymbol z) = p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)$

图5 EM算法 E-step

对于M-Step，我们是需要得到新的估计参数，这个时候，固定 $q(\boldsymbol z)$ ,重新获得ELOB的最大值，这个时候的ELOB是什么样的呢？
$\mathcal L(q,\theta)=\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{q(\boldsymbol z)}\right)q(\boldsymbol z)d\boldsymbol z$
$=\int_{\boldsymbol z}log\left(p_\theta(\boldsymbol x,\boldsymbol z)\right)p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z-\int_{\boldsymbol z}log\left(p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)\right)p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
右边第二项没有参数 $\theta$ ，所以固定 $q$ 最大化ELOB，就等于最大化第一项，而第一项是什么呢？就是 $Q$ 函数。

图6 EM算法M-step

如图6所示，我们最大化 $Q$ 函数，也就是最大化此次迭代的ELOB。在新的参数下， $\mathcal L(q,\theta^{(g+1)}) \geq \mathcal L(q,\theta^{(g)})$ ，此刻 $q$ 不变，所以和新的 $p_{\theta^{(g+1)}}(\boldsymbol z|\boldsymbol x)$ 在没有达到局部（或者全局）极大值的情况下，是两个不同的概率分布，即二者KL散度是大于0的，那此刻的似然函数等于此刻的KL散度加上此刻的ELOB，自然是有 $p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$ 。

从ELOB和KL散度的层面可以更好的理解EM算法的迭代过程。

PRML一书中，有图7所示的示意图，能比较形象的说明，EM算法的迭代过程。

图7 EM算法迭代示意

图7中，红色曲线表示(不完整数据)对数似然函数，它的最大值是我们想要得到的。我们首先选择某个初始的参数值 $\theta^{old}$ ,然后在第一个E步骤中,我们计算潜在变量上的后验概率分布,得到了 $\mathcal L(q, \theta^{old})$ 的一个更小的下界，它的值等于在 $\theta^{old}$ 处的对数似然函数值，用蓝色曲线表示。注意，下界与对数似然函数在 $\theta^{old}$ 处以切线的方式连接，因此两条曲线的梯度相同。这个界是一个凹函数，对于指数族分布的混合分布来说，有唯一的最大值。在M步骤中,下界被最大化，得到了新的值 $\theta^{new}$ ,这个值给出了比 $\theta^{old}$ 处更大的对数似然函数值。接下来的E步骤构建了一个新的下界,它在 $\theta^{new}$ 处与对数似然函数切线连接,用绿色曲线表示。以这样的方式不停迭代，直到满足条件。