从 GMM 到 EM 算法

首先需要声明的是，GMM是Gaussian Mixture Model，混合高斯模型，是一个模型。EM算法，Expection Maximization期望最大是一套计算框架（framework），一系列的计算流程（迭代式的）。一般地，我们可以使用EM算法来求解GMM问题，即 put GMM into EM framework。

EM迭代流程为：

$$\theta^{(g+1)}=\underbrace{\arg\max_\theta}_{\textrm{M-Step}}\underbrace{\int_z\log p(X,z|\theta)p(z|X,\theta^{(g)})dz}_{\textrm{E-Step}}$$

$X=\{x_1, x_2,\dots,x_N\}$
其 joint density 或者叫 log joint density：

$$\log P(X)=\sum_{i=1}^N\log\sum_{\ell=1}^k\alpha_k\mathcal{N}(x_i|\mu_\ell,\sigma_\ell)$$

混合高斯模型（GMM）参数集 $\Theta=\{\mu_1, \dots,\mu_k,\Sigma_1,\dots,\Sigma_k,\alpha_1,\dots,\alpha_{k-1} \}$

$$P(x_i|\Theta)=\sum_{\ell=1}^k\alpha_\ell\mathcal{N}(x_i|\mu_\ell,\Sigma_\ell)$$

$$\Theta_{MLE}=\arg\max_\Theta\sum_{i=1}^N\log\sum_{\ell=1}^k\alpha_\ell\mathcal{N}(x_i|\mu_\ell,\Sigma_\ell)$$

其中 $\mathcal{L}(\Theta|X)=\sum\limits_{i=1}^N\log\sum\limits_{\ell=1}^k\alpha_{\ell}\mathcal{N}(x_i|\mu_{\ell},\Sigma_{\ell})$，混合的情形下，不再像单高斯（single mode Gaussian）的情况，不存在一个显示的解析解，一种常规的替代方案是使用迭代的方式去寻找（像MCMC算法那样，收敛到稳态？:-D），而这一方式正是著名的EM算法。

所谓迭代的方式即是提供如下的一种递归关系：

$$\Theta^{(g+1)}=f(\Theta^{(g)})$$

EM算法给出的 $f(\cdot)$ 是：

$$\Theta^{(g+1)}=\arg\max_{\Theta}\int_zP(X,z|\Theta)P(z|X,\Theta^{(g)})dz$$

该递推关系还应至少满足，$\log P(X|\Theta^{(g+1)})>\log P(X|\Theta^{(g)})$（Log Likelihood）

其中引变量 $z_i$（$z_i=\{1,\dots,k\}$） 标识样本 $x_i$ 所属的高斯号（哪一个高斯），这样就将GMM的fitting问题转换为了single mode Gauss的fitting问题了（将每个高斯对应的数据摘出来）。引变量的存在使得问题得以简化。对所添加的隐藏变量的要求是不能改变边缘分布（marginal distribution）：

$$p(x_i)=\int_{z_i}\underbrace{p(x_i|z_i)}_{\mathcal{N}(x_i|\mu_{z_i},\Sigma_{z_i})}\underbrace{p(z_i)}_{\alpha_{z_i}}dz_i$$

又因为 $z_i$是离散型随机变量，取值为 $z_i=\{1,\dots,k\}$，又可将积分符号改造为求和符号。也即：

$$p(x_i)=\sum_{z_i=1}^kp(x_i|z_i)p(z_i)=\sum_{z_i=1}^k \alpha_{z_i}\mathcal{N}(x_i|\mu_{z_i},\Sigma_{z_i})$$

此时关于 $x_i$的边缘分布，刚好就是 $\mathcal{L}(\Theta|X)=\sum\limits_{i=1}^N\log\underbrace{\sum\limits_{\ell=1}^k\alpha_\ell\mathcal{N}(x_i|\mu_{\ell},\Sigma_{\ell})}$

也即添加引变量之后并未改变数据的边缘分布。

