机器学习面试必知：高斯混合模型GMM和期望最大化EM算法

通过将更基本的概率分布（例如高斯分布）进行线性组合的这样的叠加方法，可以被形式化为概率模型，被称为混合模型。通过使用足够多的高斯分布，并且调节它们的均值和方差以及线性组合的系数，几乎所有的连续概率密度都能以任意的精度去近似。

我们考虑K个高斯概率密度的叠加，形式为$$p(x)=\sum _{k=1}^K{\pi }_{k}N(x|u_k,{\Sigma }_k)$$这就被称为混合高斯。参数${\pi }_k$被称为混合系数并且${\sum }_{k=1}^K{\pi }_k=1$。根据概率的加和规则和乘积规则，边缘概率密度为$$p(x)=\sum _{k=1}^{K}p(k)p(x|k)$$其中我们把${\pi }_k=p(k)$看成选择第k个成分的先验概率，把密度$N(x|u_k,{\Sigma }_k)=p(x|k)$看成以k为条件的概率。后验概率$p(k|x)$起着一个重要作用，它也被成为责任。根据贝叶斯定理，我们可以写出后验概率$${\gamma }_k(x)=p(k|x)=\frac{p(k)p(x|k)}{{\sum }_{l}p(l)p(x|l)}=\frac{{\pi }_{k}N(x|u_k,{\Sigma }_k)}{{\sum }_l{\pi }_{l}N(x|u_l,{\Sigma }_l)}$$很自然地对数似然函数可以表示为$$\mathrm{l}\mathrm{n}p(X|\pi ,u,\Sigma )=\sum _{n=1}^N\mathrm{l}\mathrm{n}\left\{\sum _{k=1}^K{\pi }_{k}N(x_n|u_k,{\Sigma }_k)\right\}$$其中 X = { x 1 , . . . , x n } X=\left \{ x_{1},...,x_{n}\right \} X={x1​,...,xn​}, π = { π 1 , . . . , π k } \pi=\left \{ \pi_{1},...,\pi_{k}\right \} π={π1​,...,πk​}, u = { u 1 , . . . , u k } u=\left \{ u_{1},...,u_{k}\right \} u={u1​,...,uk​}, Σ = { Σ 1 , . . . , Σ k } \Sigma=\left \{ \Sigma_{1},...,\Sigma_{k}\right \} Σ={Σ1​,...,Σk​}。因为对数中存在一个求和式，这会导致参数的最大似然解不再有一个封闭式的解析解。一个方法就是期望最大化（EM）的强大框架。
我们先来看看高斯混合模型是如何得到的
我们引入一个K维二值随机变量z，这个变量采用了1-of-K表示方法，其中一个特定的元素$z_k=1$，其余所有的元素等于0。于是$z_k$的值满足 z k ∈ { 0 , 1 } z_{k}\in \left\{ 0,1\right\} zk​∈{0,1}。我们根据边缘概率分布$p(z)$和条件概率分布$p(x|z)$来定义联合概率分布。z的边缘概率分布根据混合系数${\pi }_k$进行赋值$$p(z_k=1)={\pi }_k$$$$0\le {\pi }_k\le 1,\sum _{k=1}^K{\pi }_k=1$$由于$z$使用了1-of-K表示方法，因此我们可以将这个概率分布写成$$p(z)=\prod _{k=1}^K{\pi }_k^{z_k}$$类似地，给定z一个特定的值，x的条件概率分布是一个高斯分布$$p(x|z_k=1)=N(x|u_k,{\Sigma }_k),p(x|z)=\prod _{k=1}^{K}N(x|u_k,{\Sigma }_k{)}^{z_k}$$x的边缘概率分布可以通过将联合概率分布对所有可能的z求和的方式得到$$p(x)=\sum _{z}p(z)p(x|z)=\sum _z\prod _{k=1}^K({\pi }_{k}N(x|u_k,{\Sigma }_k){)}^{z_k}=\sum _{j=1}^K\prod _{k=1}^K({\pi }_{k}N(x|u_k,{\Sigma }_k){)}^{I_{kj}}=\sum _{k=1}^K{\pi }_{k}N(x|u_k,{\Sigma }_k)$$因此x的边缘概率分布就是高斯混合分布。另一个起着重要作用的变量是$\gamma (z_k)$表示的是在给定x的条件下z的条件概率$p(z_k=1|x)$$$\gamma (z_k)=p(z_k=1|x)=\frac{p(z_k=1)p(x|z_k=1)}{{\sum }_{j=1}^{K}p(z_j=1)p(x|z_j=1)}=\frac{{\pi }_{k}N(x|u_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_{j}N(x|u_j,{\Sigma }_j)}$$类似地我们将${\pi }_k$看成$z_k=1$的先验概率，将$\gamma (z_k)$看成观测到x之后，对应的后验概率。正如我们将看到的那样，$\gamma (z_k)$也可以看成是分量k对于解释观测值x的责任。从上面可以看到似然函数为$$\mathrm{l}\mathrm{n}p(X|\pi ,u,\Sigma )=\sum _{n=1}^N\mathrm{l}\mathrm{n}\left\{\sum _{k=1}^K{\pi }_{k}N(x_n|u_k,{\Sigma }_k)\right\}$$根据以往一般的经验，我们接下来需要最大化这个似然函数。在这之前我们需要强调一下由于奇异性造成的应用于高斯混合模型的最大似然框架中的一个大问题。
假设我们有这样一个数据集 X = { x 1 , . . . , x n } X=\left \{ x_{1},...,x_{n}\right \} X={x1​,...,xn​}N×D的矩阵，其中第n行是$x_n^T$。对应的隐含变量表示成N×K的矩阵Z，第n行为$z_n^T$。我们考虑这样一个高斯混合模型，它的分量协方差矩阵为${\Sigma }_k={\sigma }_k^{2}I$,$I$是一个单位矩阵。假设混合模型的第j个分量的均值$u_j$与某个数据点完全相同即$u_j=x_n$,这样这个数据点会为似然函数贡献一项$$N(x_n|x_n,{\sigma }_j^{2}I)=\frac{1}{(2\pi {)}^{0.5}}\frac{1}{{\sigma }_j^D}$$较为致命的是一旦${\sigma }_j\to 0$,那么这一项将会趋于无穷大，因此对数似然函数也会趋向于无穷大。

因此，对数似然函数的最大化不是一个具有良好定义的问题，因为这种奇异性总会出现，会发生在任何一个“退化”到一个具体的数据点上的高斯分量上。而这个问题为什么不会在单一的高斯分布中出现呢，我们注意到如果单一的高斯分布退化到一个数据点上，那么它总会给其他数据点产生的似然函数贡献可乘的因子，这些因子会以指数的速度趋向于零，从而使得整体的似然函数趋向于零而不是无穷大。然而，一旦我们混合概率分布中存在两个分量，其中一个分量会具有有限的方差，因此对所有的数据点都会赋予一个有限的概率值，而另一个分量会收缩到一个具体的数据点，因此会给对数似然函数贡献一个不断增加的值。更一般的，我们可以把这种奇异性当做是最大似然方法中不可避免的过拟合现象。可以参照下图理解，右边的对应着即其中一种分量会收缩到一个具体的数据点，这样会形成一种病态解。

而运用贝叶斯的方法能有效地解决过拟合问题即能解决这种奇异困难。而一般的方法是如果检测到某种分量收缩到一个点时，我们就重新随机设置它的均值，并且将它的方差设置为某个较大的值，再继续优化。
不仅如此，寻找最大似然解还存在另外一种问题，我们想象得到K个高斯分量混合而成的概率分布总共存在K！种等价的解，因为任意交换其中两种分量的解还是会得到和未交换前一样的效果。
总的来说最大化高斯混合模型的对数似然函数很复杂，从公式上体现在求和项出现在对数计算中，所以不会像单一的高斯分布那么好求解。如果我们令其导数为0，也无法得到一个解析解。
一种优雅的并且强大的寻找带有潜在变量的模型的最大似然解的方法被称为期望最大化（EM）

首先我们先使得似然函数关于高斯分量的均值$u_k$的导数等于零即$$\frac{\partial \mathrm{l}\mathrm{n}p}{\partial u_k}=0\to \sum _{n=1}^N\frac{{\pi }_{k}N(x_n|u_k,{\Sigma }_k)}{{\sum }_j{\pi }_{j}N(x_n|u_j,{\Sigma }_j)}{\Sigma }_k^{-1}(x_n-u_k)=0$$其中$$N(x_n|u_k,{\Sigma }_k)=\frac{1}{(2\pi {)}^{\frac{D}{2}}}\frac{1}{|{\Sigma }_k{|}^{0.5}}exp\left\{-\frac{1}{2}(x-u{)}^T{\Sigma }^{-1}(x-u)\right\}$$$$\frac{{\pi }_{k}N(x_n|u_k,{\Sigma }_k)}{{\sum }_j{\pi }_{j}N(x_n|u_j,{\Sigma }_j)}=\gamma (z_{nk})$$$$\sum _{n=1}^N\gamma (z_{nk}){\Sigma }_k^{-1}(x_n-u_k)=0\to u_k=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})x_n$$$$N_k=\sum _{n=1}^N\gamma (z_{nk})$$
我们可以将$N_k$看作是分配到聚类k的数据点的有效数量。均值$u_k$通过对数据集里所有的数据点求加权平均的方式得到，其中$x_n$的权因子由后验概率$\gamma (z_{nk})$给出。可以与单一高斯分布的最大似然估计进行对比$$u_{ML}=\frac{1}{N}\sum _{n=1}^N{x}_n$$$${\Sigma }_{ML}=\frac{1}{N}\sum _{n=1}^N(x_n-u_{ML})(x_n-u_{ML}{)}^T$$
类似地我们先使得似然函数关于高斯分量的协方差矩阵${\Sigma }_k$的导数等于零$${\Sigma }_k=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})(x_n-u_k)(x_n-u_k{)}^T$$
最后我们对混合系数${\pi }_k$求导数，同时要考虑它的限制条件${\sum }_{k=1}^K{\pi }_k=1$，我们使用拉格朗日乘数法，最大化下面这个函数$$L=\mathrm{l}\mathrm{n}p(X|\pi ,u,\Sigma )+\lambda (\sum _{k=1}^K{\pi }_k-1)$$$$\frac{\partial L}{\partial {\pi }_k}=0\to \sum _{n=1}^N\frac{N(x_n|u_k,{\Sigma }_k)}{{\sum }_j{\pi }_{j}N(x_n|u_j,{\Sigma }_j)}+\lambda =0$$
$$\sum _{n=1}^N\frac{{\pi }_{k}N(x_n|u_k,{\Sigma }_k)}{{\sum }_j{\pi }_{j}N(x_n|u_j,{\Sigma }_j)}=-\lambda {\pi }_k\to \sum _{n=1}^N\frac{{\sum }_k{\pi }_{k}N(x_n|u_k,{\Sigma }_k)}{{\sum }_j{\pi }_{j}N(x_n|u_j,{\Sigma }_j)}=-\lambda \sum _j{\pi }_k\to \lambda =-N$$我们最后会得到$${\pi }_k=\frac{N_k}{N}$$ 从而第k个分量的混合系数是那个分量对于解释数据点的“责任”的平均值。
我们可以看到上述这三个参数$u_k$,${\Sigma }_k$,${\pi }_k$并没有给出一个解析解，因为责任$\gamma (z_{nk})$也依赖于这些参数。那么知道了$\gamma (z_{nk})$不就能求出这三个参数了吗。所以我们能使用一种简单的迭代方法来寻找问题的最大似然解。而这个迭代过程就是EM算法应用于高斯混合模型的一个实例。我们先为这三个参数设置一个初始值，然后用上面得到的关系式进行交替更新，被称为E步骤和M步骤。每次进行E步骤和接下来的M步骤都保证了对数似然函数的增大。实际应用中，当对数似然函数的变化量或者参数的变化量低于阈值时，我们就认为算法收敛。通常运行K均值算法找到高斯混合模型的一个合适的初始值，接下来使用EM算法进行调节。协方差矩阵通过K均值算法找到的聚类样本协方差得到，混合系数可以设置为分配到对应类别中的数据点所占比例。算法必须避免似然函数带来的奇异性，即高斯分量退化到一个具体的数据点，并且一般情况下对数似然函数会有多个局部极大值，EM并不能保证找到其中最大的一个。

