EM基础

EM算法主要是求解隐变量模型，其中最基础的一类方法叫做混合模型，最典型的就是GMM（高斯混合模型）。隐变量模型中，观测数据是基于隐变量为基础生成的结果，所以也算是生成模型。

通常来说，对混合模型的定义是，我们有一个分布 $p(z)$ ，用来从一批模型中挑选一个，然后 $p(x|z)$ 代表选取了 $z$ 之后对应的那个模型。比如GMM中 $p(x|z)$ 就是高斯分布。我们实际拿到的观测数据 $x$ 就是从分布 $p(z)$ 先选定一个模型，再基于 $p(x|z)$ 采样出来的结果。

接下来主要引用基于PRML的叙述来解释。基于观测数据我们希望解出GMM参数，那么要从计算似然函数出发。由于无法观测隐变量，似然函数相当于把所有隐变量可能的取值都过一遍做加和：

这里 $\theta$ 是所有的模型参数，包括 $p(z|\theta)$ （选哪个模型时的权重分布）以及 $p(x|z, \theta)$ （每个模型各自的参数），因此有 $\sum_{z}p(x,z| \theta) = \sum_{z}p(z| \theta)p(x|z, \theta)$ 。

这样一个似然函数，对于GMM这样的简单模型来说，用基于梯度的方法优化也是可行的（PRML 9.2末尾，以及这个材料都提到过），但是用EM解会更加高效。

为了理解EM，我们需要给隐变量引入一个不确定的分布 $q(z)$ 来把似然函数拆解成两部分：

前面一项为ELBO，后面一项为KL：

基于上面的分拆，我们就可以使用EM了。EM框架也是一个迭代逻辑，每一轮迭代基于当前的参数 $\theta^{(old)}$ 算出一个让似然函数更大的 $\theta^{(new)}$ 。每一轮迭代具体步骤分为E跟M两步：

对于E步，其核心目的是要构造一个最逼近似然函数的ELBO，这个我们可以通过让KL这项为0来完成。为了让KL为0，我们就需要令 $q(z) = p(z|x, \theta^{(old)})$ ，也就是 $z$ 关于 $x$ 跟参数 $\theta^{(old)}$ 的后验概率分布。所以这步的主要工作量就是算出 $p(z|x, \theta^{(old)})$ 。
对于M步，其核心目的就是让ELBO最大，且ELBO中的 $q(z)$ 已经在E步给出。我们把在E步算出来的结果带入到ELBO就有下面的结果，其中有一部分是常数项，可以忽略，真正要优化的是 $Q(\theta, \theta^{(old)})$ 。也就是我们要算 $\theta^{(new)} = argmax_\theta Q(\theta, \theta^{(old)})$ 。

以上就是EM基本的迭代过程，下图是一个可视化的解释。对于E步来说，构造了一个ELBO（蓝线），因为在 $\theta^{(old)}$ 处KL为0，所以这时的ELBO跟似然在此处相等，其他地方则是小于似然（红线），因为其他地方KL>0。然后在M步我们取到一个让本轮ELBO最大的位置 $\theta^{(new)}$ ，那对应位置的似然值必然更大。下一轮迭代基于 $\theta^{(new)}$ 重新构造新一轮的ELBO（绿线），以此类推，直到逼近最大的似然值。

我们用GMM作为具体的例子做一下计算。我们假设GMM有 $K$ 个高斯模型混合而成，第 $k$ 个高斯的参数为 $\mu_k, \Sigma_k$ ，并且第k个高斯在选取模型时候被选到的概率为 $\pi_k$ ，也是参数的一部分：

然后我们先从E步开始，给定每一条观测样本 $x$ ，计算对应后验概率分布，我们记为 $\gamma(z_k)$ ，也就是根据这条样本，判断他是从第 $k$ 个模型生成的概率：

这里利用 $p(z|x) = p(x|z)p(z)/\sum_zp(x,z)$ 展开得到结果，这里书上简化了参数没写上标old，但是我们确切知道他们的当前值。

