机器学习－由浅入深学EM算法

1. 最大似然估计到EM算法

之前在介绍《最大似然估计》中提到，MLE估计的目标是产生样本集合可能性最大的参数。其核心思想是在已知结果情况下，然后寻求使该结果出现的可能性最大的条件，并以此作为估计值。假设：要根据iris数据中的花蕊长度来估计花瓣的长度，我们可以利用最小二乘来估计参数。根据极大似然估计的定义，我们要获得产生样本可能性最大的参数$\theta$，则似然函数可以记为：

$$\ell(\theta)=\prod^m_{i=1} \frac{1}{\sqrt{2\pi}\sigma} exp(- \frac{(y^{(i)}-\theta x)^2}{2 \sigma})$$

由于线性回归服从独立同分布的假设，因此，我们可以忽略$\sigma$的影响，直接估计$\theta$。

但是，我们假设如果知道不仅花蕊的长度与花瓣的长度有关，同时，花的种类对花瓣的长度有影响，即不同种类花的花瓣长度可能服从不同分布，我们仍然在独立同分布的假设下采用极大似然估计求解是不合适的。因为，不同种类花的花瓣与花蕊长度的误差项不服从同一正态分布，即不同分布的$sigma$不同，因此，我们无法忽略该参数对目标函数的影响，需要同时估计两个参数$\theta$和$\sigma$。但是假设我们如果无法获知各个样本隶属的花的类别，即无法通过参数$\sigma$的先验知识估计参数$\sigma$，那么我们如何来估计花瓣的长度呢？

总结上面的问题，我们发现，（1）我们不知道每个样本隶属于的花的类别；（2）我们也不知道花瓣与花蕊的线性关系是什么。从数学角度考虑这些问题，可以映射成：（1）抽样样本不知道是从哪个分布抽取的；（2）每个分布下的参数是什么。

可以确认的是，只有当我们知道每个样本属于哪个分布，我们才会准确的估计每个分布下花瓣与花蕊长度的线性关系。但是，各个类别之间的样本混在一起，我们无法获知哪个样本属于哪个类别，所以就无法准确估计各分布的参数。反之，只有当我们获得各个分布参数的准确估计下，我们才能计算出各样本隶属于哪个类别。因此，无法直接采用极大似然估计求解。这个问题在数学上是一个典型的循环依赖问题，即“鸡生蛋，蛋生鸡”的问题。想要解决循环依赖问题，我们可以先给定任意一个参数$\theta$初始值，这样，通过初始值来估计参数$\sigma$，然后再用估计出的参数$\sigma$再次估计参数$\theta$，通过多次的循环迭代，最终收敛到一个最优解。

EM算法（Expectation Maximization）的核心思想就是这样，假设我们要估计初始状态未知的两个参数，同时两个参数间存在依赖关系，那么我们就可以先赋予其中一个参数初始值，得到另一个参数的估计值，然后在另一个参数估计值的基础上，重新估计第一个参数，直至收敛。在上面的问题中，如果我们先给出每个花的类别下花瓣和花蕊长度间的线性关系（$\theta$），那么我们可以根据$\theta$值计算每个样本可能属于的花的类别（花瓣长度和花蕊长度的线性误差小的类别），这个过程在EM算法中属于Expectation步骤。通过Expectation步，每个样本都有所属类别，这样我们就可以将样本分成几部分，然后采用极大似然估计的思想，对每个部分中的样本重新估计参数，即估计各类别下花瓣和花蕊长度间的线性关系，这个就是EM算法中的Maximization步骤。如此迭代，直至各个参数不在发生变化即可停止迭代。

一般情况下，我们把每个样本隶属于哪种花对应的参数称为隐含变量，因为，如果没有隐含变量的话，我们可以很直观的求解待估计参数，当存在隐含变量后，本来很简单求解的问题变得复杂。因此，EM算法存在的意义就在于此，可用于包含隐含变量的参数估计。但是，按照最大似然的思想，我们每次估计会使得目标函数上升，最终达到极值，而EM算法怎样保证算法收敛？怎样保证可以获得目标函数的极值？这些问题就必须从数学角度进行推导了～

