最大期望算法

统计学学术名词

最大期望算法（Expectation-Maximization algorithm, EM），或Dempster-Laird-Rubin算法，是一类通过迭代进行极大似然估计（Maximum Likelihood Estimation, MLE）的优化算法，通常作为牛顿迭代法（Newton-Raphson method）的替代用于对包含隐变量（latent variable）或缺失数据（incomplete-data）的概率模型进行参数估计。

历史

对EM算法的研究起源于统计学的误差分析（error analysis）问题。1886年，美国数学家Simon Newcomb在使用高斯混合模型（Gaussian Mixture Model, GMM）解释观测误差的长尾效应时提出了类似EM算法的迭代求解技术。在极大似然估计（Maximum Likelihood Estimation, MLE）方法出现后，英国学者Anderson McKendrick在1926年发展了Newcomb的理论并在医学样本中进行了应用。1956年，Michael Healy和Michael Westmacott提出了统计学试验中估计缺失数据的迭代方法，该方法被认为是EM算法的一个特例。1970年，B. J. N. Blight使用MLE对指数族分布的I型删失数据（Type I censored data）进行了讨论。Rolf Sundberg在1971至1974年进一步发展了指数族分布样本的MLE并给出了迭代计算的完整推导。

EM算法的正式提出来自美国数学家Arthur Dempster、Nan Laird和Donald Rubin，其在1977年发表的研究对先前出现的作为特例的EM算法进行了总结并给出了标准算法的计算步骤，EM算法也由此被称为Dempster-Laird-Rubin算法。1983年，美国数学家吴建福（C.F. Jeff Wu）给出了EM算法在指数族分布以外的收敛性证明。

此外，在二十世纪60-70年代对隐马尔可夫模型（Hidden Markov Model, HMM）的研究中，Leonard E. Baum提出的基于MLE的HMM参数估计方法，即Baum-Welch算法（Baum-Welch algorithm）也是EM算法的特例之一。

理论

EM算法是基于极大似然估计（Maximum Likelihood Estimation, MLE）理论的优化算法。给定相互独立的观测数据，和包含隐变量、参数的概率模型，根据MLE理论，的最优单点估计在模型的似然取极大值时给出：。考虑隐变量，模型的似然有如下展开：

隐变量可以表示缺失数据，或概率模型中任何无法直接观测的随机变量，式中第一行是隐变量为连续变量的情形，第二行是隐变量为离散变量的情形，积分/求和的部分也被称为的联合似然（joint liklihood）。不失一般性，这里按离散变量为例进行说明。由MLE的一般方法，对上式取自然对数后可得：

上述展开考虑了观测数据的相互独立性。引入与隐变量有关的概率分布，即隐分布（可认为隐分布是隐变量对观测数据的后验，参见标准算法的E步推导），由Jensen不等式，观测数据的对数似然有如下不等关系：

当使不等式右侧取全局极大值时，所得到的至少使不等式左侧取局部极大值。因此，将不等式右侧表示为后，EM算法有如下求解目标：式中的等效于MM算法（Minorize-Maximization algorithm）中的代理函数（surrogate function），是MLE优化问题的下限，EM算法通过最大化代理函数逼近对数似然的极大值。

算法

标准算法

计算框架

EM标准算法是一组迭代计算，迭代分为两部分，即E步和M步，其中E步“固定”前一次迭代的，求解使取极大值；M步使用求解使取极大值。EM算法在初始化模型参数后开始迭代，迭代中E步和M步交替进行。由于EM算法的收敛性仅能确保局部最优，而不是全局最优。因此通常对EM算法进行随机初始化并多次运行，选择对数似然最大的迭代输出结果。以下给出EM算法E步和M步的推导。

1. E步（Expectation-step, E-step）

由EM算法的求解目标可知，E步有如下优化问题：

考虑先前的不等关系，这里首先对进行展开：

注意到，推导上式时考虑了：。由贝叶斯定理（Bayes' theorem），上式可化为：

式中为Kullback-Leibler散度（Kullback-Leibler divergence, KL）或相对熵（relative entropy），表示吉布斯自由能（Gibbs free energy），即由Jensen不等式得到的代理函数等价于隐分布的自由能。求解的极大值等价于求解隐分布自由能的极大值，即隐分布对隐变量后验的KL散度的极小值。由KL散度的性质可知，其极小值在两个概率分布相等时取得，因此当时，取极大值，对EM算法的第次迭代，E步有如下计算：

2. M步（Maximization step, M-step）

在E步的基础上，M步求解模型参数使取极大值。该极值问题的必要条件是：

式中表示联合似然对隐分布的数学期望。在为凸函数时（例如隐变量和观测服从指数族分布），上述推导也是充分的。由此得到M步的计算：

计算步骤

将统计模型带入EM算法的计算框架即可得到其计算步骤。这里以高斯混合模型（Gaussian Mixture Model, GMM）为例进行说明。由GMM的一般定义可知，其似然和参数有如下表示：

根据学习数据的维度，式中表示均值为，方差/协方差为的正态分布/联合正态分布。为正态分布的混合比例，为参与混合的分布总数。定义与观测数据有关的隐变量：，令隐分布表示GMM聚类的软指定（soft assignment），即每个数据来源于第个分布的概率，则隐变量有离散取值。

将上述内容带入EM算法的计算框架后，E步有如下展开：

GMM中有：，因此E步的计算步骤为：

