最大期望算法

编辑 锁定 讨论999
最大期望算法(Expectation-Maximization algorithm, EM),或Dempster-Laird-Rubin算法 [1]  ,是一类通过迭代进行极大似然估计(Maximum Likelihood Estimation, MLE)的优化算法 [2]  ,通常作为牛顿迭代法(Newton-Raphson method)的替代用于对包含隐变量(latent variable)或缺失数据(incomplete-data)的概率模型进行参数估计 [2-3] 
EM算法的标准计算框架由E步(Expectation-step)和M步(Maximization step)交替组成,算法的收敛性可以确保迭代至少逼近局部极大值 [4]  。EM算法是MM算法(Minorize-Maximization algorithm)的特例之一,有多个改进版本,包括使用了贝叶斯推断的EM算法、EM梯度算法、广义EM算法等 [2] 
由于迭代规则容易实现并可以灵活考虑隐变量 [3]  ,EM算法被广泛应用于处理数据的缺测值 [1-2]  ,以及很多机器学习(machine learning)算法,包括高斯混合模型(Gaussian Mixture Model, GMM) [5] 隐马尔可夫模型(Hidden Markov Model, HMM) [6]  的参数估计。
中文名
最大期望算法
外文名
Expectation Maximization algorithm, EM
类    型
优化算法
提出者
Arthur Dempster,Nan Laird,
 
Donald Rubin 等
提出时间
1977年
应    用
数据分析,机器学习

最大期望算法历史

编辑
对EM算法的研究起源于统计学的误差分析(error analysis)问题。1886年,美国数学家Simon Newcomb在使用高斯混合模型(Gaussian Mixture Model, GMM)解释观测误差的长尾效应时提出了类似EM算法的迭代求解技术 [7]  。在极大似然估计(Maximum Likelihood Estimation, MLE)方法出现后,英国学者Anderson McKendrick在1926年发展了Newcomb的理论并在医学样本中进行了应用 [8]  。1956年,Michael Healy和Michael Westmacott提出了统计学试验中估计缺失数据的迭代方法 [9]  ,该方法被认为是EM算法的一个特例 [2]  。1970年,B. J. N. Blight使用MLE对指数族分布的I型删失数据(Type I censored data)进行了讨论 [10]  。Rolf Sundberg在1971至1974年进一步发展了指数族分布样本的MLE并给出了迭代计算的完整推导 [11-12] 
EM算法的正式提出来自美国数学家Arthur Dempster、Nan Laird和Donald Rubin,其在1977年发表的研究对先前出现的作为特例的EM算法进行了总结并给出了标准算法的计算步骤,EM算法也由此被称为Dempster-Laird-Rubin算法 [1]  [2]  。1983年,美国数学家吴建福(C.F. Jeff Wu)给出了EM算法在指数族分布以外的收敛性证明 [4] 
此外,在二十世纪60-70年代对隐马尔可夫模型(Hidden Markov Model, HMM)的研究中,Leonard E. Baum提出的基于MLE的HMM参数估计方法,即Baum-Welch算法(Baum-Welch algorithm)也是EM算法的特例之一 [6]  [13-14] 

最大期望算法理论

编辑
EM算法是基于极大似然估计(Maximum Likelihood Estimation, MLE)理论的优化算法。给定相互独立的观测数据
,和包含隐变量
、参数
的概率模型
,根据MLE理论,
的最优单点估计在模型的似然取极大值时给出:
。考虑隐变量,模型的似然有如下展开 [15]  [3] 
隐变量可以表示缺失数据,或概率模型中任何无法直接观测的随机变量,式中第一行是隐变量为连续变量的情形,第二行是隐变量为离散变量的情形,积分/求和的部分也被称为
的联合似然(joint liklihood)。不失一般性,这里按离散变量为例进行说明。由MLE的一般方法,对上式取自然对数后可得 [3] 
上述展开考虑了观测数据的相互独立性。引入与隐变量有关的概率分布
,即隐分布(可认为隐分布是隐变量对观测数据的后验,参见标准算法的E步推导),由Jensen不等式,观测数据的对数似然有如下不等关系 [3] 
使不等式右侧取全局极大值时,所得到的
至少使不等式左侧取局部极大值。因此,将不等式右侧表示为
后,EM算法有如下求解目标 [3] 
式中的
等效于MM算法(Minorize-Maximization algorithm)中的代理函数(surrogate function),是MLE优化问题的下限,EM算法通过最大化代理函数逼近对数似然的极大值。

