GMM高斯混合模型原理推导(二)

前言

上一篇GMM高斯混合模型原理推导(一)，我们的连乘已经变成了连加号，只需要求出$P(z,x)$相关的概率即可
数学基础：【概率论与数理统计知识复习-哔哩哔哩】

原理推导

对于$P(z,x)$，为什么不是$P(z,x|\theta )$？因为$\theta$是参数，不是随机变量
$$P(z_i=Ck,x_i)=P(x_i|z_i=Ck)P(z_i=Ck)=p_k\ast N(x_i|{\mu }_k,{\Sigma }_k)$$
而$P(z|x)$我们目前没有必要求出来，因为他的参数是给定${\theta }^t$，我们所需要的是变量，因为后续要求导求极值。而$P(z|x,{\theta }^t)$已经相当于一个常数

因此
$$\begin{array}{c}\begin{array}{rl}{E}_{P(Z|X,{\theta }^t)}\left[logP(Z,X|\theta )\right]=& \sum _{k=1}^K\sum _{i=1}^{n}log[p_k\ast N({\mu }_k,{\Sigma }_k)]P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{k=1}^K\sum _{i=1}^n[log{p}_k+logN({\mu }_k,{\Sigma }_k)]P(z_i=Ck|x_i,{\theta }^t)\end{array}\end{array}$$
先求出$p_k$，对于$p_k$，我们知道有约束条件$\sum _{k=1}^K{p}_k=1$，所以构造拉格朗日函数
$$L(\theta ,\lambda )=\sum _{k=1}^K\sum _{i=1}^n[log{p}_k+logN(x_i|{\mu }_k,{\Sigma }_k)]P(z_i=Ck|x_i,{\theta }^t)+\lambda \left[\sum _{k=1}^K{p}_k-1\right]$$
让其对关于$p_k$求导
$$\begin{array}{c}\begin{array}{rl}& \frac{\partial L(\theta ,\lambda )}{\partial p_k}=\sum _{i=1}^n\frac{1}{p_k}P(z_i=Ck|x_i,{\theta }^t)+\lambda =0\\ & 等式左右乘以p_k\\ & \sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)+\lambda p_k=0\end{array}\end{array}$$
因此，当$k=1,2,\cdots ,K$时
$$\begin{gathered} \sum _{i=1}^{n}P(z_i=C1)+\lambda p_1=0 \\ \sum _{i=1}^{n}P(z_i=C2)+\lambda p_2=0 \\ ⋮ \end{gathered}$$
所以
$$\sum _{i=1}^{n}P(z_i=C1)+\lambda p_1+\sum _{i=1}^{n}P(z_i=C2)+\lambda p_2+\cdots +\sum _{i=1}^{n}P(z_i=Ck)+\lambda p_k=0$$
即
$$\begin{array}{c}\begin{array}{rl}& \sum _{k=1}^K\sum _{i=1}^n[P(z_i=Ck)+\sum _{k=1}^K\lambda p_k\\ =& \sum _{i=1}^n\sum _{k=1}^{K}P(z_i=Ck)+\lambda \sum _{k=1}^K{p}_k\\ =& 0\end{array}\end{array}$$
因为${\sum }_{k=1}^K{p}_k=1$，$\sum _{k=1}^{K}P(z_i=Ck)=1$

所以最终变成
$$\sum _{i=1}^{n}1+\lambda =0\to n+\lambda =0\to \lambda =-n$$
将$\lambda =-n$代入之前的$\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)+\lambda p_k=0$得
$$p_k=\frac{1}{n}\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)$$
$p$有了，接下来就是求解$\mu ,\Sigma$。

正态分布的概率密度函数
$$f(x)=\frac{1}{(2\pi {)}^{\frac{d}{2}}|\Sigma {|}^{\frac{1}{2}}}exp\left\{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )\right\}$$
其中$d$代表x的维度。

