EM算法通用形式（ESL 8.5.2）

8.5.2的EM算法的描述比较晦涩，这里总结一下EM算法的general形式。EM算法主要用于最大化似然函数。某些似然函数要最大化可能比较困难，但如果可以引入某些隐变量（latent data），那么这个最大化问题可能会变得简单一些。在The Element of Statistical Learning这本书里面这叫做data augmentation。需要注意的是，一般机器学习里面（尤其深度学习）说的data augmentation指的是对原来的数据进行一定的变换（比如计算机视觉里面把图片旋转拉伸放缩）再加入到原始数据里面，这样数据量大了，一定程度上面可以减少泛化误差，或者增加训练样本的量。可以参见巨著Deep Learning里面7.4节。还要注意，EM算法本质上是一个优化算法。从机器学习算法的三大构成（模型，目标函数，优化方法）来看，假定的数据分布是模型，最大似然是目标函数，而EM算法是用以优化出结果的方法。用形式化的方法来说，如果假定的数据分布的具有参数$\theta$，拟合的数据是$x_i$，那么似然函数则为

$$\sum_{i} log\ p(x_i;\theta).$$
很多人可能对于那个分号有点懵逼，其实就是表示给定参数 $\theta$下 $x_i$的概率。

8.5.2对EM算法通用形式的描述有点不太好懂，先从吴恩达的比较好懂的版本开始说起。理解EM算法需要一个数据工具，Jensen不等式。

Theorem 如果$f$是一个凸函数，$X$是一个随机变量，那么

$$E[f(X)] \geq f(EX)$$
如果 $f$是严格凸的，那么仅当$X=EX$的时候上式取等号。

证明这个不等式也不难，需要用到凸函数的性质。如果$x,c$在$f$的定义域上，那么

$$f(x)-f(c)\geq f'(c)(x-c)。$$
即$f(x)-f(c)-f'(c)(x-c)\geq 0$。假设$c=EX$，那么等式两边对$X$的密度函数$g(x)$积分（如果是离散变量就是乘以对应的概率函数求和了，或者说对概率测度求积分），可以得到

$$\begin{equation} \begin{aligned} \int f(x)g(x)dx-f(c)-f'(c)(\int xg(x)dx -c ) \geq 0 \\ E[f(x)]-f(c)\geq 0 \end{aligned} \end{equation}$$
因为 $c=EX$，第三项消掉。继续代入，就可以证明出Jensen不等式了。因为凸函数和凹函数是对称的，因此如果 $f$是一个凹函数，那么$E[f(x)] \leq f(EX)$

回到最大似然的问题，假设我们额外增加的latant variable隐变量是$z$，对于混合高斯模型来说，这个隐变量就是指定某个$x_i$属于哪一个高斯分布。对于不同的$x_i$，$z$会有不同的分布。我们记对于$x_i$来说$z$的分布是$Q_i$（\sum Q_i(z) = 1, Q_i(z)\geq 0，z是连续变量时要积分）。另外，对于某个数据点$x_i$，假设$z$的取值是$z_i$。

$$\begin{equation} \begin{aligned} \mathcal{l}(\theta) &= \sum_i log\ p(x_i;\theta) \\ &= \sum_i log \sum_{z_i} p(x_i,z_i;\theta) \\ &= \sum_i log \sum_{z_i} Q_i(z_i) \frac{p(x_i,z_i;\theta)}{Q_i(z_i)} \\ & \geq \sum_i \sum_{z_i} Q_i(z_i) log \frac{p(x_i,z_i;\theta)}{Q_i(z_i)} \end{aligned} \end{equation}$$

最后一行用了Jensen不等式，$log E_{z_i \sim Q_i} \frac{p(x_i,z_i;\theta)}{Q_i(z_i)} \geq E_{z_i \sim Q_i} log \frac{p(x_i,z_i;\theta)}{Q_i(z_i)}$。我们把最后一行的公式写成$J(Q,\theta)$，那么这就是原来似然函数的一个下界lower bound。如果给定$\theta_0$，存在$Q$使得$J(Q,\theta_0)=\mathcal{l}(\theta_0)$，那么最大化$J(Q,\theta)$相当于逐渐地增加$\mathcal{l}(\theta)$。假设$\theta' = \arg\max J(Q,\theta)$，那么$\mathcal{l}(\theta_0)=J(Q,\theta_0) \leq J(Q,\theta') \leq \mathcal{l}(\theta')$。
根据Jensen不等式，如果$J(Q,\theta)=\mathcal{l}(\theta)$成立，那么$\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}=c$是一个常数，但是$Q_i(z_i)$是一个概率分布，因此

$$\sum_{z}Q_i(z) = \frac{1}{c} \sum_{z} p(x_i,z;\theta) =1$$

因此

$$\begin{equation} \begin{aligned} Q_i(z_i) &= \frac{p(x_i,z_i;\theta)}{\sum_{z} p(x_i,z;\theta)} \\ &= \frac{p(x_i,z_i;\theta)}{p(x_i;\theta)} \\ &= p(z_i \mid x_i ; \theta) \\ \end{aligned} \end{equation}$$