M步通过E步输出的隐变量后验计算模型参数，在GMM中，M步计算框架的优化问题有如下表示：

不失一般性，带入单变量正态分布的解析形式后对模型参数求偏导数可得M步的计算步骤：

根据以上计算步骤，这里给出一个在Python 3环境使用EM算法求解GMM的编程实现：

改进算法

基于贝叶斯推断（Bayesian inference）的EM算法

参见：变分贝叶斯EM

在MLE理论下，EM算法仅能给出模型参数的单点估计，引入贝叶斯推断方法后，EM算法能够给出模型参数的后验（posterior）分布避免过度拟合，其中常见的例子是极大后验估计（Maximum A Posteriori estimation, MAP）的EM算法。MAP-EM在标准EM算法的基础上引入了模型参数的先验（prior）：，此时MAP-EM的优化目标由模型的似然转变为后验，其离散形式可表示为：

类比标准EM算法，考虑隐分布后，由Jensen不等式可得到对数后验的代理函数，即隐变量的自由能：

由此可得MAP-EM的迭代步骤：

MAP-EM在Dempster et al. (1977)就已被提出，但不同于标准EM，MAP-EM的隐分布是隐变量和模型参数的联合分布，其对应的隐变量后验往往没有解析形式。在贝叶斯体系下，求解该隐变量后验的方法包括马尔可夫链蒙特卡罗（Markov Chain Monte Carlo, MCMC）和变分贝叶斯估计（Variational Bayesian Inference, VBI），对前者，可证明由MCMC求解的MAP-EM等价于吉布斯采样（Gibbs sampling）算法；对后者，由VBI求解的MAP-EM被称为变分贝叶斯EM算法（Variational Bayesian EM algorithm, VBEM）。

VBEM使用平均场理论（Mean Field Theory, MFT）将隐分布近似为其在各个维度上分布的乘积：并由此得到以下迭代步骤：

使用VBEM的常见例子是语言建模问题中的隐含狄利克雷分布（latent dirichlet allocation）。

EM梯度算法（EM gradient algorithm）和广义EM算法（Generalized EM algorithm, GEM）

参见：广义期望最大算法

EM算法的M步通过计算偏导数求解代理函数的极大值，EM梯度算法（EM gradient algorithm）将该过程替换为牛顿迭代法（Newton-Raphson method）以加速迭代收敛。更进一步地，当代理函数不是凸函数或无法有效地对求解极大值时，可以使用广义EM算法（GEM）。GEM有两种实现方式，一是在M步使用非线性优化策略，例如共轭梯度算法（conjugate gradients algorithm），二是将原M步的求导计算分解为多个条件极值问题逐个计算模型参数，后者也被称为最大条件期望算法（Expectation Conditional Maximization algorithm, ECM）。

EM算法的E步也可以按ECM的方法分解为条件极值问题，由先前推导可知，E步的优化问题仅有一个全局极大值，即隐分布，因此在E步将MLE的优化目标：联合似然对观测样本按因子展开并对每个展开分别使用EM算法，可以得到同样的优化结果。对于M步，如果观测样本来自指数族分布，则M步也可以在每次迭代仅对有限个样本的展开进行。在指数族问题中使用EM算法的上述推广，可以避免在迭代中反复处理整个观测样本集，降低计算开销。

α-EM算法（α-EM algorithm）

α-EM算法是对标准算法的隐变量概率分布引入权重系数的改进版本。标准的EM算法是α-EM算法在时的特例。给定恰当的超参数，α-EM能够比标准EM算法更快收敛。有研究将α-EM算法应用于神经网络的概率学习和隐马尔可夫模型的参数估计。

性质

收敛性与稳定性：EM算法必然收敛于对数似然的局部极大值或鞍点（saddle point），其证明考虑如下不等关系：

由上式可知EM算法得到的对数似然是单调递增的，即从次迭代到次迭代，EM算法至少能维持当前的优化结果，不会向极大值的相反方向运动，因此EM算法具有数值稳定性（numerical stablility）。上述不等关系也被用于EM算法迭代终止的判定，给定计算精度，当时迭代结束。EM算法收敛性的具体证明参见Wu (1983)。

计算复杂度：在E步具有解析形式时，EM算法是一个计算复杂度和存储开销都很低的算法，可以在很小的计算资源下完成计算。在E步不具有解析形式或使用MAP-EM时，EM算法需要结合其它数值方法，例如变分贝叶斯估计或MCMC对隐变量的后验分布进行估计，此时的计算开销取决于问题本身。

与其它算法的比较：相比于梯度算法，例如牛顿迭代法和随机梯度下降（Stochastic Gradient Descent, SGD），EM算法的优势是其求解框架可以加入求解目标的额外约束，例如在高斯混合模型的例子中，EM算法在求解协方差时可以确保每次迭代的结果都是正定矩阵。EM算法的不足在于其会陷入局部最优，在高维数据的问题中，局部最优和全局最优可能有很大差异。

应用

EM算法及其改进版本被用于机器学习算法的参数求解，常见的例子包括高斯混合模型（Gaussian Mixture Model, GMM）、概率主成份分析（probabilistic Principal Component Analysis）、隐马尔可夫模型（Hidden Markov Model, HMM）等非监督学习算法。EM算法可以给出隐变量，即缺失数据的后验，因此在缺失数据问题（incomplete-data probelm）中有应用。

全国各地天气预报查询

上海市

云南省

云南省

云南省

云南省

云南省