要求均值和协方差，先把拉格朗日函数里面的正态分布写成概率密度函数的形式
$$\begin{array}{c}\begin{array}{rl}L(\theta ,\lambda )=& \sum _{k=1}^K\sum _{i=1}^n[log{p}_k+logN(x_i|{\mu }_k,{\Sigma }_k)]P(z_i=Ck|x_i,{\theta }^t)+\lambda \left[\sum _{k=1}^K{p}_k-1\right]\\ =& \sum _{k=1}^K\sum _{i=1}^n\left[log{p}_k+log\left[\frac{1}{(2\pi {)}^{\frac{d}{2}}|{\Sigma }_k{|}^{\frac{1}{2}}}exp\left\{-\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k)\right\}\right]\right]P(z_i=Ck|x_i,{\theta }^t)+\lambda \left[\sum _{k=1}^K{p}_k-1\right]\\ =& \sum _{k=1}^K\sum _{i=1}^n\left[log{p}_k+log\frac{1}{(2\pi {)}^{\frac{d}{2}}|{\Sigma }_k{|}^{\frac{1}{2}}}+log\left[exp\left\{-\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k)\right\}\right]\right]P(z_i=Ck|x_i,{\theta }^t)+\lambda \left[\sum _{k=1}^K{p}_k-1\right]\\ =& \sum _{k=1}^K\sum _{i=1}^n\left[log{p}_k+log\frac{1}{(2\pi {)}^{\frac{d}{2}}|{\Sigma }_k{|}^{\frac{1}{2}}}-\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k)\right]P(z_i=Ck|x_i,{\theta }^t)+\lambda \left[\sum _{k=1}^K{p}_k-1\right]\\ =& \sum _{k=1}^K\sum _{i=1}^n\left[log{p}_k-\frac{d}{2}2\pi -\frac{1}{2}log|{\Sigma }_k|-\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k)\right]P(z_i=Ck|x_i,{\theta }^t)+\lambda \left[\sum _{k=1}^K{p}_k-1\right]\end{array}\end{array}$$
对拉格朗日函数关于${\mu }_k$求导，以下直接给出用到得矩阵求导公式
$$\frac{\partial (x^{T}Ax)}{\partial x}=2Ax(假设A为对称阵)$$
矩阵求导依然满足链式求导法则。所以可以将$(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k)$中的$(x_i-{\mu }_k)$当作x，求完外层的导数后再求里面的相乘即可。
$$\begin{array}{c}\begin{array}{rl}\frac{\partial L(\theta ,\lambda )}{\partial {\mu }_k}=& \sum _{i=1}^n{\Sigma }_k^{-1}(x_i-{\mu }_k)P(z_i=Ck|x_i,{\theta }^t)\\ =& {\Sigma }_k^{-1}\sum _{i=1}^n(x_i-u_k)P(z_i=Ck|x_i,{\theta }^t)\\ =& 0\\ 即：& \sum _{i=1}^n(x_i-u_k)P(z_i=Ck|x_i,{\theta }^t)=0\\ =& \sum _{i=1}^n{x}_{i}P(z_i=Ck|x_i,{\theta }^t)-\sum _{i=1}^n{\mu }_{k}P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{i=1}^n{x}_{i}P(z_i=Ck|x_i,{\theta }^t)-{\mu }_k\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)\\ =& 0\end{array}\end{array}$$
移项得
$$u_k=\frac{\sum _{i=1}^n{x}_{i}P(z_i=Ck|x_i,{\theta }^t)}{\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)}$$
再对${\Sigma }_k$求导。

对于${\Sigma }_k$，我们知道，它是一个矩阵，标量对矩阵求导可以对每一个分量求导求解，或者利用迹技巧直接求解。

本文两种都讲一次吧，读者对哪种感兴趣就用哪种

第①种：分量求导