$$\log p(X|\Theta)=\log p(X,z|\Theta)-\log p(z|X,\Theta)\\ \Downarrow（大名鼎鼎的EM期望最大算法）\\ E_{p(z|X,\Theta^{(g)})}[\log p(X|\Theta)]=E_{p(z|X,\Theta^{(g)})}[\log p(X,z|\Theta)]-E_{p(z|X,\Theta^{(g)})}[\log p(z|X,\Theta)]\\ \Downarrow\\ \log p(X|\Theta)=\underbrace{\int_z \log p(X,z|\Theta)p(z|X,\Theta^{(g)})dz}_{Q(\Theta,\Theta^{(g)})}-\underbrace{\int_z\log p(z|X,\Theta)p(z|X,\Theta)dz}_{H(\Theta,\Theta^{(g)})}\\ \Downarrow\\ \log p(X|\Theta)=Q(\Theta, \Theta^{(g)})-H(\Theta, \Theta^{(g)})\\ $$

put GMM into EM framework

EM framework：

$$\Theta^{(g+1)}=\arg\max_\Theta\int_z\log p(X,Z|\Theta)p(Z|X,\Theta^{(g)})dz$$

• 如何定义 $p(X,Z|\Theta)$（$z_i$ 是对高斯的指定）：

$$p(X,Z|\Theta)=\prod_{i=1}^Np(x_i,z_i|\Theta)=\prod_{i=1}^N\underbrace{p(x_i|z_i,\Theta)}_{\mathcal{N}(\mu_{z_i},\Sigma_{z_i})}\underbrace{p(z_i|\Theta)}_{\alpha_{z_i}}=\prod_{i=1}^N\alpha_{z_i}\mathcal{N}(\mu_{z_i},\Sigma_{z_i})$$

还记得 $p(X|\Theta)$ 的形式吗？

$$p(X|\Theta)=\sum_{\ell=1}^k\alpha_\ell\mathcal{N}(X|\mu_\ell,\sigma_\ell)=\prod_{i=1}^N\sum_{\ell=1}^k\alpha_\ell\mathcal{N}(x_i|\mu_\ell,\sigma_\ell)$$
可见 $p(X,Z|\Theta)$ 是 $p(X|\Theta)$ 关于高斯的指定。

• 再来看 $p(Z|X,\Theta)$ 的定义

因为 $(x_i,z_i)$彼此独立：

$$p(Z|X,\Theta)=\prod_{i=1}^Np(z_i|x_i,\Theta)=\prod_{i=1}^N\frac{p(z_i)p(x_i|z_i)}{\sum_{z_i}p(z_i)p(x_i|z_i)}=\prod\frac{\alpha_{z_i}\mathcal{N}(x_i|\mu_{z_i},\Sigma_{z_i})}{\sum \alpha_{z_i}\mathcal{N}(x_i|\mu_{z_i},\Sigma_{z_i})}$$

代入到EM的框架下：

E-Step:

$$=\sum_{i=1}^N\sum_ {z_i=1}^k\left ( \log \alpha_{\ell}+\log \mathcal{N}(x_i|\mu_{\ell},\Sigma_{\ell})\right )p(\ell|x_i,\Theta^{(g)})\\ =\sum_{i=1}^N\sum_ {\ell=1}^k\left ( \log \alpha_{\ell}+\log \mathcal{N}(x_i|\mu_{\ell},\Sigma_{\ell})\right )p(\ell|x_i,\Theta^{(g)})$$
M-Step:
$$\Rightarrow \alpha_\ell=\frac1N\sum_{i=1}^Np(\ell|x_i,\Theta^{(g)})$$

补充

GMM问题求解的困难在于，$\mathcal{L}(\Theta|X)=\sum_{i=1}^N\log\sum_{\ell=1}^k\alpha_\ell\mathcal{N}(x_i|\mu_\ell,\Sigma_\ell)$，对和式求对数。（numpy提供了np.logaddexp()，:-D）

>>> import numpy as np
>>> np.logaddexp(np.log(1), np.log(1))
0.69314718055994529
                            # np.logaddexp(x, y) = log(e^x+e^y)
>>> np.log(1+1)
0.69314718055994529
                            # np.logaddexp(log(x), log(y)) = log(x+y)