高斯混合模型GMM聚类的步骤和推导

0. 引言

  由于最近要做聚类算法方面的内容，看了很多资料，在高斯混合模型(GMM)这里一直没有一个让我完全推导清楚的、理解的文章。经过三天打鱼两天晒网 不懈努力，总算是有一点自己的理解，我希望尽量通俗地把GMM讲明白，同时也希望尽量详细地对公式进行推导和解释。因此，我会先给出GMM算法的总体步骤，保证拿上先可以直接使用，然后再进行具体的推导和解释。文中可能有一些自己理解不严谨的地方，还请大家指正。

1. 算法初窥

  已知样本集是 D = { x 1 , x 2 , . . . , x m } D=\{x_1,x_2,...,x_m\} D={x1​,x2​,...,xm​}，要将这些样本聚成 $k$ 类。我们认为样本服从混合高斯分布：
$$p_M(\mathbf{x})=\sum _{i=1}^k{\alpha }_i\cdot p(\mathbf{x}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})$$
  其中 p ( x ∣ μ i , Σ i ) = 1 ( 2 π ) n / 2 ∣ Σ i ∣ 1 / 2 e x p { − 1 2 ( x − μ i ) T Σ i − 1 ( x − μ i ) } p(\bm{x}|\bm{\mu_i}, \bm{\Sigma_i})=\frac{1}{_{(2\pi)^{n/2}|\bm{\Sigma_i}|^{1/2}}}exp\{-\frac{1}{2}(\bm{x}-\bm{\mu_i})^T\bm{\Sigma_i}^{-1}(\bm{x}-\bm{\mu_i})\} p(x∣μi​,Σi​)=(2π)n/2∣Σi​∣1/2​1​exp{−21​(x−μi​)TΣi​−1(x−μi​)}是一个多元高斯分布，即一个混合成分；
  ${\alpha }_i$表示混合系数，即选择第 $i$ 个混合成分的概率。

第一步 初始化高斯混合分布的模型参数 ${\alpha }_i,{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}}$
第二步 计算$x_j$由各混合成分生成的后验概率，即观测数据$x_j$由第 $i$ 个分模型生成的概率$p(z_j=i|{\mathbf{x}}_{\mathbf{j}})$，并记为${\gamma }_{ji}$
    ${\gamma }_{ji}=\frac{{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})}{{}_{{\sum }_{l=1}^{{}^k}{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})}}$

第三步 计算新的模型参数：

    ${\mathbf{\mu }}_{\mathbf{i}}^{\prime }=\frac{{\sum }_{{}_{j=1}}^m{\gamma }_{ji}{{\mathbf{x}}_{\mathbf{j}}}^{|}}{{\sum }_{{}_{j=1}}^m{\gamma }_{ji}^{|}}$

    ${\mathbf{\Sigma }}_{\mathbf{i}}^{\prime }=\frac{{\sum }_{{}_{j=1}}^m{\gamma }_{ji}({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}}^{\prime })({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}}^{\prime }{)}^{T^{|}}}{{\sum }_{{}_{j=1}}^m{\gamma }_{ji}^{|}}$

    ${\alpha }_i^{\prime }=\frac{{\sum }_{{}_{j=1}}^m{\gamma }_{ji}^{|}}{m}$

第四步 按照新的模型参数重复2，3步，直到满足停止条件
第五步 将每个样本按照 λ j = arg max ⁡ i ∈ { 1 , 2 , . . . , k } γ j i \lambda_j=\argmax\limits_{i\in\{1,2,...,k\}} \gamma_{ji} λj​=i∈{1,2,...,k}argmax​γji​划入对应的簇。即对每个样本来自哪个分模型的概率大就划入哪个分模型的簇中，最终就得到了 $k$ 个聚类

2. 高斯混合模型的引入

  与k-means聚类不同，高斯混合聚类是采用概率模型来刻画聚类结构。实际上我们可以采用任意不同的概率分布模型来进行刻画，高斯分布是最普遍的一种，如下：
  高斯分布：
$$p(\mathbf{x})=\frac{1}{(2\pi {)}^{n/2}|\mathbf{\Sigma }{|}^{1/2}}exp[-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })]$$
  而单高斯分布模型有其局限性，不能完全反映数据分布的特点，因此我们用多个高斯分布的线性叠加来刻画实际样本，其中一个高斯分模型称为一个混合成分。

理论上来说，当叠加的高斯分模型数量足够多时，可以表征任意一种分布。（这其实很好理解，类比足够多微小线段可以逼近任意一条曲线、足够多复指数信号可以描述任意信号…是一样的道理）

  高斯混合分布：