先认识一下下面两个求导常用公式(A为矩阵)，此处不作推导，感兴趣可以百度或者看书
$$\begin{gathered} (In|A|{)}^{{}^{\prime }}=(A^{-1}{)}^T; \\ (A^{-1}{)}^{\prime }=-A^{-1}{A}^{\prime }{A}^{-1}; \end{gathered}$$
对${\Sigma }_k$求导
$$\begin{array}{c}\begin{array}{rl}\frac{\partial L(\theta ,\lambda )}{\partial {\Sigma }_k}=& \sum _{i=1}^n\left[-\frac{1}{2}({\Sigma }_k^{-1}{)}^T+\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}{\Sigma }_k^{{}^{\prime }}{\Sigma }_k^{-1}(x_i-{\mu }_k)\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{i=1}^n\left[-\frac{1}{2}{\Sigma }_k^{-1}+\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}{\Sigma }_k^{{}^{\prime }}{\Sigma }_k^{-1}(x_i-{\mu }_k)\right]P(z_i=Ck|x_i,{\theta }^t)\end{array}\end{array}\text{\hspace{1em}(1)}$$

对于里面的$(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}{\Sigma }_k^{{}^{\prime }}{\Sigma }_k^{-1}(x_i-{\mu }_k)$，我们知道${\Sigma }_k$是协方差矩阵，我们分别对里面的分量进行求导。我们令$A={\Sigma }_k^{-1}(x_i-{\mu }_k)$，则$(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}{\Sigma }_k^{{}^{\prime }}{\Sigma }_k^{-1}(x_i-{\mu }_k)=A^T{\Sigma }_k^{{}^{\prime }}A$，所以
$$\frac{\partial A^T{\Sigma }_{k}A}{\partial {\Sigma }_{ij}}=A_i\ast A_j=(A\ast A^T{)}_{ij}$$
为啥等于这个呢？来看**(以下省略掉${\Sigma }_k^{-1}$,不影响最终结果，后面运算的时候再加回去即可，现在只是证明上面所写的合理性)**
$$A^T{\Sigma }_{k}A=\left(\begin{array}{cc}(x^1-{\mu }^1)& (x^2-{\mu }^2)\end{array}\right)\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right)\left(\begin{array}{c}(x^1-{\mu }^1)\\ (x^2-{\mu }^2)\end{array}\right)$$
再看
$$A{A}^T=\left(\begin{array}{c}(x^1-{\mu }^1)\\ (x^2-{\mu }^2)\end{array}\right)\left(\begin{array}{cc}(x^1-{\mu }^1)& (x^2-{\mu }^2)\end{array}\right)=\left(\begin{array}{cc}(x^1-{\mu }^1)(x^1-{\mu }^1)& (x^1-{\mu }^1)(x^2-{\mu }^2)\\ (x^2-{\mu }^2)(x^1-{\mu }^1)& (x^2-{\mu }^2)(x^2-{\mu }^2)\end{array}\right)$$
对${\Sigma }_{ij}$求导相当于矩阵的每一个元素对${\Sigma }_{ij}$求导，那么理论上也只有对应位置的数值是1，其余为0。因为其余元素被视为常数，而对应位置的求导就是标量对标量的求导，所以直接等于1。比如对${\Sigma }_{11}$求导，所得
$$\begin{array}{c}\begin{array}{rl}\frac{\partial A^T{\Sigma }_{k}A}{\partial {\Sigma }_{11}}=& \left(\begin{array}{cc}(x^1-{\mu }^1)& (x^2-{\mu }^2)\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}(x^1-{\mu }^1)\\ (x^2-{\mu }^2)\end{array}\right)\\ =& (x^1-{\mu }^1)(x^1-{\mu }^1)\\ =& (A{A}^T{)}_{11}\end{array}\end{array}$$
所以，以此类推
$$\frac{\partial A^T{\Sigma }_{k}A}{\partial {\Sigma }_k}=A{A}^T$$
所以公式(1)等于
$$\begin{array}{c}\begin{array}{rl}\frac{\partial L(\theta ,\lambda )}{\partial {\Sigma }_k}=& \sum _{i=1}^n\left[-\frac{1}{2}{\Sigma }_k^{-1}+\frac{1}{2}A{A}^T\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{i=1}^n\left[-\frac{1}{2}{\Sigma }_k^{-1}+\frac{1}{2}{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& -\sum _{i=1}^n\frac{1}{2}{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)+\sum _{i=1}^n\frac{1}{2}{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)\\ =& 0\end{array}\end{array}$$
移项
$$\begin{array}{c}\begin{array}{rlr}& \sum _{i=1}^n\frac{1}{2}{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)=\sum _{i=1}^n\frac{1}{2}{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)& \\ & 即\sum _{i=1}^n{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)=\sum _{i=1}^n{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)& \\ & 等式左右，都左乘以{\Sigma }_k得& \\ & \sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)=\sum _{i=1}^n(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}P(z_i=Ck|x_i,{\theta }^t)& \\ & 等式左右，都右乘以{\Sigma }_k得& \\ & {\Sigma }_k\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)=\sum _{i=1}^n(x_i-{\mu }_k)(x_i-{\mu }_k{)}^{T}P(z_i=Ck|x_i,{\theta }^t)& \end{array}\end{array}$$
所以
$${\Sigma }_k=\frac{\sum _{i=1}^n(x_i-{\mu }_k)(x_i-{\mu }_k{)}^{T}P(z_i=Ck|x_i,{\theta }^t)}{{\sum }_{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)}$$
第②种：迹技巧

对于迹技巧，先来看两个微分公式（A是矩阵，并且可逆），此处不作推导，感兴趣可以百度或者看书
$$\begin{gathered} d|A|=|A|tr(A^{-1}dA); \\ d(A^{-1})=-A^{-1}(dA)A^{-1} \end{gathered}$$
对于迹技巧，对原函数求微分,原函数与${\Sigma }_k$相关的只有两项

第一项$log|{\Sigma }_k|$
$$d(log|{\Sigma }_k|)=\frac{1}{|{\Sigma }_k|}d|\Sigma |=\frac{1}{|{\Sigma }_k|}|{\Sigma }_k|tr({\Sigma }_k^{-1}d{\Sigma }_k)=tr({\Sigma }_k^{-1}d{\Sigma }_k)$$
第二项$(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k)$
$$\begin{array}{c}\begin{array}{r}d((x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(x_i-{\mu }_k))=-(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(d{\Sigma }_k){\Sigma }_k^{-1}(x_i-{\mu }_k)\end{array}\end{array}$$
所以，原函数的微分
$$dL(\theta ,\lambda )=\sum _{i=1}^n\left[-\frac{1}{2}tr({\Sigma }_k^{-1}d{\Sigma }_k)+\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(d{\Sigma }_k){\Sigma }_k^{-1}(x_i-{\mu }_k)\right]P(z_i=Ck|x_i,{\theta }^t)$$
给其套入迹
$$\begin{array}{c}\begin{array}{rl}tr(dL(\theta ,\lambda ))=& tr\left(\sum _{i=1}^n\left[-\frac{1}{2}tr({\Sigma }_k{|}^{-1}d{\Sigma }_k)+\frac{1}{2}(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(d{\Sigma }_k){\Sigma }_k^{-1}(x_i-{\mu }_k)\right]P(z_i=Ck|x_i,{\theta }^t)\right)\\ =& \sum _{i=1}^n\left[\frac{1}{2}tr({\Sigma }_k^{-1}d{\Sigma }_k)+\frac{1}{2}tr((x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}(d{\Sigma }_k){\Sigma }_k^{-1}(x_i-{\mu }_k))\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{i=1}^n\left[-\frac{1}{2}tr({\Sigma }_k^{-1}d{\Sigma }_k)+\frac{1}{2}tr({\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}d{\Sigma }_k)\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{i=1}^n\left[\frac{1}{2}tr(-{\Sigma }_k^{-1}d{\Sigma }_k+{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}d{\Sigma }_k)\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& \sum _{i=1}^n\left[\frac{1}{2}tr\left(-{\Sigma }_k^{-1}+{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1})d{\Sigma }_k\right)\right]P(z_i=Ck|x_i,{\theta }^t)\\ =& tr\left(\sum _{i=1}^n\frac{1}{2}(P(z_i=Ck|x_i,{\theta }^t)(-{\Sigma }_k^{-1}+{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}))d{\Sigma }_k\right)\end{array}\end{array}$$
去掉迹得
$$\begin{array}{c}\begin{array}{rl}& dL(\theta ,\lambda )=\sum _{i=1}^n\frac{1}{2}(P(z_i=Ck|x_i,{\theta }^t)(-{\Sigma }_k^{-1}+{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}))d\Sigma \\ & 即\frac{dL(\theta ,\lambda )}{d{\Sigma }_k}=\sum _{i=1}^n\frac{1}{2}P(z_i=Ck|x_i,{\theta }^t)(-{\Sigma }_k^{-1}+{\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1})=0\end{array}\end{array}$$
移项
$$\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t){\Sigma }_k^{-1}=\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t){\Sigma }_k^{-1}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T{\Sigma }_k^{-1}$$
和上面第一种方法一样，因为$P(z_i=Ck|x_i,{\theta }^t)$是标量，故等式左右，先都左乘${\Sigma }_k$，再都右乘${\Sigma }_k$，得
$${\Sigma }_k\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)=\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T$$
最终
$${\Sigma }_k=\frac{\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T}{\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)}$$

结果

$$\begin{gathered} p_k=\frac{1}{n}\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t); \\ {u}_k=\frac{\sum _{i=1}^n{x}_{i}P(z_i=Ck|x_i,{\theta }^t)}{\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)}; \\ {\Sigma }_k=\frac{\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)(x_i-{\mu }_k)(x_i-{\mu }_k{)}^T}{\sum _{i=1}^{n}P(z_i=Ck|x_i,{\theta }^t)} \end{gathered}$$

