【机器学习算法推导】高斯混合模型GMM与EM算法

1.高斯混合模型

  极大似然估计是一种应用很广泛的参数估计方法。在已有某个地区身高数据以及知道身高服从高斯分布的情况下，利用极大似然估计的方法可以估计出高斯分布$\mu ,\sigma$两个参数。
  如果是多组数据，多个模型呢？获取现在我们有全国多个省份的身高数据，但并不知道它们具体属于哪个省份，只知道每个省之间服从不同的高斯分布，此时的模型称为高斯混合模型（GMM），其公式为
$$P(x;\theta )=\sum _{i=1}^k{p}_{k}N(x;{\mu }_k,{\Sigma }_k)..........(1)$$
  此时用极大似然估计的方法并不能很好地求出模型里面的参数，因为模型中存在一个隐变量——样本的label。令其为Z，则可以写出似然函数的形式：
$$L(\theta )=\sum _{i=1}^m\log p(x_i;\theta )=\sum _{i=1}^m\log \int p(x_i,z;\theta )dz..........(2)$$
  如果根据极大似然估计的算法对式子(2)进行求导，意味着我们要对$\log (y1+y2+...)$的形式进行求导，这计算量简直爆炸。

2.Jensen不等式

  这里我们需要用到EM算法来进行求解，在介绍EM算法前，我们先来了解一个不等式——Jensen不等式，它的定义如下：
  如果X是一个随机变量，f(X)是一个凸函数（二阶导数大或等于0），那么有：
$$E(f(x))\ge f[E(x)]$$
  如果f（X）是凹函数，不等号反向
$$E(f(x))\le f[E(x)]$$
  当且仅当X是常数的时候等号成立。

  结合图像我们可以对Jensen不等式有一个直观的理解，对于x=a和x=b两个点，E(f(x))可以看成是$\frac{1}{2}(f(a)+f(b))$，f[E(x)]可以看成是$f(\frac{1}{2}(f(a)+f(b)))$，即ab中点对应的函数值。由于图像是凹图像，所以不难得出$E(f(x))\ge f[E(x)]$。

3.EM算法及推导过程

  利用Jensen不等式，由于log函数是凹函数，式子(2)可以改写成：
$$\begin{array}{rl}\sum _{i=1}^m\log \int p(x_i,z;\theta )dz& =\sum _{i=1}^m\log \int Q(z)\frac{P(x_i,z;\theta )}{Q(z)}dz\\ & \ge \sum _{i=1}^m\int Q(z)\log \frac{P(x_i,z;\theta )}{Q(z)}dz..........(3)\end{array}$$
  由于$Q(z)$是概率分布，所以满足
$$\sum _{z}Q(z)=1..........(4)$$
  根据Jenson不等式，等式成立时，有
$$\frac{p(x_i,z;\theta )}{Q(z)}=C..........(5)$$
  根据4和5，有
$$\begin{gathered} \sum _{z_i}Q(z)=\sum _{z_i}\frac{p(x_i,z;\theta )}{C}=1..........(6) \\ Q(z)=\frac{p(x_i,z;\theta )}{C}=\frac{p(x_i,z;\theta )}{{\sum }_{z}p(x_i,z;\theta )}=\frac{p(x_i,z;\theta )}{p(x_i)}=p(z|x_i;\theta )..........(7) \end{gathered}$$
  将(7)回代到(3)，有
$$\sum _{i=1}^m\log \int p(x_i,z;\theta )dz\ge \sum _{i=1}^m\int p(z|x_i;\theta )\log \frac{P(x_i,z;\theta )}{p(z|x_i;\theta )}dz$$
  通过每次令
$$\hat{\theta }=\arg \underset{\theta }{\max }\int p(z|x_i;\theta )\log \frac{P(x_i,z;\theta )}{p(z|x_i;\theta )}dz=\int p(z|x_i;\theta )\log P(x_i,z;\theta )dz$$
  就可以不断迭代使得似然函数最大。因此EM算法本质上是一种令下界不断提高，从而使得似然函数也不断提高的迭代算法。
  因此，EM算法的流程如下：

1. 初始化参数${\theta }_0$，开始迭代
2. E-Step：计算$$Q(\theta ,{\theta }_i)=E_{z|x,{\theta }_i}[\log P(x,z|\theta )]=\int \log P(x,z|\theta )P(z|x,{\theta }_i)dz$$
3. M-Step：根据计算得到的$Q(\theta ,{\theta }_i)$，求使$Q(\theta ,{\theta }_i)$极大化的$\theta$,得到新的参数${\theta }_{i+1}$。
4. 重复2和3直至收敛