$$p_M(\mathbf{x})=\sum _{i=1}^k{\alpha }_i\cdot p(\mathbf{x}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})$$
我们认为，手里拿到的样本就是根据这个概率分布抽取得到的（或者说“生成的”）
例如，对于第$j$个样本${\mathbf{x}}_{\mathbf{j}}$，就根据 $p_M({\mathbf{x}}_{\mathbf{j}})={\sum }_{i=1}^k{\alpha }_{ji}\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})$ 得到。

$\mathbf{T}\mathbf{I}\mathbf{P}\mathbf{S}:$

1. 这里的$p(x)$、$p_M(x)$指的是概率密度函数，不是概率，在有些概率书上为了区别，用$f(x)$表示，这里都用$p(x)$表示，但心里要清楚其含义。
2. 接上条，所以$p(\mathbf{x}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})$不是条件概率，而是概率密度，"$|{\mu }_i,{\Sigma }_i$"只是明确一下这个概率密度函数包含的参变量。实际上它表示的就是上面单高斯分布的$p(x)$。
3. $\mathbf{x}$是一条样本，但是有$n$个维度，因此是一个$n$维向量。
4. ${\alpha }_i>0$是在生成这条样本时，选择通过第$i$个分模型来生成的概率，且${\sum }_{i=1}^k{\alpha }_i=1$。（不能说成"样本来自第$i$个分模型的概率"，因为这里是一个先验的情况，如果这样说就成了后验了）
5. ${\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i$是第$i$个分模型的参数。其中，${\mathbf{\mu }}_i$表示均值，是一个$n$维向量，${\mathbf{\Sigma }}_i$表示协方差矩阵，是一个$n\times n$方阵。

3. 按照高斯混合模型进行聚类划分

  上面说了我们认为手里拿到的样本就是通过高斯混合模型抽取得到的，那么反过来我们要怎么把这些样本用高斯混合模型划分成不同的类别呢？
  一个很直接的想法自然是按照模型的混合成分划成 $k$ 类，一个数据最可能从哪个分模型得来就认为属于哪一类。
  在这里，我们要引入一个隐变量 z j ∈ { 1 , 2 , . . . , k } z_j\in\{1,2,...,k\} zj​∈{1,2,...,k} 表示得到样本 ${\mathbf{x}}_{\mathbf{j}}$的高斯分模型。

注：

1. 有的书上用一维向量来表示，即若认为样本${\mathbf{x}}_{\mathbf{j}}$来自第2个高斯分模型，则$z_j=[0,1,0,0,...,0]$。 这里直接用数字来表示来自第几个分模型。
2. 根据$z_j$的含义很容易看出，$P(z_j=i)$表示${\mathbf{x}}_{\mathbf{j}}$是通过第$i$个分模型生成的概率，就是高斯混合模型中的参数${\alpha }_{ji}$

  前面我们说了，$\alpha$是一个先验概念，是从模型到样本的过程。而现在我们在已经拿到了样本的情况下反推其来自哪个分模型，是逆向过程，因此我们用$p_M(z_j=i|{\mathbf{x}}_{\mathbf{j}})$来表示样本${\mathbf{x}}_{\mathbf{j}}$来自第$i$个分模型的后验概率，并简记为${\gamma }_{ji}$。有：
$$\begin{array}{rl}{p}_M(z_j=i|{\mathbf{x}}_{\mathbf{j}})& =\frac{P(z_j=i)\cdot p_M({\mathbf{x}}_{\mathbf{j}}|z_j=i)}{p_M({\mathbf{x}}_{\mathbf{j}})}\\ & =\frac{{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})}{{\sum }_{l=1}^k{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})}\end{array}$$

注：

1. $p_M({\mathbf{x}}_{\mathbf{j}}|z_j=i)$表示按照第$i$个高斯分模型生成${\mathbf{x}}_{\mathbf{j}}$的概率密度，第$i$个高斯分模型的参数是${\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}}$，故而就等于$p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})$
2. $p_M({\mathbf{x}}_{\mathbf{j}})$表示综合所有的混合成分后总的概率密度
3. 上述等式第一行由贝叶斯公式得到
    贝叶斯公式： $p(A|B)=\frac{p(A)p(B|A)}{{}_{p(B)}}$

那么显而易见地，每个样本${\mathbf{x}}_{\mathbf{j}}$的簇标记${\lambda }_j$如下确定：
λ j = arg max ⁡ i ∈ { 1 , 2 , . . . , k } γ j i \lambda_j=\argmax_{i \in \{1,2,...,k\}}\gamma_{ji} λj​=i∈{1,2,...,k}argmax​γji​
即，${\mathbf{x}}_{\mathbf{j}}$来自哪个分模型的概率最大，就认为属于哪一类。