那么，接下来只需要计算出$P(z_i=Ck|x_i,{\theta }^t)$，$\theta$是参数，下面省略掉
$$P(z_i=Ck|x_i)=\frac{P(z_i=Ck,x_i)}{P(x_i)}$$
对于$P(x_i)$
$$P(x_i)=\sum _{z_i}P(x_i,z_i)=\sum _{k=1}^{K}P(x_i,z_i=Ck)$$
前面我们算出来过$P(x_i,z_i=Ck)=p_k\ast N(x_i|{\mu }_k,{\Sigma }_k)$

所以
$$P(z_i=Ck|x_i)=\frac{p_k\ast N(x_i|{\mu }_k,{\Sigma }_k)}{\sum _{k=1}^K{p}_k\ast N(x_i|{\mu }_k,{\Sigma }_k)}$$
请务必注意式子中，分子处得k来自左边的Ck，而分母的k是来自求和符号

算法流程

①随机初始化模型参数$p^t,{\mu }^t,{\Sigma }^t$。

②计算出$P(z_i=Ck|x_i,\theta )$

③依据公式计算出$p^{t+1},{\mu }^{t+1},{\Sigma }^{t+1}$

④计算$p^{t+1},{\mu }^{t+1},{\Sigma }^{t+1}$和$p^t,{\mu }^t,{\Sigma }^t$的差值，如果差值小于$ϵ$(自己设定的值)。如果小于则说明变化太小，证明收敛，结束算法。否则循环②，③步骤

代码实现

GMM高斯混合模型代码实现

结束

至此推导和代码已都全部完成。很多地方推导并不严谨，如有问题，还请指出。阿里嘎多