EM算法学习

前言

$EM$算法是一种迭代优化算法，常用于具有隐变量的混合模型参数估计，在机器学习中有极为广泛的用途。说到参数估计，一般首先想到的是概率论上的极大似然估计，$EM$算法其实就是为了解决那些直接用极大似然函数求导无法解决的参数估计问题。

极大似然估计

假设现在有一个概率分布函数$p(x|\theta )$，$\theta$代表这个分布函数的参数(例如，对于一维高斯分布函数来说，$\theta$就代表的是 θ = { μ , σ 2 } \theta = \{\mu,\sigma^2\} θ={μ,σ2})。然后有$N$个样本数据 X = { x 1 , x 2 , ⋯   , x N } X=\{x_1,x_2,\cdots,x_N\} X={x1​,x2​,⋯,xN​}是从这个概率分布中采样得到的，而且这些样本数据都是相互独立的。那么极大似然估计函数定义为：
$$L(\theta )=\prod _i^{N}p(x_i|\theta )=L(\theta |X)$$
极大似然估计，通俗理解来说，就是利用已知的样本信息，反推最具有可能（最大概率）导致这些样本出现的模型参数值。换句话说，极大似然估计提供了一种给定观察数据来评估模型参数的方法，即：“模型已定，参数未知”。

极大似然问题就是希望找到一个$\theta$使似然函数最大，用公式来表达即：
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}L(\theta )$$
为了计算的方便一般是最大化$lnL(\theta )$，即
$$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}lnL(\theta )$$
然而计算参数$\theta$的难易往往取决于$L(\theta )$的形式。如果$L(\theta )$的形式相对简单，例如$p(x|\theta )$代表的是一个一维高斯分布，那它的参数 θ = { μ , σ 2 } \theta = \{\mu,\sigma^2\} θ={μ,σ2}求起来就相对来说比较容易。只需要对极大似然函数求参数的偏导，然后令偏导数为$0$即可。
$$\begin{array}{rl}lnL(\theta )& =ln\prod _{i=1}^n[\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x_i-\mu {)}^2}{2{\sigma }^2})]\\ & =\sum _{i=1}^{n}ln[\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{(x_i-\mu {)}^2}{2{\sigma }^2})]\\ & =\sum _{i=1}^n[-ln\sqrt{2\pi }\sigma -\frac{(x_i-\mu {)}^2}{2{\sigma }^2}]\\ & \\ \frac{\partial lnL(\theta )}{\partial \mu }& =\frac{\partial ({\sum }_{i=1}^n[-ln\sqrt{2\pi }\sigma -\frac{(x_i-\mu {)}^2}{2{\sigma }^2}])}{\partial \mu }\\ & =\frac{\partial ({\sum }_{i=1}^n[-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}])}{\partial \mu }\\ & =\sum _{i=1}^n\frac{(x_i-\mu )}{{\sigma }^2}\end{array}$$
令$\frac{\partial lnL(\theta )}{\partial \mu }=0$得：
$${\mu }_{MLE}=\frac{1}{n}\sum _{i=1}^n{x}_i$$
对$\sigma$的估计这里略过。

但是，如果$L(\theta )$的形式比较复杂，例如$GMM$，其中含有隐变量，这种直接求偏导的方式是没法求得解析解的，这种情况就可以借助$EM$算法来求解。

$EM$算法

假设我们的模型有两个随机变量$X,Z$，两个随机变量的联合分布函数为$p(x,z|\theta )$，其中$x$和$z$代表的是$X,Z$中的观测值，$\theta$是联合分布函数的未知参数。假设我们现在只可观测到$X$中的$x$，$Z$中的$z$是无法被观测到的隐变量，我们只能假设它是服从某种分布，但我们仍然希望只通过观测到的$X$来估计出模型中的参数$\theta$。因为没有观测到$Z$，所以我们的似然函数就要通过$X$的边缘密度函数来进行计算了：
$$L(\theta )=logp(x|\theta )={\int }_{z}logp(x,z|\theta )dz$$
在这样的模型中，由于积分的存在，最大化这个函数可能得不到解析解。 对于这种情况，$EM$算法可以提供一种计算该函数参数$\theta$ 局部最大值的方法。

$EM$算法的步骤

$EM$算法是一种当数据不完整或有隐变量时，从给定的数据集中找到潜在分布参数的最大似然估计的通用方法。$EM$算法通过$E—Step$和$M—Step$迭代交替求解极大似然参数估计问题。$EM$算法整个求解的过程中都是在更新参数$\theta$，使之收敛到一个局部极大值。