4. 确定高斯混合模型参数

  上面已经说了当已知高斯混合模型时，就可以进行聚类的划分，那么如何求解这个模型，得到它的三个参数呢？
  我们在这里要用到的是EM算法（期望最大算法）。其实原理很简单：为什么我们能拿到手中的样本，而不是其他数据呢？我们认为这是由于选出这样一组样本的概率最大，所以才运气爆表，被我们拿到手。
  由上文知，按照高斯混合模型选出一个样本${\mathbf{x}}_{\mathbf{j}}$的概率密度
$$p_M({\mathbf{x}}_{\mathbf{j}})=\sum _{i=1}^k{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})$$
  对于手中的$m$个样本，选到任意一个都是一个独立事件，最终的概率自然是全部相乘，即
$$\prod _{j=1}^m{p}_M({\mathbf{x}}_{\mathbf{j}})$$
但是，连乘不好处理，因此一般习惯对它取对数，于是样本集$D$的最大化对数似然函数就定义如下：
$$\begin{array}{rl}LL(D)& =\ln (\prod _{j=1}^m{p}_M({\mathbf{x}}_{\mathbf{j}}))\\ & =\sum _{j=1}^m\ln (p_M({\mathbf{x}}_{\mathbf{j}}))\\ & =\sum _{j=1}^m\ln (\sum _{i=1}^k{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}}))\end{array}$$

  只要能求出使$LL(D)$最大的参数就可以了。
  那么怎么求满足要求的参数呢？
  我们设参数 θ i = { ( α i , μ i , Σ i ) } \theta_i=\{(\alpha_i,\mu_i,\Sigma_i)\} θi​={(αi​,μi​,Σi​)}能使$LL(D)$最大化，那么$LL(D)$对每个参数的偏导数应该为0，但是偏导数为0求出的参数有可能只是局部最优解（$LL(D)$取极大值或驻点），而不是全局最优解（$LL(D)$取最大值）。
  经过后面的推导，我们可以发现求出的每个参数，都可以用${\gamma }_{ji}$表示。所以，我们在求出了一组模型参数后，按照这种模型得到对应的${\gamma }_{ji}$，再用得到的${\gamma }_{ji}$继续按照偏导数为0的方式求出新的参数。如此循环迭代，直到我们认为足够为止。

至于为什么每次迭代都可以使求得的参数更优，这个问题就不在本文展开叙述了，具体可以参考EM算法的相关资料。

  现在我们来具体求解每个参数：

① $\mathbf{\mu }:$

$$\begin{array}{rl}& \frac{\partial LL(D)}{\partial {\mathbf{\mu }}_{\mathbf{i}}}=0\\ & \Rightarrow \frac{\partial }{\partial {\mathbf{\mu }}_{\mathbf{i}}}\sum _{j=1}^m\ln (\sum _{i=1}^k{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}}))=0\\ & \Rightarrow \sum _{j=1}^m\frac{1}{{\sum }_{l=1}^k{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})}\cdot \frac{\partial }{\partial {\mathbf{\mu }}_{\mathbf{i}}}[\sum _{l=1}^k{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})]=0\end{array}$$

  (这里因为对${\mathbf{\mu }}_{\mathbf{i}}$求偏导，为了避免混淆，将求和变量写成 $l$)
  对$\frac{\partial }{\partial {\mathbf{\mu }}_{\mathbf{i}}}[{\sum }_{l=1}^k{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})]$来说，只有当 $l=i$时，包含${\mu }_i$的内容，其余对${\mu }_i$求偏导均为0，可以舍去，则继续推导如下：

$$\Rightarrow \sum _{j=1}^m\frac{1}{{\sum }_{l=1}^k{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})}\cdot \frac{\partial }{\partial {\mathbf{\mu }}_{\mathbf{i}}}[{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})]=0$$

  其中，

