期望最大化算法 (Expectation-Maximization Algorithm)
期望最大化算法(Expectation-Maximization Algorithm,简称 EM 算法)是一种通过迭代方式求解含有隐变量(Latent Variables)或缺失数据的概率模型极大似然估计(MLE)和最大后验估计(MAP)的通用方法。由 Arthur Dempster、Nan Laird 和 Donald Rubin 于 1977 年在《皇家统计学会杂志》上发表的经典论文《通过 EM 算法从不完整数据中求极大似然估计》中系统提出并命名。事实上,EM 的核心思想早已以特例形式存在于多个领域:Baum-Welch 算法(隐马尔可夫模型,1970)、McKendrick 的基因频率估计(1926)乃至 Newcomb 在 1886 年对缺失数据问题的处理均可视为 EM 的雏形。Dempster 等人的贡献在于统一了这些分散的特例,给出了收敛性的严格证明,并正式命名为 EM。如今 EM 算法已广泛应用于聚类分析、混合模型、隐马尔可夫模型、因子分析、项目反应理论和缺失数据处理等计算统计学领域。
问题设定与直观动机
给定观测数据 X={x1,x2,…,xN} 和未观测的隐变量 Z={z1,z2,…,zN},目标是极大化关于参数 θ 的观测数据的对数似然:
ℓ(θ;X)=logp(X∣θ)=log∫p(X,Z∣θ)dZ
由于积分(或求和)出现在对数内部,直接优化上式通常难以得到解析解,而数值优化又面临高维积分的计算困难。然而,若隐变量 Z 已知,完整数据的对数似然 logp(X,Z∣θ) 往往属于指数族分布,形式简洁且易于优化——这正是 EM 算法的出发点:既然真实隐变量未知,就以迭代方式利用当前参数推断隐变量的条件分布,用概率权重填补缺失信息,然后在此"完整数据"假设下优化参数。这种"先猜测、再优化"的迭代策略使得原本不可解的极大化问题变得可计算。
算法步骤
EM 算法每次迭代包含两个步骤:
E 步:期望步 (Expectation Step)
基于当前参数估计值 θ(t),计算隐变量 Z 的后验分布 p(Z∣X,θ(t)),并构造完整数据对数似然关于此后验的期望,即定义 Q 函数:
Q(θ∣θ(t))=EZ∣X,θ(t)[logp(X,Z∣θ)]
这一步本质是"软分配"——用当前参数推断每个隐变量的归属概率,而非硬性指派。在指数族模型中,E 步只需计算隐变量的充分统计量的条件期望,计算复杂度通常可控。
M 步:最大化步 (Maximization Step)
极大化 Q 函数以获得新的参数估计:
θ(t+1)=argθmaxQ(θ∣θ(t))
若 Q 函数的极大化有解析解(如高斯混合模型),则直接更新;否则可借助数值优化方法。当 M 步仅使 Q 上升而非完全极大化时,称为广义EM算法(Generalized EM, GEM),同样保证似然单调性。
交替执行 E 步和 M 步直至 θ 或对数似然收敛于预设阈值。
理论保证:Jensen 不等式推导
EM 算法的单调性可由Jensen 不等式严格证明。对于任意隐变量的概率密度 q(Z),观测对数似然可分解为:
logp(X∣θ)=log∫p(X,Z∣θ)dZ=log∫q(Z)q(Z)p(X,Z∣θ)dZ≥∫q(Z)logq(Z)p(X,Z∣θ)dZ≡L(q,θ)
上述不等号来自 Jensen 不等式(因 log 是严格凹函数),L(q,θ) 称为证据下界(Evidence Lower Bound, ELBO)。对数似然与 ELBO 之间的差距恰好是 q(Z) 与真实后验 p(Z∣X,θ) 之间的KL 散度:
logp(X∣θ)−L(q,θ)=DKL(q(Z)∥p(Z∣X,θ))≥0
EM 算法的 E 步在固定 θ(t) 时,令 q(Z)=p(Z∣X,θ(t)) 使 KL 散度归零,ELBO 紧贴对数似然;随后的 M 步固定此 q 极大化 ELBO 关于 θ,从而保证原始对数似然单调不减:
ℓ(θ(t+1);X)≥ℓ(θ(t);X)
这一单调性保证了 EM 的数值稳定性,是其相对于直接使用梯度下降方法的重要优势——每次迭代至少不使目标函数变差。
经典应用:高斯混合模型
高斯混合模型(Gaussian Mixture Model, GMM)是 EM 算法的典型教学案例,也是聚类分析的核心工具。假设数据由 K 个高斯分量混合生成:
p(xi∣θ)=k=1∑KπkN(xi∣μk,Σk),k=1∑Kπk=1
其中隐变量 zik∈{0,1} 指示第 i 个样本属于第 k 个分量。完整数据似然为:
p(X,Z∣θ)=i=1∏Nk=1∏K[πkN(xi∣μk,Σk)]zik
E 步: 计算责任度(responsibility),即样本 i 来自分量 k 的后验概率:
γik(t)=∑j=1Kπj(t)N(xi∣μj(t),Σj(t))πk(t)N(xi∣μk(t),Σk(t))
M 步: 以责任度为权重更新参数:
πk(t+1)=N1i=1∑Nγik(t),μk(t+1)=∑i=1Nγik(t)∑i=1Nγik(t)xi
Σk(t+1)=∑i=1Nγik(t)∑i=1Nγik(t)(xi−μk(t+1))(xi−μk(t+1))⊤
直观上,M 步的 μk 是所有数据点的加权平均,权重即为责任度——样本越可能属于分量 k,对更新 μk 的贡献越大。这也解释了为何 EM 被称为"软 K-means":K-means 聚类是 GMM 在协方差矩阵趋近于 ϵI(ϵ→0)、各分量等权、责任度退化为最近邻硬分配时的极限特例。相比于 K-means 的硬聚类,GMM 的软分配能表达归属的不确定性,并提供聚类的概率解释和模型选择准则(如BIC)。
收敛性与局限性
收敛性: EM 确保对数似然单调不减,这是其最重要的实用优势。收敛速率取决于缺失信息比例——Wu(1983)指出,当缺失信息量较大时,EM 仅以线性速率收敛,远慢于牛顿法的二次收敛。Murray(1977)进一步证明,若将 EM 视为一种不动点迭代,其收敛速率由完整信息矩阵与缺失信息矩阵之比的谱半径决定。
主要局限:
- 局部最优: 似然函数非凹时,EM 仅保证收敛到局部极大值或鞍点,结果高度依赖初始值。实践中常运行多次随机初始化,选取对数似然最高者作为最终解。确定性退火(deterministic annealing)等启发式策略也可改善全局搜索能力。
- 慢收敛: 在高维或信息大量缺失场景下,线性收敛速率可能极慢。此时可考虑使用加速EM(如 Aitken 加速)或在收敛后期切换至 Newton-Raphson 方法。
- E 步不可解: 当后验 p(Z∣X,θ) 无解析形式时,需借助蒙特卡洛EM(MCEM)采样近似期望,或使用变分推断(Variational Inference)以可处理的分布族逼近真实后验。
- 奇异解: 在 GMM 等模型中,单一分量可能坍缩到单个数据点,导致协方差矩阵奇异(似然发散至无穷大)。实践中加入协方差矩阵的正则化项、使用共轭先验的 MAP-EM 框架,或在 M 步约束协方差的最小特征值,均可有效缓解此问题。
变体与推广
EM 算法衍生出多种变体以应对不同场景:ECM(Expectation Conditional Maximization)在 M 步中将参数分块条件优化,将高维极大化分解为一系列低维子问题,显著降低每次 M 步的计算复杂度;ECME 进一步允许某些条件最大化步直接优化原始对数似然而非 Q 函数,以牺牲部分理论简洁性换取更快的收敛速度;MCEM 在 E 步积分不可解析时使用 Markov Chain Monte Carlo(MCMC)采样近似期望,适用于贝叶斯后验推断场景;变分EM(Variational EM)在 E 步使用一族可处理的分布 q(Z∣ϕ) 近似真实后验,通过优化 ELBO 而非精确求解后验来回避计算瓶颈,广泛应用于主题模型(如隐含狄利克雷分布 LDA)和贝叶斯深度生成模型(如变分自编码器 VAE)。此外,随机EM(Stochastic EM, SEM)在 E 步中对每个隐变量进行随机采样而非计算期望,通过多次随机模拟的均值获得参数估计,在小样本和离散隐变量场景中表现稳健。