2. EM算法原理

假设给定训练集$X={x^{(1)}, x^{(2)}, \ldots, x^{(m)}}$，我们希望求解包含隐含变量$Z＝{z^{(1)}, z^{(2)}, \ldots, z^{(m)}}$的目标函数$f(x,z;\theta)$中的参数$\theta$。根据最大似然估计，我们可以获得：

$$\theta=arg max_{\theta} \ell (\theta;X)= arg max_{\theta} log \sum_Z f(X,Z;\theta)＝arg max_{\theta} \sum_{i=1}^m log \sum_{z^{(i)}} f(x^{(i)},z^{(i)};\theta)$$

对于上述公式，想要求解$log \sum_Z f(X,Z;\theta)$形式的极值往往非常困难，如果我们可以观察到隐藏变量$z$的话，则可以消除$\sum$符号，这样求解极值就容易的多（采用梯度下降算法）。但实际上我们没法获得$z$的分布或者$z$值，只可能通过给定$\theta$的情况下，通过求解后验概率$P(z｜x;\theta)$获得参数$z$的估计值。

我们令$L(\theta)=log \sum_Z P(X,Z;\theta)$，其中$Z$为隐含变量。很明显，求解该log嵌套$\sum$形式的极值很麻烦，因此不妨引入概率分布$Q(Z)$，其中$\sum_{Z} Q(Z)=1$，且$Q(Z) \geq 0$，得到：

$$L(\theta)=log \sum_{Z} P(X,Z;\theta)=log \sum_{Z} Q(Z) \frac{P(X,Z;\theta)}{Q(Z)}$$

a. 离散变量的期望

对于$\frac{P(X,Z;\theta)}{Q(Z)}$而言，隐含变量$Z$为离散变量，如果我们假定$Q(Z)$正好是离散随机变量$Z$的概率分布函数，且$Y=g(Z)=\frac{P(X,Z;\theta)}{Q(Z)}$为连续函数，那么：

$$E(Y)=E[g(Z)]=\sum_Z g(Z)Q(Z)$$

如果$Z$是连续随机变量，那么：

$$E(Y)=E[g(Z)]=\int^{\infty}_{-\infty} g(Z)Q(Z)dZ$$

所以：

$$L(\theta)=log \sum_{Z} Q(Z) \frac{P(X,Z;\theta)}{Q(Z)}=log E_Z[g(Z)]=log E_Z [\frac{P(X,Z;\theta)}{Q(Z)}]$$

b. Jensen不等式

如果任意函数$f$为凸函数，$X$为随机变量，则$E[f(X)] \geq f(E[X])$，若上式取等号当且仅当$X$为常数，若任意函数$f$为凹函数，则不等式反向。因为log函数为严格凹函数，所以对于任意分布$Q$，

$$L(\theta)=log E_Z [\frac{P(X,Z;\theta)}{Q(Z)}] \geq E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}]$$

c. 凹函数的下界函数

对于任意凹函数$f$，假设我们想求$arg max_{\theta} f(\theta)$，对于一般优化算法而言，我们需要保证每次迭代得到的$f(\theta_{t+1})$一定比$f(\theta_t)$要大，这样才能保证优化过程呈单调逼近于极大值。

因此，如果对于任意$t$次迭代，我们能找到$f(\theta)$的下界函数$g_{\theta_t}(\theta)$，满足$f(\theta) \geq g_{\theta_t}(\theta)$，且$f(\theta_t) \geq g_{\theta_t}(\theta_t)$那么令$\theta_{t+1}:=arg max_{\theta} g_{\theta_t}(\theta)$，此时，$g_{\theta_t}(\theta_{t+1}) ＝ max_{\theta} g_{\theta_t}(\theta_t)$,即会满足每次迭代成单调趋势：

$$f(\theta_{t+1}) \geq g_{\theta_t}(\theta_{t+1}) \geq g_{\theta_t}(\theta_t)=f(\theta_t)$$