$E—Step$ ： 在固定${\theta }^t$(${\theta }^t$代表第$t$次迭代的参数)情况下计算隐变量的条件密度函数$p(z|x,{\theta }^t)$ ，然后计算联合密度似然函数在$p(z|x,{\theta }^t)$分布下的期望$Q(\theta ,{\theta }^t)$(后面简称$Q$函数)。$E—step$就是$Expectation—Step$，就是求期望。
$$\begin{array}{rl}Q(\theta ,{\theta }^t)& =E_{p(z|x,{\theta }^t)}[logp(x,z|\theta )]\\ & ={\int }_z[logp(x,z|\theta )]p(z|x,{\theta }^t)dz\end{array}$$
$M—Step$: 找出$\theta$使$Q(\theta ,{\theta }^t)$最大，$M—Step$就是$Maximization—Step$，最大化$Q(\theta ,{\theta }^t)$，而求得极值
$${\theta }^{t+1}=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^t)$$
算法在$E—Step$和$M—Step$交替，直到$\theta$的估计收敛到一个局部最大值。

$EM$算法本质

在深入了解$EM$算法之前，让我们先了解两个重要的概念，一个是$ELBO$，一个是坐标上升法[3]。$EM$算法本质上就是$ELBO$的坐标轴上升算法。

$ELBO$是什么

要了解什么是$evidencelowerbound(ELBO)$证据下界，我们首先要了解$evidence$是什么意思。$evidence$指的是在固定模型参数$\theta$情况下得到的似然估计函数[4]：
$$evidence:=logp(x|\theta )$$
为什么这样的量会被称为$evidence$？我们可以这样想，如果我们选择模型的分布$p$和参数$\theta$比较准确，那根据观测到的样本数据计算出来的似然函数$logp(x|\theta )$应该会比较大。因此，一个更大的似然函数$logp(x|\theta )$在某种程度上“证明”了我们现在选择的模型对于生成现有数据有更大的可能。

$ELBO$其实就是似然函数的下界，如果我们知道隐变量$Z$是服从某种分布(用$q(z)$表示)，而且满足$p(x,z|\theta )=p(x|z,\theta )q(z)$，那么$ELBO$就可以表示成：
$$ELBO:=E_{Z\sim q}[log\frac{p(x,z|\theta )}{q(z)}]$$

$ELBO$的推导：
$$\begin{array}{rl}logp(x|\theta )& =log{\int }_{z}p(x,z|\theta )\\ & =log{\int }_{z}p(x,z|\theta )\frac{q(z)}{q(z)}dz\\ & =log{E}_{Z\sim q}[\frac{p(x,z|\theta )}{q(z)}]\\ & \ge \underset{ELBO}{\underbrace{E_{Z\sim q}[log\frac{p(x,z|\theta )}{q(z)}]}}(Jensen不等式)\end{array}$$
事实上，$evidence$与$ELBO$还有一个非常有意思的关系，那就是它们之间的差正好是$p(z|x,\theta )$和$q(z)$的$KL$散度(即相对熵，表示两个分布之间差异的度量，但这里不讨论)