∂ ∂ μ i [ α i ⋅ p ( x j ∣ μ i , Σ i ) ] = ∂ ∂ μ i { α i 1 ( 2 π ) n / 2 ∣ Σ i ∣ 1 / 2 exp ⁡ [ − 1 2 ( x j − μ i ) T Σ i − 1 ( x j − μ i ) ] } = α i exp ⁡ [ − 1 2 ( x j − μ i ) T Σ i − 1 ( x j − μ i ) ] ( 2 π ) n / 2 ∣ Σ i ∣ 1 / 2 ∂ ∂ μ i [ − 1 2 ( x j − μ i ) T Σ i − 1 ( x j − μ i ) ] = α i ⋅ p ( x j ∣ μ i , Σ i ) ⋅ ( x j − μ i ) \begin{aligned} & \frac{\partial}{\partial \bm{\mu_i}}[ \alpha_i \cdot p(\bm{x_j}|\bm{\mu_i},\bm{\Sigma_i})] \\ & =\frac{\partial}{\partial \bm{\mu_i}}\{\alpha_i\frac{1}{(2\pi)^{n/2}|\bm{\Sigma_i}|^{1/2}} \exp[-\frac{1}{2}(\bm{x_j}-\bm{\mu_i})^T\bm{\Sigma_i}^{-1}(\bm{x_j}-\bm{\mu_i})]\} \\ & =\alpha_i\frac{\exp[-\frac{1}{2}(\bm{x_j}-\bm{\mu_i})^T\bm{\Sigma_i}^{-1}(\bm{x_j}-\bm{\mu_i})]}{(2\pi)^{n/2}|\bm{\Sigma_i}|^{1/2}} \frac{\partial}{\partial\bm{\mu_i}}[-\frac{1}{2}(\bm{x_j}-\bm{\mu_i})^T\bm{\Sigma_i}^{-1}(\bm{x_j}-\bm{\mu_i})] \\ & =\alpha_i\cdot p(\bm{x_j}|\bm{\mu_i},\bm{\Sigma_i})\cdot(\bm{x_j}-\bm{\mu_i}) \end{aligned} ​∂μi​∂​[αi​⋅p(xj​∣μi​,Σi​)]=∂μi​∂​{αi​(2π)n/2∣Σi​∣1/21​exp[−21​(xj​−μi​)TΣi​−1(xj​−μi​)]}=αi​(2π)n/2∣Σi​∣1/2exp[−21​(xj​−μi​)TΣi​−1(xj​−μi​)]​∂μi​∂​[−21​(xj​−μi​)TΣi​−1(xj​−μi​)]=αi​⋅p(xj​∣μi​,Σi​)⋅(xj​−μi​)​

这里是向量/矩阵对另一个向量求导，不是标量求导，具体可以参考矩阵求导相关资料

  因此继续推导如下：
$$\begin{array}{rl}& \Rightarrow \sum _{j=1}^m\frac{{\alpha }_i\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{i}},{\mathbf{\Sigma }}_{\mathbf{i}})}{{\sum }_{l=1}^k{\alpha }_l\cdot p({\mathbf{x}}_{\mathbf{j}}|{\mathbf{\mu }}_{\mathbf{l}},{\mathbf{\Sigma }}_{\mathbf{l}})}({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}})=0\\ & \Rightarrow \sum _{j=1}^m{p}_M(z_j=i|{\mathbf{x}}_{\mathbf{j}})({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}})=0\\ & \Rightarrow \sum _{j=1}^m{\gamma }_{ji}\cdot ({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}})=0\\ & \Rightarrow \sum _{j=1}^m{\gamma }_{ji}{\mathbf{x}}_{\mathbf{j}}=\sum _{j=1}^m{\gamma }_{ji}{\mathbf{\mu }}_{\mathbf{i}}\\ & \Rightarrow {\mathbf{\mu }}_{\mathbf{i}}=\frac{{\sum }_{j=1}^m{\gamma }_{ji}{\mathbf{x}}_{\mathbf{j}}}{{\sum }_{j=1}^m{\gamma }_{ji}}\end{array}$$

  至此，参数${\mathbf{\mu }}_{\mathbf{i}}$迭代公式得到。

② $\mathbf{\Sigma }:$

  同理，由
$$\frac{\partial LL(D)}{\partial {\mathbf{\Sigma }}_{\mathbf{i}}}=0$$
  推得：
$${\mathbf{\Sigma }}_{\mathbf{i}}=\frac{{\sum }_{j=1}^m{\gamma }_{ji}({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}})({\mathbf{x}}_{\mathbf{j}}-{\mathbf{\mu }}_{\mathbf{i}}{)}^T}{{\sum }_{j=1}^m{\gamma }_{ji}}$$

③ $\alpha :$

  求$\alpha$的过程略有不同，因为除了要使$LL(D)$最大化以外，$\alpha$还要满足它自身的条件：${\alpha }_i\ge 0,{\sum }_{i=1}^k{\alpha }_i=1$。
  这是一个有条件的极值问题，我们要用拉格朗日乘数法来求解（具体可以参考拉格朗日乘数法求极值的相关资料）
  相当于将$LL(D)$求极值问题转化为$LL(D)+\lambda ({\sum }_{i=1}^k{\alpha }_i-1)$求极值的问题，然后依然对${\alpha }_i$求导为0，由此得到：
$${\alpha }_i=\frac{1}{m}\sum _{j=1}^m{\gamma }_{ji}$$

  至此，高斯混合模型聚类的所有参数公式均已得到，之后只要不断迭代，并按照文章第3节中的划分方式来进行聚类划分即可。

  最后，可以再回头看看文章第1节的算法总结。

---

参考资料

1. 周志华，机器学习，清华大学出版社，2016
2. 李航，统计学习方法，清华大学出版社，2012

---