最大期望算法算法

编辑

最大期望算法标准算法

计算框架
EM的计算框架:对数似然(蓝),E步(绿),M步(实心点) EM的计算框架:对数似然(蓝),E步(绿),M步(实心点) [16]
EM标准算法是一组迭代计算,迭代分为两部分,即E步和M步,其中E步“固定”前一次迭代的
,求解
使
取极大值;M步使用
求解
使
取极大值。EM算法在初始化模型参数后开始迭代,迭代中E步和M步交替进行。由于EM算法的收敛性仅能确保局部最优,而不是全局最优 [3-4]  。因此通常对EM算法进行随机初始化并多次运行,选择对数似然最大的迭代输出结果 [3]  。以下给出EM算法E步和M步的推导。
1. E步(Expectation-step, E-step)
由EM算法的求解目标可知,E步有如下优化问题 [3] 
考虑先前的不等关系,这里首先对
进行展开 [3] 
注意到,推导上式时考虑了:
。由贝叶斯定理(Bayes' theorem),上式可化为 [3] 
E步:优化代理损失(左),原优化目标(右) E步:优化代理损失(左),原优化目标(右) [15]
式中
为Kullback-Leibler散度(Kullback-Leibler divergence, KL)或相对熵(relative entropy),
表示吉布斯自由能(Gibbs free energy),即由Jensen不等式得到的代理函数等价于隐分布的自由能。求解
的极大值等价于求解隐分布自由能的极大值,即隐分布对隐变量后验
的KL散度的极小值。由KL散度的性质可知,其极小值在两个概率分布相等时取得,因此当
时,
取极大值,对EM算法的第
次迭代,E步有如下计算 [3]  [17] 
2. M步(Maximization step, M-step)
在E步的基础上,M步求解模型参数使
取极大值。该极值问题的必要条件是
[3] 
式中
表示联合似然
对隐分布
数学期望。在
凸函数时(例如隐变量和观测服从指数族分布),上述推导也是充分的 [3]  。由此得到M步的计算:
计算步骤
将统计模型带入EM算法的计算框架即可得到其计算步骤。这里以高斯混合模型(Gaussian Mixture Model, GMM)为例进行说明。由GMM的一般定义可知,其似然和参数有如下表示 [3]  [5] 
根据学习数据的维度,式中
表示均值为
,方差/协方差为
的正态分布/联合正态分布。
为正态分布的混合比例,
为参与混合的分布总数。定义与观测数据有关的隐变量:
,令隐分布
表示GMM聚类的软指定(soft assignment),即每个数据来源于第
个分布的概率,则隐变量有离散取值
将上述内容带入EM算法的计算框架后,E步有如下展开 [17] 
GMM中有:
,因此E步的计算步骤为:
M步通过E步输出的隐变量后验计算模型参数,在GMM中,M步计算框架的优化问题有如下表示 [17] 
不失一般性,带入单变量正态分布的解析形式后对模型参数求偏导数可得M步的计算步骤 [17] 
根据以上计算步骤,这里给出一个在Python 3环境使用EM算法求解GMM的编程实现:
# 导入模块
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
# 构建测试数据
N = 200; pi1 = np.array([0.6, 0.3, 0.1])
mu1 = np.array([[0,4], [-2,0], [3,-3]])
sigma1 = np.array([[[3,0],[0,0.5]], [[1,0],[0,2]], [[.5,0],[0,.5]]])
gen = [np.random.multivariate_normal(mu, sigma, int(pi*N)) for mu, sigma, pi in zip(mu1, sigma1, pi1)]
X = np.concatenate(gen)
# 初始化: mu, sigma, pi = 均值, 协方差矩阵, 混合系数
theta = {}; param = {}
theta['pi'] = [1/3, 1/3, 1/3]            # 均匀初始化
theta['mu'] = np.random.random((3, 2))   # 随机初始化
theta['sigma'] = np.array([np.eye(2)]*3) # 初始化为单位正定矩阵
param['k'] = len(pi1); param['N'] = X.shape[0]; param['dim'] = X.shape[1]
# 定义函数
def GMM_component(X, theta, c):
    '''
    由联合正态分布计算GMM的单个成员
    '''
    return theta['pi'][c]*multivariate_normal(theta['mu'][c], theta['sigma'][c, ...]).pdf(X)

def E_step(theta, param):
    '''
    E步:更新隐变量概率分布q(Z)。
    '''
    q = np.zeros((param['k'], param['N']))
    for i in range(param['k']):
        q[i, :] = GMM_component(X, theta, i)
    q /= q.sum(axis=0)
    return q

def M_step(X, q, theta, param):
    '''
    M步:使用q(Z)更新GMM参数。
    '''
    pi_temp = q.sum(axis=1); pi_temp /= param['N'] # 计算pi
    mu_temp = q.dot(X); mu_temp /= q.sum(axis=1)[:, None] # 计算mu
    sigma_temp = np.zeros((param['k'], param['dim'], param['dim']))
    for i in range(param['k']):
        ys = X - mu_temp[i, :]
        sigma_temp[i] = np.sum(q[i, :, None, None]*np.matmul(ys[..., None], ys[:, None, :]), axis=0)
    sigma_temp /= np.sum(q, axis=1)[:, None, None] # 计算sigma
    theta['pi'] = pi_temp; theta['mu'] = mu_temp; theta['sigma'] = sigma_temp
    return theta

def likelihood(X, theta):
    '''
    计算GMM的对数似然。
    '''
    ll = 0
    for i in range(param['k']):
        ll += GMM_component(X, theta, i)
    ll = np.log(ll).sum()
    return ll

def EM_GMM(X, theta, param, eps=1e-5, max_iter=1000):
    '''
    高斯混合模型的EM算法求解
        theta: GMM模型参数; param: 其它系数; eps: 计算精度; max_iter: 最大迭代次数
        返回对数似然和参数theta,theta是包含pi、mu、sigma的Python字典
    '''
    for i in range(max_iter):
        ll_old = 0
        # E-step
        q = E_step(theta, param)
        # M-step
        theta = M_step(X, q, theta, param)
        ll_new = likelihood(X, theta)
        if np.abs(ll_new - ll_old) < eps:
            break;
        else:
            ll_old = ll_new
    return ll_new, theta
# EM算法求解GMM,最大迭代次数为1e5
ll, theta2 = EM_GMM(X, theta, param, max_iter=10000)
# 由theta计算联合正态分布的概率密度
L = 100; Xlim = [-6, 6]; Ylim = [-6, 6]
XGrid, YGrid = np.meshgrid(np.linspace(Xlim[0], Xlim[1], L), np.linspace(Ylim[0], Ylim[1], L))
Xout = np.vstack([XGrid.ravel(), YGrid.ravel()]).T
MVN = np.zeros(L*L)
for i in range(param['k']):
    MVN += GMM_component(Xout, theta, i)
MVN = MVN.reshape((L, L))
# 绘制结果
plt.plot(X[:, 0], X[:, 1], 'x', c='gray', zorder=1)
plt.contour(XGrid, YGrid, MVN, 5, colors=('k',), linewidths=(2,))

最大期望算法改进算法

基于贝叶斯推断(Bayesian inference)的EM算法
在MLE理论下,EM算法仅能给出模型参数
的单点估计,引入贝叶斯推断方法后,EM算法能够给出模型参数的后验(posterior)分布避免过度拟合,其中常见的例子是极大后验估计(Maximum A Posteriori estimation, MAP)的EM算法 [2]  [18]  。MAP-EM在标准EM算法的基础上引入了模型参数的先验(prior):
,此时MAP-EM的优化目标由模型的似然转变为后验,其离散形式可表示为 [18] 
类比标准EM算法,考虑隐分布
后,由Jensen不等式可得到对数后验的代理函数,即隐变量的自由能 [18] 
由此可得MAP-EM的迭代步骤:
MAP-EM在Dempster et al. (1977)就已被提出 [1]  ,但不同于标准EM,MAP-EM的隐分布
是隐变量和模型参数的联合分布,其对应的隐变量后验
往往没有解析形式。在贝叶斯体系下,求解该隐变量后验的方法包括马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)和变分贝叶斯估计(Variational Bayesian Inference, VBI),对前者,可证明由MCMC求解的MAP-EM等价于吉布斯采样(Gibbs sampling)算法 [19]  ;对后者,由VBI求解的MAP-EM被称为变分贝叶斯EM算法(Variational Bayesian EM algorithm, VBEM) [18]  [3] 
VBEM使用平均场理论(Mean Field Theory, MFT)将隐分布近似为其在各个维度上分布的乘积:
并由此得到以下迭代步骤 [2] 
使用VBEM的常见例子是语言建模问题中的隐含狄利克雷分布(latent dirichlet allocation) [3] 
EM梯度算法(EM gradient algorithm)和广义EM算法(Generalized EM algorithm, GEM)
EM算法的M步通过计算偏导数求解代理函数的极大值,EM梯度算法(EM gradient algorithm)将该过程替换为牛顿迭代法(Newton-Raphson method)以加速迭代收敛 [2]  [20]  。更进一步地,当代理函数
不是凸函数或无法有效地对
求解极大值时,可以使用广义EM算法(GEM)。GEM有两种实现方式,一是在M步使用非线性优化策略,例如共轭梯度算法(conjugate gradients algorithm) [21]  ,二是将原M步的求导计算分解为多个条件极值问题逐个计算模型参数,后者也被称为最大条件期望算法(Expectation Conditional Maximization algorithm, ECM) [15] 
EM算法的E步也可以按ECM的方法分解为条件极值问题,由先前推导可知,E步的优化问题仅有一个全局极大值,即隐分布
,因此在E步将MLE的优化目标:联合似然
对观测样本按因子展开并对每个展开分别使用EM算法,可以得到同样的优化结果。对于M步,如果观测样本来自指数族分布,则M步也可以在每次迭代仅对有限个样本的展开进行。在指数族问题中使用EM算法的上述推广,可以避免在迭代中反复处理整个观测样本集,降低计算开销 [15]  [22] 
α-EM算法(α-EM algorithm)
α-EM算法是对标准算法的隐变量概率分布引入权重系数
的改进版本。标准的EM算法是α-EM算法在
时的特例。给定恰当的超参数
,α-EM能够比标准EM算法更快收敛。有研究将α-EM算法应用于神经网络的概率学习和隐马尔可夫模型的参数估计 [23-24] 