4.EM算法的可行性

  为什么EM算法能近似实现对观测数据的极大似然估计呢？我们需要证明我们的对数似然函数的值在迭代的过程中一直在增大。假设当前迭代到第i轮，参数为${\theta }_i$，希望新估计值$\theta$可以使$L(\theta )$增大，即$L(\theta )>L({\theta }_i)$，则：　
$$\begin{array}{rl}L(\theta )-L({\theta }_i)& =\log (\sum _{Z}p(x|z,\theta )p(z|\theta ))-\log \left(p(x|{\theta }_i)\right)\\ & =\log \sum _{Z}p(z|x,{\theta }_i)\frac{p(x|z,\theta )p(z|\theta )}{p(z|x,{\theta }_i)}-\log \left(p(x|{\theta }_i)\right)\\ & \ge \sum _{Z}p(z|x,{\theta }_i)\log \frac{p(x|z,\theta )p(z|\theta )}{p(z|x,{\theta }_i)}-\log \left(p(x|{\theta }_i)\right)\\ & =\sum _{Z}p(z|x,{\theta }_i)\log \frac{p(x|z,\theta )p(z|\theta )}{p(z|x,{\theta }_i)(p(x|{\theta }_i)}...........(8)\end{array}$$
  令$$B(\theta ,{\theta }_i)=L({\theta }_i)+\sum _{Z}p(z|x,{\theta }_i)\log \frac{p(x|z,\theta )p(z|\theta )}{p(z|x,{\theta }_i)(p(x|{\theta }_i)}...........(9)$$
  则有$L(\theta )\ge B(\theta ,{\theta }_i)$，因此可以认为$B(\theta ,{\theta }_i)$为$L(\theta )$的一个下界，且有$L({\theta }_i)=B({\theta }_i,{\theta }_i)$，为了使$L(\theta )$尽可能地变大，我们解以下式子：
$$\begin{array}{rl}\arg \underset{\theta }{\max }B(\theta ,{\theta }_i)& =\arg \underset{\theta }{\max }L({\theta }_i)+\sum _{Z}p(z|x,{\theta }_i)\log \frac{p(x|z,\theta )p(z|\theta )}{p(z|x,{\theta }_i)(p(x|{\theta }_i)}\\ & =\arg \underset{\theta }{\max }\sum _{Z}p(z|x,{\theta }_i)\log (p(x|z,\theta )p(z|\theta ))\\ & =\arg \underset{\theta }{\max }\sum _{Z}p(z|x,{\theta }_i)\log p(x,z|\theta )\\ & =\arg \underset{\theta }{\max }Q(\theta ,{\theta }_i)...........(10)\end{array}$$
  而这正是EM中我们所做的一次迭代，由此可证EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

5.EM算法的收敛性

  可以使用单调有界的思想来证明EM算法收敛。首先，$P(x|\theta )$为观测数据的似然函数，我们要证明递增性，即
$$P(x|{\theta }_{i+1})\ge P(x|{\theta }_i)...........(11)$$
  证明如下:
$$\begin{array}{rl}& P(x|\theta )=\frac{P(x,z|\theta )}{P(z|x,\theta )}⟺\log P(x|\theta )=\log P(x,z|\theta )-\log P(z|x,\theta )...........(12)\end{array}$$
  令
$$H(\theta ,{\theta }_i)=\sum _z\log P(z|x,\theta )P(z|x,{\theta }_i)...........(13)$$
  根据Q函数定义有，
$$Q(\theta ,{\theta }_i)=\sum _z\log P(x,z|\theta )P(z|x,{\theta }_i)$$
  于是对数似然函数可以写成
$$\log P(x|\theta )=Q(\theta ,{\theta }_i)-H(\theta ,{\theta }_i)..........(14)$$
  上式中分别取${\theta }_i$ 和${\theta }_{i+1}$并相减，有
$$\begin{array}{rl}& \log P(x|{\theta }_{i+1})-\log P(x|{\theta }_i)\\ & =[Q({\theta }_{i+1},{\theta }_i)-Q({\theta }_i,{\theta }_i)]-[H({\theta }_{i+1},{\theta }_i)-H({\theta }_i,{\theta }_i)]..........(15)\end{array}$$
  由于${\theta }_{i+1}$使$Q(\theta ,{\theta }_i)$极大，所以有
$$Q({\theta }_{i+1},{\theta }_i)-Q({\theta }_i,{\theta }_i)\ge 0..........(16)$$
  而由式子13以及Jensen不等式，有
$$\begin{array}{rl}& H({\theta }_{i+1},{\theta }_i)-H({\theta }_i,{\theta }_i)\\ & =\sum _z\left(\log \frac{P(z|x,{\theta }_{i+1})}{P(z|x,{\theta }_i)}P(z|x,{\theta }_i)\right)\\ & \le \log \left(\sum _z\frac{P(z|x,{\theta }_{i+1})}{P(z|x,{\theta }_i)}P(z|x,{\theta }_i)\right)\\ & =\log P(z|x,{\theta }_{i+1})=0.........(17)\end{array}$$
  可知
$$\log P(x|{\theta }_{i+1})-\log P(x|{\theta }_i)\ge 0..........(18)$$
  又因为$P(x|{\theta }_i)\le 1$有界，可知$L({\theta }_i)=\log P(x|{\theta }_i)$必定收敛。
  由于采用迭代的方式进行求解，EM算法可以保证收敛到一个稳定点，但是却不能保证收敛到全局的极大值点，因此它是局部最优的算法，当然，如果我们的优化目标$B(\theta ,{\theta }_i)$是凸的，则EM算法可以保证收敛到全局最大值，这点和梯度下降法这样的迭代算法相同。初值的选择非常重要，常用的方法是选择多个初值尝试，比较结果选取较好的。

6.EM的另一种推导

$$\begin{array}{rl}\log P(x|\theta )& =\log P(x,z|\theta )-\log P(z|x,\theta )\\ & =\log \frac{P(x,z|\theta )}{q(z)}-\log \frac{P(z|x,\theta )}{q(z)}\end{array}$$
  左右两边同时对z积分，此时有
$$left=\int \log P(x|\theta )q(z)dz=\log P(x|\theta )\int q(z)dz=P(x|\theta )$$
  右边有
$$right={\int }_{z}q(z)\log \frac{P(x,z|\theta )}{q(z)}dz-{\int }_{z}q(z)\log \frac{P(z|x,\theta )}{q(z)}dz$$
  根据相对熵（KL散度）的定义，对于p(x)和q(x)两个分布，p对q的相对熵表示为
$$D(p||q)=\sum _{i=1}^{n}p(x)\log \frac{p(x)}{q(x)}$$
  可以发现right等式右边第二项其实就是q(z)对$P(z|x,\theta )$的相对熵，即
$$P(x|\theta )={\int }_{z}q(z)\log \frac{P(x,z|\theta )}{q(z)}dz+D(q(z)||p(z|x,\theta ))$$
  相对熵有两个性质，一个是不具有对称性，即$D(p||q)̸=D(q||p)$，另一个性质是$D(p||q)\ge 0$。当p和q的概率分布相同时，相对熵为0。也就是说，第一项其实就是$p(x|\theta )$的一个下界，我们称之为ELBO。我们想要让下界ELBO增大，从而使得$p(x|\theta )$增大，则
$$\begin{array}{rl}\hat{\theta }& =\arg \underset{\theta }{\max }{\int }_{z}q(z)\log \frac{P(x,z|\theta )}{q(z)}dz\\ & =\arg \underset{\theta }{\max }{\int }_{z}P(z|x,{\theta }_i)\log \frac{P(x,z|\theta )}{P(z|x,{\theta }_i)}dz\\ & =\arg \underset{\theta }{\max }{\int }_{z}P(z|x,{\theta }_i)\log P(x,z|\theta )dz\end{array}$$
  因此得到$\theta$的迭代公式。

7.应用EM算法求解GMM

  下面我们使用EM算法来求解GMM模型，GMM模型如下：
$$p(x|\theta )=\sum _{k=1}^K{p}_{k}N(x|{\mu }_k,{\Sigma }_k).........(19)$$
  联合概率为
$$p(x,z)=p(z)p(x|z)=p_k\cdot N(x|{\mu }_z,{\Sigma }_z)$$
  后验概率为
$$p(z|x)=\frac{p(x,z)}{p(x)}=\frac{p_k\cdot N(x|{\mu }_z,{\Sigma }_z)}{{\sum }_{k=1}^K{p}_k\cdot N(x|{\mu }_k,{\Sigma }_k)}$$
  在EM算法中，有
$$\begin{array}{rl}Q(\theta ,{\theta }_t)& ={\int }_z\log P(x,z|\theta )P(z|x,{\theta }_t)dz\\ & =\sum _z\log \prod _{i=1}^{N}P(x_i,z_i|\theta )\prod _{i=1}^{N}P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_1,z_2,...,z_n}\sum _{i=1}^N\log P(x_i,z_i|\theta )\prod _{i=1}^{N}P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_1,z_2,...,z_n}(\log P(x_1,z_1|\theta )+\log P(x_2,z_2|\theta )+...+\log P(x_n,z_n|\theta ))\prod _{i=1}^{N}P(z_i|x_i,{\theta }_t)\end{array}$$
  对于其中一项
$$\begin{array}{rl}& \sum _{z_1,z_2,...,z_n}\log P(x_1,z_1|\theta )\prod _{i=1}^{N}P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_1,z_2,...,z_n}\log P(x_1,z_1|\theta )P(z_1|x_1,{\theta }_t)\prod _{i=2}^{N}P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_1}\log P(x_1,z_1|\theta )P(z_1|x_1,{\theta }_t)\sum _{z_2,...,z_n}\prod _{i=2}^{N}P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_1}\log P(x_1,z_1|\theta )P(z_1|x_1,{\theta }_t)\sum _{z_2}P(z_2|x_2,{\theta }_t)\sum _{z_3}P(z_3|x_3,{\theta }_t)...\sum _{z_n}P(z_3|x_n,{\theta }_t)\\ & =\sum _{z_1}\log P(x_1,z_1|\theta )P(z_1|x_1,{\theta }_t)\end{array}$$
  因此