然后对于M步，我们写出ELBO，让其对 $\mu_k, \Sigma_k, \pi_k$ 分别求导为0即可得到ELBO的argmax。比如对于 $\mu_k$ 求导：

就有：

对于 $\Sigma_k, \pi_k$ 的结果这里就不在赘述。

EM与变分推断

上述最基本的EM算法只能应对GMM等简单的隐变量模型。对于复杂一点的模型来说，我们很难精确的计算后验概率 $p(z|x, \theta^{(old)})$ ，这时候我们可以用一个带参数的分布 $q(z|\phi)$ 去拟合 $p(z|x, \theta^{(old)})$ 。这里引入了一组新的变分参数 $\phi$ 。这样我们在E步算出 $q(z|\phi)$ 之后，在M步仍然计算 $Q$ 函数，只不过此时由于我们用的是变分推断拟合的结果， $Q$ 函数可以写作 $Q(\theta, \phi)$ ，这个就是所谓变分EM。

变分推断的核心就是用 $q$ 拟合 $p$ 使得 $KL(q||p)$ 最小，基于我们前面推导的关于似然/KL/ELBO的三者关系，我们知道 $\theta$ 确定的时候似然值是确定的，那么求 $KL(q||p)$ 等价于求ELBO最大，这也是实际上做变分推断的核心思路。另外由于我们是用一个新的分布 $q$ 去近似，那么我们可以让这个分布尽可能的简单，其中核心要点就是让各个维度互不相关，也就是平均场理论。比如我们的 $q$ 选一个高斯分布，那么其协方差矩阵应该是个对角阵，任意两个维度之间协方差应为0。

变分EM不像原始EM能做到每一步似然都一定上涨最后到达最大值，因为每次的KL距离并不能保证等于0，所以存在KL距离会有波动的情况，也就是变分EM可能像下左图一样，但是也可能像下右图那样。我们能保证的就是ELBO一直在涨，并且似然是大于当前这轮的ELBO。因此，可以把KL看成是某种（作用在似然函数上的）正则项，是一种为了让变分推断容易实现而选取简单的 $q$ 带来的代价。

为了给变分EM加快速度，我们在每一轮迭代的M步计算 $Q(\theta, \phi)$ 的时候不要用所有的数据，而是随机抽取一个batch的数据做近似，这样我们在E步也就不用将所有数据点对应的 $q(z|\phi)$ 都算出来了。这种类似从GD演化到SGD的思路就叫做SVI。

但是以上方法中，我们的推断都是对每一条样本数据点都设置了一个对应的分布，但是实际上我们可以在整个数据集上训练一个分布 $q(z|x, \phi)$ ，输入一个数据 $x$ 输出一个对应分布，这样我们维护的参数量可以大大减少，进一步加快推断速度。这种思路就是Amortized VI。VAE就是用了这种思想，使用一个encoder网络做x到φ的映射。由于VAE这个工作的流行，Amortized VI可以认为是接受度较高的一种VI方法，这个在另一篇文章中有进一步的展开。

奔跑的兔子：VAE(变分自编码器)快速复习

MCEM

MCEM采取另一种估计 $p(z|x)$ 的方法，是直接从 $p(z|x)$ 中采样作为 $q$ 的代替，这样在M步使用Monte Carlo计算ELBO。由于这个方案完全依赖于采样也不做近似，所以碰到复杂的后验分布，需要用到MCMC采样。变分EM系列的方法，当然也有可能碰到M步需要Monte Carlo的情况，但由于 $q$ 是一个较为简单的近似分布，所以至少不需要用MCMC这么复杂收敛慢的方案。而且MCEM基本上对应了早期变分EM的缺点，就是需要每个数据点独立处理，也就是对每个数据点对应的后验分布做的采样，而变分这边已经演化到了AVI，明显快很多。因此这类方法目前看实用性较低了。