最大期望算法性质

编辑
收敛性与稳定性:EM算法必然收敛于对数似然的局部极大值或鞍点(saddle point),其证明考虑如下不等关系 [3] 
由上式可知EM算法得到的对数似然是单调递增的,即从
次迭代到
次迭代,EM算法至少能维持当前的优化结果,不会向极大值的相反方向运动,因此EM算法具有数值稳定性(numerical stablility)。上述不等关系也被用于EM算法迭代终止的判定,给定计算精度
,当
时迭代结束。EM算法收敛性的具体证明参见Wu (1983) [4] 
计算复杂度:在E步具有解析形式时,EM算法是一个计算复杂度和存储开销都很低的算法,可以在很小的计算资源下完成计算 [2]  。在E步不具有解析形式或使用MAP-EM时,EM算法需要结合其它数值方法,例如变分贝叶斯估计或MCMC对隐变量的后验分布进行估计,此时的计算开销取决于问题本身 [2]  [3] 
与其它算法的比较:相比于梯度算法,例如牛顿迭代法和随机梯度下降(Stochastic Gradient Descent, SGD),EM算法的优势是其求解框架可以加入求解目标的额外约束,例如在高斯混合模型的例子中,EM算法在求解协方差时可以确保每次迭代的结果都是正定矩阵 [3]  。EM算法的不足在于其会陷入局部最优,在高维数据的问题中,局部最优和全局最优可能有很大差异 [2] 