对于$t＋1$次迭代，我们仍然可以找到函数$g_{\theta_{t+1}}(\theta)$，满足$f(\theta_{t+2}) \geq g_{\theta_{t+1}}(\theta_{t+2}) \geq g_{\theta_{t+1}}(\theta_{t+1})=f(\theta_{t+1})$，据此计算出$\theta_{t+2}$。

再回到EM算法中，如果我们可以找到$L(\theta)$的下界函数，并且让不等式在$\theta_t$处取等号，即可使得每次迭代的目标函数更优。我们已知$L(\theta) \geq E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}]$，所以可以看出$E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}]$为下界函数，但是我们还要让下界函数在$\theta=\theta_t$处取等号，即要求$\frac{P(X,Z;\theta)}{Q(Z)}$为常数。若要使得$Q(Z) \propto P(X,Z;\theta_t)$，且$\sum Q(Z) = 1$，则根据全概率公式可知，仅当$Q(Z)=\frac{P(X,Z;\theta_t)}{\sum_Z P(X,Z;\theta_t)}=P(Z｜X;\theta_t)$时，满足$Q(Z)$与$P(X,Z;\theta_t)$呈正比关系。

因此，我们令$Q(Z)=P(Z｜X;\theta_t)$，此时：

$$L(\theta) \geq E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}] ＝ E_{Z \sim Q = Z|X;\theta_t} [log \frac{P(X,Z;\theta)}{P(Z|X;\theta_t)}]$$

显然，上式在$\theta = \theta_t$时，$L(\theta_t) ＝ E_{Z｜X;\theta_t} [log \frac{P(X,Z;\theta_t)}{P(Z｜X;\theta_t)}] = E_{Z｜X;\theta_t} [log \sum_Z P(X,Z;\theta_t)]$可以取等号。

综上所述，$E_{Z｜X;\theta_t} [log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)}]$为$L(\theta)$的下界函数，我们可以通过每次迭代计算找到满足：

$$\theta_{t+1} := arg max_{\theta} E_{Z｜X;\theta_t} [log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)}]$$

d. $\theta$更新公式推导

找到目标函数的下界函数$g_{\theta_t}(\theta)$后，我们可以通过更新$\theta_{t+1}:=arg max_{\theta} g_{\theta_t}(\theta)$实现每次的迭代，并可以保证迭代过程呈单调变化趋势，最终收敛。

我们可以再对更新公式进行变换，使得计算更为简便：

$$\theta_{t+1} = arg max g_{\theta_t} (\theta) \\ = arg max E_{Z｜X;\theta_t} \, [log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)}] \\ = arg max \sum_Z P(Z｜X;\theta_t) log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)} \\ = arg max \sum_Z [P(Z｜X;\theta_t)logP(X,Z;\theta) - P(Z｜X;\theta_t)log P(Z｜X;\theta_t)] \, (\because 与 \theta无关) \\ = arg max \sum_Z [P(Z｜X;\theta_t)logP(X,Z;\theta) \\ = arg max E_{Z｜X;\theta_t} log P(X,Z;\theta) 。$$

3. EM算法过程

根据上面的推导过程，我们可以定义EM算法的Expectation步和Maximization步。

Repeat until converge｛

1. E－step：计算$P(Z｜X;\theta_t)$得到下界函数

$$L(\theta)=E_{Z｜X;\theta_t} \, \, log P(X,Z;\theta)=\sum_Z P(Z｜X;\theta_t) log P(X,Z;\theta)$$

2. M－step：更新$\theta$

$$\theta_{t+1} := arg max_{\theta} E_{Z｜X;\theta_t} \, \, log P(X,Z;\theta)$$

｝

4. 总结

总的来说，Expectation Maximization算法是一个比较难懂的算法，中间的推导过程较为复杂，可能学习一遍两遍很难把这个算法吃透，需要比较深入的统计学知识才能搞明白，这篇博文也写了三天才把算法原理过程写完。

EM算法是用于对含有隐含变量的目标函数做最大似然估计的一种方法。EM算法的应用特别广泛，例如，Kmeans聚类，高斯混合模型，隐马尔可夫模型等。后面的博文我会专门撰写一篇EM算法实例，用于具体说明EM算法求解过程。