$$\begin{array}{rl}Q(\theta ,{\theta }_t)& =\sum _{z_1,z_2,...,z_n}(\log P(x_1,z_1|\theta )+\log P(x_2,z_2|\theta )+...+\log P(x_n,z_n|\theta ))\prod _{i=1}^{N}P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_1}\log P(x_1,z_1|\theta )P(z_1|x_1,{\theta }_t)+...+\sum _{z_n}\log P(x_n,z_n|\theta )P(z_n|x_n,{\theta }_t)\\ & =\sum _{i=1}^N\sum _{z_i}\log P(x_i,z_i|\theta )P(z_i|x_i,{\theta }_t)\\ & =\sum _{i=1}^N\sum _{z_i}\log (p_{z_i}\cdot N(x_i|{\mu }_{z_i},{\Sigma }_{z_i}))\cdot P(z_i|x_i,{\theta }_t)\\ & =\sum _{z_i}\sum _{i=1}^N\log (p_{z_i}\cdot N(x_i|{\mu }_{z_i},{\Sigma }_{z_i}))\cdot P(z_i|x_i,{\theta }_t)\\ & =\sum _{k=1}^K\sum _{i=1}^N\log (p_k\cdot N(x_i|{\mu }_k,{\Sigma }_k))\cdot P(z_i=c_k|x_i,{\theta }_t)\\ & =\sum _{k=1}^K\sum _{i=1}^N(\log p_k+\log N(x_i|{\mu }_k,{\Sigma }_k))\cdot P(z_i=c_k|x_i,{\theta }_t)\end{array}$$
  接下来求解${\theta }_{t+1}=(p_k,{\mu }_k,{\Sigma }_k)$
$$\begin{gathered} p_k^{t+1}=\arg \underset{p_k}{\max }\sum _{k=1}^K\sum _{i=1}^N\log p_k\cdot P(z_i=c_k|x_i,{\theta }_t),s.t.\sum _{k=1}^K{p}_k=1 \\ {\mu }_k^{t+1}=\arg \underset{{\mu }_k}{\max }\sum _{k=1}^K\sum _{i=1}^N\log N(x_i|{\mu }_k,{\Sigma }_k)\cdot P(z_i=c_k|x_i,{\theta }_t) \\ {\Sigma }_k^{t+1}=\arg \underset{{\Sigma }_k}{\max }\sum _{k=1}^K\sum _{i=1}^N\log N(x_i|{\mu }_k,{\Sigma }_k)\cdot P(z_i=c_k|x_i,{\theta }_t) \end{gathered}$$
  通过展开之后求导得到最优解，其中$p_k^{t+1}$的求解使用拉格朗日乘子法进行求解。重复以上E-step和M-step直至对数似然值不再有明显变化为止。