最大期望算法应用

编辑
EM算法及其改进版本被用于机器学习算法的参数求解,常见的例子包括高斯混合模型(Gaussian Mixture Model, GMM) [5]  、概率主成份分析(probabilistic Principal Component Analysis) [25] 隐马尔可夫模型(Hidden Markov Model, HMM) [6]  等非监督学习算法。EM算法可以给出隐变量,即缺失数据的后验
,因此在缺失数据问题(incomplete-data probelm)中有应用 [1]  [2] 
参考资料
  • 1.    Dempster, A.P., Laird, N.M. and Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. Journal of the royal statistical society. Series B (methodological), pp.1-38.
  • 2.    McLachlan, G. and Krishnan, T..The EM algorithm and extensions (Vol. 382):John Wiley & Sons.,2007:pp.1-73, 276-277
  • 3.    Polykovskiy, D. and Novikov, A., Bayesian Methods for Machine Learning  .Coursera and National Research University Higher School of Economics.2017[引用日期2018-12-12]
  • 4.    Wu, C.J., 1983. On the convergence properties of the EM algorithm. The Annals of statistics, pp.95-103.
  • 5.    Xu, L. and Jordan, M.I., 1996. On convergence properties of the EM algorithm for Gaussian mixtures. Neural computation, 8(1), pp.129-151.
  • 6.    Bilmes, J.A., 1998. A gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models. International Computer Science Institute, 4(510), p.126.
  • 7.    Newcomb, S., 1886. A generalized theory of the combination of observations so as to obtain the best result. American journal of Mathematics, pp.343-366.
  • 8.    Dietz, K., 1997. Introduction to McKendrick (1926) Applications of mathematics to medical problems. In Breakthroughs in Statistics (pp. 17-57). Springer, New York, NY.
  • 9.    Healy, M. and Westmacott, M., 1956. Missing values in experiments analysed on automatic computers. Applied statistics, pp.203-206.
  • 10.    Blight, B.J.N., 1970. Estimation from a censored sample for the exponential family. Biometrika, 57(2), pp.389-395.
  • 11.    Sundberg, R., 1972. Maximum Likelihood Theory and Applications for Distributions Generated when Observing a Function an Exponential Family Variable. Institute of Mathematical Statics, Stockholm University.
  • 12.    Sundberg, R., 1974. Maximum likelihood theory for incomplete data from an exponential family. Scandinavian Journal of Statistics, pp.49-58.
  • 13.    Baum, L.E. and Petrie, T., 1966. Statistical inference for probabilistic functions of finite state Markov chains. The annals of mathematical statistics, 37(6), pp.1554-1563.
  • 14.    Baum, L.E. and Eagon, J.A., 1967. An inequality with applications to statistical estimation for probabilistic functions of Markov processes and to a model for ecology. Bulletin of the American Mathematical Society, 73(3), pp.360-363.
  • 15.    Bishop, C.M..Pattern recognition and machine learning. In information science and statistics.New York, NJ:Springer,2006:450-455
  • 16.    Do, C.B. and Batzoglou, S., 2008. What is the expectation maximization algorithm?. Nature biotechnology, 26(8), supplementary note, p.897.
  • 17.    Expecation Maximization, STA-663: Computational Statistics and Statistical Computing  .Duke University.2016[引用日期2018-12-13]
  • 18.    Beal, M.J., 2003. Variational algorithms for approximate Bayesian inference. Doctor of Philosophy Thesis. University of London.
  • 19.    Mitchell, T., Gibbs and Metropolis sampling (MCMC methods) and relations of Gibbs to EM. In: 10-701 Machine Learning  .School of Computer Science, Carnegie Mellon University.2011[引用日期2019-01-22]
  • 20.    Lindstrom, M.J. and Bates, D.M., 1988. Newton—Raphson and EM algorithms for linear mixed-effects models for repeated-measures data. Journal of the American Statistical Association, 83(404), pp.1014-1022.
  • 21.    Meng, X.L. and Rubin, D.B., 1993. Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80(2), pp.267-278.
  • 22.    Neal, R.M. and Hinton, G.E., 1998. A view of the EM algorithm that justifies incremental, sparse, and other variants. In Learning in graphical models (pp. 355-368). Springer, Dordrecht.
  • 23.    Matsuyama, Y., 1997. The α-EM algorithm: A block connectable generalized leaning tool for neural networks. In International Work-Conference on Artificial Neural Networks (pp. 483-492). Springer, Berlin, Heidelberg.
  • 24.    Matsuyama, Y., 2011. Hidden Markov model estimation based on alpha-EM algorithm: Discrete and continuous alpha-HMMs. In Neural Networks (IJCNN), The 2011 International Joint Conference on (pp. 808-816). IEEE.
  • 25.    Tipping, M.E. and Bishop, C.M., 1999. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3), pp.611-622.
展开全部 收起