给定一个高斯混合模型，目标是关于参数(均值，协方差，混合系数)最大化似然函数

• 初始化均值$u_k$,协方差${\Sigma }_k$,混合系数${\pi }_k$,计算对数似然函数的初始值
• E步骤，使用当前参数计算“责任”$$\gamma (z_{nk})=\frac{{\pi }_{k}N(x_n|u_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_{j}N(x_n|u_j,{\Sigma }_j)}$$
• M步骤，使用当前的责任更新估计参数$$u_k^{new}=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})x_n$$$${\Sigma }_k^{new}=\frac{1}{N_k}\sum _{n=1}^N\gamma (z_{nk})(x_n-u_k^{new})(x_n-u_k^{new}{)}^T$$$${\pi }_k^{new}=\frac{N_k}{N}，N_k=\sum _{n=1}^N\gamma (z_{nk})$$
  -计算对数似然函数，看似然函数和参数的变化量是否收敛，没有收敛的话继续进行E步骤和M步骤。$$\mathrm{l}\mathrm{n}p(X|\pi ,u,\Sigma )=\sum _{n=1}^N\mathrm{l}\mathrm{n}\left\{\sum _{k=1}^K{\pi }_{k}N(x_n|u_k,{\Sigma }_k)\right\}$$

更一般的EM算法推导如下：
假设我们有n个样本数据 X = { x 1 , . . . , x n } X=\left \{ x_{1},...,x_{n}\right \} X={x1​,...,xn​}，找出样本的模型参数$\theta$极大化模型分布的对数似然函数$$\theta =arg\underset{\theta }{\max }\sum _{i=1}^{n}logP(x_i;\theta )$$如果我们有未观测到的隐含数据 Z = { z 1 , . . . , z n } Z=\left \{ z_{1},...,z_{n}\right \} Z={z1​,...,zn​}，那么相应的对数似然函数如下$$\theta =arg\underset{\theta }{\max }\sum _{i=1}^{n}log\sum _{z_i=1}^{n}P(x_i,z_i;\theta )$$由于无法直接求解上式，我们需要一种特殊的方法，由于是求的最大值，所以我们可以舍弃一些模型较小的值，所以可以先求出下界$$\sum _{i=1}^{n}log\sum _{z_i=1}^{n}P(x_i,z_i;\theta )=\sum _{i=1}^{n}log\sum _{z_i=1}^n{Q}_i(z_i)\frac{P(x_i,z_i;\theta )}{Q_i(z_i)}\ge \sum _{i=1}^n\sum _{z_i=1}^n{Q}_i(z_i)log\frac{P(x_i,z_i;\theta )}{Q_i(z_i)}$$这里用了引入了一个新的分布Q,并且用了Jessen不等式求下界，满足等号的情况如下$$\frac{P(x_i,z_i;\theta )}{Q_i(z_i)}=c,c=常数$$因为Q是一个分布那么很自然地想到$$\sum _{z_i}{Q}_i(z_i)=1\to \sum _{z_i}P(x_i,z_i;\theta )=c\sum _{z_i}{Q}_i(z_i)=c$$$$Q_i(z_i)=\frac{P(x_i,z_i;\theta )}{c}=\frac{P(x_i,z_i;\theta )}{{\sum }_{z_i}P(x_i,z_i;\theta )}=\frac{P(x_i,z_i;\theta )}{P(x_i;\theta )}=P(z_i|x_i;\theta )$$此时我们只需要优化$$arg\underset{\theta }{\max }\sum _{i=1}^n\sum _{z_i=1}^n{Q}_i(z_i)logP(x_i,z_i;\theta ),其中Q_i(z_i)=P(z_i|x_i;\theta )$$上式的前半部分求最大化就是M步骤，E部分就是$$先求条件概率分布Q_i(z_i)=P(z_i|x_i;\theta )$$$$再求损失函数L(\theta )=\sum _{i=1}^n\sum _{z_i=1}^n{Q}_i(z_i)logP(x_i,z_i;\theta )$$其中${\sum }_{z_i=1}^n{Q}_i(z_i)logP(x_i,z_i;\theta )$可以理解成$logP(x_i,z_i;\theta )$基于条件概率分布$Q_i(z_i)$的期望。

来个实际例题试试2008年的文章What is the expectation maximization algorithm?

可以看到结果基本相似。