公式推导：
$$\begin{array}{rl}KL(q(z)||p(z|x,\theta ))& =E_{Z\sim q}[log\frac{q(z)}{p(z|x,\theta )}]\\ & =E_{Z\sim q}[logq(z)]-E_{Z\sim q}[logp(z|x,\theta )]\\ & =E_{Z\sim q}[logq(z)]-E_{Z\sim q}[log\frac{p(x,z|\theta )}{p(x|\theta )}]\\ & =E_{Z\sim q}[logq(z]-E_{Z\sim q}[logp(x,z|\theta )]+E_{Z\sim q}[logp(x|\theta )]\\ & =logp(x|\theta )-E_{Z\sim q}[log\frac{p(x,z|\theta )}{q(z)}]\\ & =evidence-ELBO\end{array}$$

坐标上升法

坐标上升法是一个相对简单的最大化函数的迭代优化策略，每次通过更新函数中的一维，通过多次的迭代以达到优化函数的目的。 例如，给定一个有两个变量$a,b$的函数$f(a,b)$，先固定变量$a$维为$\hat{a}$，然后找到$b$维的$\hat{b}$使得$f(a,b)$取得最大值，然后固定$b$维为$\hat{b}$，寻找$a$维函数的最大值。这两步不断循环迭代，直到收敛到最大值。

​ 坐标上升法

$EM$：$ELBO$的坐标上升算法

在我们假设的模型中，参数$\theta$是未知的，而且隐变量$Z$的分布$q(z)$也是未知的，那么$ELBO$就可以看成是一个关于$\theta$和$q(z)$的函数:
$$F(q,\theta ):=E_{Z\sim q}[log\frac{p(x,z|\theta )}{q(z)}]$$
$EM$算法就是关于$F(q,\theta )$函数的坐标上升算法。$E—Step$相当于固定参数$\theta$最大化$F(q,\theta )$去求解$q(z)$，$M—Step$相当于固定$q(z)$去最大化$F(q,\theta )$求解参数$\theta$。

首先，让我们把参数$\theta$固定为${\theta }^t$最大化$F(q,\theta )$来求解隐变量$Z$的分布$q(z)$
$$\begin{array}{rl}\underset{q}{\mathrm{argmax}}F(q,{\theta }^t)& =\underset{q}{\mathrm{argmax}}{E}_{Z\sim q}[log\frac{p(x,z|{\theta }^t)}{q(z)}]\\ & =\underset{q}{\mathrm{argmin}}KL(q(z)||p(z|x,{\theta }^t))(ELBO与KL(q(z)||p(z|x,{\theta }^t))的关系）\\ & =p(z|x,{\theta }^t)(KL散度的性质，当q(z)=p(z|x,{\theta }^t)时，KL散度最小)\end{array}$$
注意，计算$p(z|x,{\theta }^t)$正是在$E—Step$ 中构造$Q$函数所需要的！

接下来，固定$q^t(z)=p(z|x,{\theta }^t)$最大化$F(q^t,\theta )$，由于$q^t(z)$固定，那其实最大化$F(q^t,\theta )$就等同于最大化$Q$函数：
$$\begin{array}{rl}\underset{q}{\mathrm{argmax}}F(q^t,\theta )& =\underset{q}{\mathrm{argmax}}{E}_{Z\sim q}[log\frac{p(x,z|\theta )}{p(z|x,{\theta }^t)}]\\ & =\underset{q}{\mathrm{argmax}}{E}_{Z\sim q}[logp(x,z|\theta )]-\underset{constant}{\underbrace{E_{Z\sim q}[p(z|x,{\theta }^t)]}}\\ & =\underset{q}{\mathrm{argmax}}\underset{Q(\theta ,{\theta }^t)}{\underbrace{E_{p(z|x,{\theta }^t)}[logp(x,z|\theta )]}}\end{array}$$
最大化$Q$函数就是$EM$算法中的$M—Step$。

$EM$算法图示

首先，让我们重申一下$F(q,\theta )$函数和似然函数$L(\theta )$的一些关系。

第一，在上面介绍$ELBO$(即$F(q,\theta )$)的时候已经介绍了$ELBO$和似然函数的关系，就是它们之间的差正好是$KL(q(z)||p(z|x,\theta ))$。$KL$散度有一个性质，当$p(z|x,\theta )=q(z)$时，$KL(q(z)||p(z|x,\theta ))=0$。所以$EM$在第$t$次迭代的时候，如果有$KL(q(z)||p(z|x,{\theta }^t))=0$，那$F(q,{\theta }^t)$与$L({\theta }^t)$相等，故这两个函数应该在${\theta }^t$相交。

第二，固定$q(z)$为$q_t(z)$，函数$F(q_t,\theta )$在${\theta }^{t+1}$处一定大于或等于在${\theta }^t$处，即
$$F(q_t,{\theta }^{t+1})\ge F(q_t,{\theta }^t)$$
这是因为${\theta }^{t+1}$是算法中的中最大化$F(q_t,\theta )$得到的。

第三，由$ELBO$(即$F(q,\theta )$)和似然函数的关系，则必有
$$L(\theta )\ge F(q_t,\theta )$$
上面的三个性质可由下图表示：

​ $F(q,\theta )$与$L(\theta )$的关系

基于上面的三个性质，我们可以推断$EM$算法通过在似然函数$L(\theta )$的“攀爬”最终会收敛到$L(\theta )$的一个局部极值。

$EM$算法的过程：在第$t$次迭代中，$F(q_t,\theta )$ 低于$L(\theta )$函数，且在${\theta }_t$与之相交， 此时，$F(q_t,{\theta }^t)=L({\theta }^t)$。 由于$F(q_t,{\theta }_{t+1})\ge F(q_t,{\theta }_t)$，此时$L({\theta }^{t+1})\ge L({\theta }^t)$。在下一次迭代中，我们可以通过$F(q_{t+1},{\theta }^{t+1})$计算${\theta }^{t+1}$ 处的似然函数。 因此，每次迭代$F(q_{t+1},{\theta }^{t+1})$ 后，在似然函数上都有一个新的、比当前值更高的值。 这个过程如下图所示：

​ EM算法的迭代过程

参考文献或网页

[1] The EM Algorithm Explained. https://medium.com/@chloebee/the-em-algorithm-explained-52182dbb19d9

[2] Expectation-Maximization(EM) 算法 http://blog.tomtung.com/2011/10/em-algorithm/

[3] Expectation-maximization: theory and intuition https://mbernste.github.io/posts/em/

[4] The evidence lower bound (ELBO) https://mbernste.github.io/posts/elbo/