因此，当$Q_i(z_i)$是给定$x_i$的后验概率，$J(Q,\theta)=\mathcal{l}(\theta)$。接下来只要最大化$J(Q,\theta)$，那么就可以迭代新的$\theta'$使得$\mathcal{l}(\theta) \leq \mathcal{l}(\theta')$。计算$Q_i(z_i) = p(z_i \mid x_i ; \theta)$就是EM的E-step，之后最大化$J(Q,\theta)$就是EM的M-step。所以本来只是最大化$\mathcal{l}(\theta)$的问题变成了$J(Q,\theta)$的坐标上升方法。可能有人会发现$Q_i(z_i)=p(z_i \mid x_i ; \theta)$，此处的$\theta$是前一次迭代得到的$\theta$，而后面优化$J(Q,\theta)$得到的$\theta$是本次迭代得到的$\theta$。这个在下面的版本中会体现得更好。

接下来是ESL上面8.5.2的版本，稍微不一样的使得问题复杂的就是一开始上述的$Q_i$就固定了。为了化简符号，我们以$X$代表要拟合的数据，$Z$是隐变量。因此要拟合的最大似然就是$\mathcal{l}(\theta) = P(X;\theta)$。然而

$$P(Z \mid X ; \theta) = \frac{ P( Z, X ; \theta) }{ P(X ; \theta) }$$
因此
$$P(X ; \theta) = \frac{ P( Z, X ; \theta) }{ P(Z \mid X ; \theta) }$$
取对数，那么就是
$$log P(X ; \theta) = log \frac{ P( Z, X ; \theta) }{ P(Z \mid X ; \theta) }$$
如果对于上式取对 $P(Z \mid X ; \theta)$的期望，因为左边不涉及 $Z$，那么可以得到
$$log P(X ; \theta) = \sum_{Z} P(Z \mid X ; \theta) log \frac{ P( Z, X ; \theta) }{ P(Z \mid X ; \theta) } ＝ E_{Z} log \frac{ P( Z, X ; \theta) }{ P(Z \mid X ; \theta) }$$
这个已经很像是上面的 $J(Q,\theta)$了。但是这边一开始就用 $P(X ; \theta) = \frac{ P( Z, X ; \theta) }{ P(Z \mid X ; \theta) }$构造了 $\mathcal{l}(\theta) = J(Q,\theta)$的情形。

然而只有两式相等的情形是不够的，需要知道如何迭代使得$\mathcal{l}(\theta)$不断上升。 接下来，假设我们已经知道了在第t次迭代$\theta_t$的值，对于上面取期望的操作，我们取的是$P(Z \mid X ; \theta_t)$的期望，因此

$$\begin{equation} \begin{aligned} log P(X ; \theta) &= \sum_{Z} P(Z \mid X ; \theta_t) log \frac{ P( Z, X ; \theta) }{ P(Z \mid X ; \theta) } \\ &= \sum_{Z} P(Z \mid X ; \theta_t) log P( Z, X ; \theta) - \sum_{Z} P(Z \mid X ; \theta_t) log P(Z \mid X ; \theta)\\ &= J_1(\theta_t, \theta) - J_2(\theta_t, \theta) \end{aligned} \end{equation}$$
接下来是证明只需要优化 $J_1(\theta_t, \theta)$当中的 $\theta$就相当于优化了 $\mathcal{l}(\theta)$。原因是 $J_2(\theta_t, \theta) \geq J_2(\theta_t,\theta_t)$。证明这个也是用Jensen不等式。
$$\begin{equation} \begin{aligned} J_2(\theta_t, \theta) - J_2(\theta_t, \theta_t) &= \sum_{Z} P(Z \mid X ; \theta_t) log \frac{ P( Z \mid X ; \theta) }{ P(Z \mid X ; \theta_t) } \\ &\geq log \sum_{Z} P(Z \mid X ; \theta_t) \frac{ P( Z \mid X ; \theta) }{ P(Z \mid X ; \theta_t) } \\ &＝ log \sum_{Z} P( Z \mid X ; \theta) \\ &=0 \end{aligned} \end{equation}$$
因此对于任意 $\theta$必有 $J_2(\theta_t, \theta) \geq J_2(\theta_t, \theta_t)$。因此如果 $\theta'$最大化了 $J_1(\theta_t, \theta)$，那么可以令 $\theta_{t+1}=\theta'$使得 $\mathcal{l}(\theta_t) \leq \mathcal{l}(\theta_{t+1})$。因此在ESL的书上算法8.2，E-step是计算 $P(Z \mid X ; \theta_t)$形成 $J_1(\theta_t, \theta)$的确切形式，M-step是优化 $J_1(\theta_t, \theta)$。而且这里可以看到最大化 $J_1(\theta_t, \theta)$是不必要的，只需要升高一点就可以了。

以上两个版本是对EM的不同解读。共同之处都是需要额外的隐变量使得问题简化。第一个版本是用Jensen不等式求出一个过当前解的下界函数，通过优化这个相对简单的函数来迭代解。而第二个版本则是直接把似然函数变成两个部分来看，第二部分必然减少但是因为负号而增加，因而只需要优化第一部分。相对来说第二版本，即ESL书上的版本，相对比较晦涩。
可能有人会发现在这个版本里面，在对数log里面分母的一项是$\theta$，是参与M-step优化的，而第一个版本里面的是$\theta_t$，是不参与M-step优化的。然而因为第二个版本只需要优化$J_1(\theta_t, \theta)$，即log里面分子的一项，所以和第一个版本殊途同归。

ESL上面8.5.3节的本质和第一个版本是一回事。这里也间接解决了课本习题8.1和8.2。吴恩达版本也是习题8.7的一个说明。

熟悉非负矩阵分解算法的人可能一下子就觉得很熟悉，其实NMF的paper里面也说到用了EM-like的方法。而且NMF里面也是通过不断放缩找到过当前解的一个上界函数来优化以迭代，思路本质上是一样的。