初学者看EM算法最大的障碍往往不在数学推导而在“这台机器到底在解决什么问题转起来的逻辑是什么”。我之前学的时候先啃了很多遍公式感觉每一步都能看懂合上书一问“这套操作凭什么能把参数估出来”还是答不上来。后来换了个思路从带隐变量的数据入手把整个流程手写一遍才真正明白EM算法是什么说白了它就是一套在数据中有隐藏变量、无法直接做极大似然估计时通过反复“猜测隐变量 更新参数”来完成参数估计的通用框架。而且它适用的场景远不止面试题像高斯混合模型、隐马尔可夫模型、带缺失数据的回归底层都会见到它的影子。如果你是正在学机器学习的读者或者工作中需要处理混合模型、聚类、缺失数据这类问题这篇内容应该能帮你把EM算法从“听说过”变成“能自己写出来”。整个思路我会尽量用大白话推进每一段配合实际项目里需要注意的细节来讲不会只丢公式就完事。1. 先看一个“数据丢了标记”的小例子EM算法在解决什么问题1.1 两枚硬币实验里的隐变量不妨先想象一个特别朴素的实验我有两枚硬币A和B它们各自的正面概率分别记作θ_A和θ_B。现在每一轮我随机挑一枚硬币然后连续抛10次记录这一轮出现的正面次数。假如做了很多轮之后我手上有一个表格写着第1轮正面6次、第2轮正面8次、第3轮正面3次……全部记录都在唯独少了一列每一轮到底用的是A还是B。这个“不知道用哪枚硬币”的状态就是一个典型的隐变量用Z来表示。观测到的正面次数可以叫X是我们看得见的部分。现在的问题是在Z缺失的情况下怎么估计两个硬币的正面概率θ_A和θ_B这个例子虽然简单但它和很多实际问题在逻辑上是同构的。比如用户聚类时我们只能看到用户的购买行为但看不到他到底属于哪个消费群体做语音识别时只能观测到声学特征却不知道上一时刻处于哪个发音状态文本建模时能看到单词却看不到文章背后的主题。这些“看不见的类别”就是隐变量。1.2 如果能看见隐变量问题会变成什么样先绕一步假设数据没有缺失每一轮我清清楚楚知道用的是A还是B。那估计参数会非常轻松。统计所有用A硬币的轮次累计正面数除以累计抛掷数就是θ_A同理算B。用公式说就是最大似然估计写下观测数据的对数似然函数求导令其等于零解出来的闭合解往往就是“数数然后求比例”。问题就出在“不知道用哪枚硬币”的时候。如果不知道隐变量观测数据的边缘似然就要写成对所有可能的隐变量求和。比如某一轮正面次数是h那么这一轮数据出现的概率就是“选A且抛出了h次正面”加上“选B且抛出了h次正面”两项求和。整个实验的对数似然就变成log L(θ_A, θ_B) Σ log [ p(ZA) * C(10,h) * θ_A^h (1-θ_A)^(10-h) p(ZB) * C(10,h) * θ_B^h (1-θ_B)^(10-h) ]原本对数里是单个概率密度求导很容易现在对数里面套了一个“和的求和”每个参数项都被其他参数纠缠在一起。直接对θ求导并令其为零得到的方程不是线性的也没有一眼能看出的闭合解。这种情况下EM算法就该上场了。2. E步和M步到底怎么迭代从猜硬币到公式落地2.1 E步先猜一轮每个样本有多大概率属于某个“组”EM算法的核心思想可以概括成一句非常朴素的话既然隐变量看不见那就先用当前参数猜测它大概是什么然后把这个猜测当作已知重新估计参数估计完参数以后再回头更新猜测循环往复。第一次循环总得有个起点所以我们会先给参数瞎猜一个初值比如θ_A0.6、θ_B0.4。这里的“猜”不是让样本硬性分给A或B而是计算后验概率。还是用抛硬币例子某轮正面次数为h那么在当前参数下这一轮来自A的概率可以算出来γ_A θ_A^h (1-θ_A)^(10-h) / [ θ_A^h (1-θ_A)^(10-h) θ_B^h (1-θ_B)^(10-h) ]γ_B就是1-γ_A。这个γ就是隐变量的责任度术语叫responsibility它是EM算法里E步的核心产物。注意一个关键点我们不是直接猜测“这一轮用A”而是给这一轮数据打了一个例如“70%来自A、30%来自B”的软标签。这种软化处理比硬归类保留了更多不确定性这也是为什么EM在很多问题上比某些启发式聚类更稳。2.2 M步再估计一次把加权后的样本当成完整数据拿到每个样本的责任度后就可以用加权数据重新估计参数。假设所有轮的正面次数为h_1, h_2, …, h_m。第一次迭代第i轮属于A的责任度标记为γ_i^A那么θ_A的新估计就是θ_A_new Σ γ_i^A * h_i / (Σ γ_i^A * 10)这个式子和“知道隐变量后直接数正面数除以总次数”长得几乎一样区别只是每个样本都乘以了一个权重。如果某一轮有80%的可能来自A那它对参数估计的贡献就打八折。这一点非常关键M步其实是在做加权最大似然估计。由于每一轮的总抛掷次数相同分母直接乘10如果每轮次数不一样只需要改成对抛掷次数加权求和逻辑一致。E步和M步轮流执行每一轮迭代都会得到一组新参数。新的θ_A和θ_B会改变所有轮次的责任度责任度变化后又会影响下一轮M步的参数估计。如此反复迭代直到参数变化或者似然值变化足够小才停止。我第一次跑这种迭代的时候看着参数从0.6、0.4慢慢往0.8、0.3附近移动才真正明白教科书里说的“期望最大化”不是在解决什么高深玄学就是不断用概率猜测补齐缺失信息再重新计算。2.3 为什么这样更新能保证似然不减有个问题一定会在实践中冒出来凭什么这几行代码能保证最终收敛我最初的担心是每次M步更新都基于当前责任度算的局部最优会不会换了参数之后整体观测似然反而下降EM算法之所以可靠是因为它的每一次更新都不会让观测数据的对数似然下降。标准证明会用到Jensen不等式这里提供一个直觉我们真正想最大化的是log p(X | θ)但隐变量存在导致表达式里是“对隐变量求和以后再取对数”不好处理。EM换个思路每一步都构造一个当前θ_t条件下的对数似然下界学术上叫ELBO然后去最大化这个下界。在做完一次最大化后下界会顶到当前θ_t时对数似然的位置下一次E步会把下界重新贴紧新的后验分布M步再把下界往上抬一截周而复始。因为每一步都是“贴紧—上移—再贴紧—再上移”所以最终观测似然会被一下一下顶上去。也可以换一种工程上的理解把每轮迭代想象成爬山虽然不知道山的全貌但每一步都能保证自己往上走一小步不会往后退。只要没有数值bug这个保证是严格的。也正因如此如果你在某次EM迭代后看到似然函数显著下降基本可以断定代码里出了实现问题或者数据里有NaN之类的脏东西。2.4 另一种好用的视角EM是在做坐标上升接触贝叶斯方向以后我越发觉得EM还可以理解为在坐标上升中交替更新两部分变量一部分是隐变量对应分布的参数另一部分是模型参数。如果我们给目标函数加上一个KL散度项构建出ELBO那么每次E步是在固定模型参数的情况下优化关于隐变量分布的这部分每次M步是在固定隐变量分布的情况下优化模型参数。这就是一个标准的坐标上升过程。这个视角不只是为了看得更“高维”它其实帮我们理解了很多变体算法。比如有些复杂模型里E步无法求得精确后验那就用变分推断去近似这一步于是就有了变分EM。如果M步也无法求出解析解可以用数值优化做一步上升就有了广义EM(GEM)。先把这些名字放下理解核心是“两块变量交替更新每一步都往不下降的方向走”后面再接触高级变体就不会懵。3. 最经典落地场景手写一个高斯混合模型的EM实现3.1 高斯混合模型要估计的参数如果说硬要选一个“EM算法练手最好的项目”我一定会推荐高斯混合模型。GMM假设数据由K个高斯分布混合生成每个样本先按概率π_k选择一个组件再从第k个高斯分布里采样。现在我只观测到一堆样本点不知道它们具体来自哪个组件这个组别就是隐变量。估计的目标参数包括三组每个组件的混合权重π_k、均值向量μ_k、协方差矩阵Σ_k。和高斯分布的直观区别是GMM可以拟合一个分布里有多个“团”的数据而单个高斯很难描述多峰结构。实践里也常拿它做聚类K-Means给的是硬分类每个点只能归属到最近的簇GMM给的是软分类每个点分别有概率属于各个簇。这个特性在需要输出不确定性置信度的场景特别有用。3.2 E步的完整计算过程假设现在模型参数有一套初始值我们要计算每个样本x_i属于第k个高斯组件的后验概率也就是责任度γ_ik。公式是γ_ik π_k * N(x_i | μ_k, Σ_k) / (Σ_j π_j * N(x_i | μ_j, Σ_j))其中N(x | μ, Σ)是多维高斯分布的概率密度函数。这个公式很直观分子是当前参数下“样本从第k个组件来的概率”分母是对所有组件的总概率做归一化保证每个样本对所有责任度的和为1。工程上有一个立刻会踩到的坑如果直接用概率密度函数计算当数据维度较高或者某个组件离样本很远时N(x_i | μ_k, Σ_k)会非常小甚至下溢到0。分母一旦出现一堆0后面算比值就会出现NaN。解决办法是全程在log空间里计算log密度也就是用对数概率相加再通过log-sum-exp技巧做归一化。下面这段是我经常用的一种实现骨架import numpy as np from scipy.stats import multivariate_normal from scipy.special import logsumexp def e_step(X, pi, mu, Sigma): n, d X.shape k len(pi) log_phi np.zeros((n, k)) for j in range(k): # 算log密度用scipy自带函数即可 log_phi[:, j] np.log(pi[j]) multivariate_normal.logpdf(X, meanmu[j], covSigma[j]) # 用log-sum-exp归一化顺便算观测对数似然 log_lik logsumexp(log_phi, axis1).sum() gamma np.exp(log_phi - logsumexp(log_phi, axis1, keepdimsTrue)) return gamma, log_lik这里有两个值得讲清楚的点。第一对每个组件logpdf时均值向量和协方差矩阵会直接决定这一项的数值计算前一定要确保协方差矩阵是合法的对称正定矩阵。第二归一化时我没直接在原始空间除以分母而是先减掉每行的log-sum-exp得到的responsibility数值会稳定得多。很多人第一次自己写GMM跑不出正确结果八成就是在这一步丢了精度。3.3 M步更新参数的三件事E步结束之后手里有每个样本对每个组件的责任度γ_ik。M步要做的事情就轻松了用这些责任度作为权重重新计算三组参数。第k个组件的混合权重π_k直接是所有样本对第k个组件的平均责任度π_k (1/N) * Σ γ_ik均值向量μ_k是按责任度加权的样本均值μ_k (Σ γ_ik * x_i) / (Σ γ_ik)协方差矩阵Σ_k是按责任度加权的样本协方差Σ_k (Σ γ_ik * (x_i - μ_k)(x_i - μ_k)^T) / (Σ γ_ik)这三个公式的共同点是把责任度当作这个样本“属于第k类”的权重。所有贡献做完加权求和后再用权重的总和做归一化。如果第k个组件的总责任度过低说明这个组件几乎没有样本支持更新出来会很不稳定这也是后面的坑点。M步更新完之后用新参数回到E步继续循环直到观测对数似然的变化量小于给定阈值。3.4 每次迭代都要盯着对数似然写GMM的时候一定要保留一个记录把每一步的观测对数似然log p(X | θ)打印出来。它不仅仅是监控收敛的指标也是检查代码有没有写错的试金石。理论上EM算法每一步更新后似然只增不减一旦你在打印日志里发现下降比如从-1200跳到-1350那就说明某个责任度计算错了或协方差矩阵更新时用了上一轮的均值而不是本轮的均值。有些实现还会额外检查参数更新的幅度比如当前后两次参数变化小于1e-6就提前终止。但在实际数据上我更推荐看对数似然的增量因为参数空间不同维度取值范围差异很大直接用参数差做阈值可能过早或过晚停止。我常用的停止条件是连续两次对数似然差小于1e-6或者达到最大迭代次数。4. 初值与收敛性几个让我折腾很久的坑4.1 EM算法只能保证找到局部最优EM算法有一个很扎心的现实它只能保证收敛到一个局部最优解不是全局最优。原因是GMM对参数的似然函数通常不是凸函数不同初值走出来的结果可能差很多。我做过一次测试同一个数据集用不同随机初值跑GMM有的收敛到两个聚类中心恰好分开有的收敛到两个中心重合、输出一个看似正常但明显没分开数据的模型后者的对数似然也更低。这在实际项目里意味着不能指望一次随机初始化就得出好模型。一个很常见的解决方案是“多次重启”用不同初值跑十次、二十次每次保留最终对数似然值最高的那次作为最终模型。多次重启看起来笨但在EM这种迭代算法里非常实用。除了随机初始化更稳的做法是先用K-Means聚类的结果作为初始中心。K-Means本身也是一种迭代但相对快速的硬聚类方法它的聚类中心比随机点更接近数据真实分布用这个输出当GMM的均值初始值通常能让EM收敛得又快又稳。4.2 K-Means初始化的具体配方具体操作可以这样先用K-Means把样本分成K类记录每类的中心点center_k和标签。GMM启动时的μ_k直接用center_kπ_k用每类样本数量除以总样本数Σ_k则用同一类内样本的协方差矩阵对样本极少的情况可以加一个很小的正则项防止协方差阵奇异。如果连K-Means都不想调也可以直接随机抽K个样本点作为初始μ但效果稍差尤其数据维度高时可能把大量组件初始化到空白区域导致收敛后出现空组件。多次随机启动配合其他初始化手段通常能规避大部分局部最优。实际项目中我会这样组合前几次启动用K-Means初始化如果发现了明显更好的似然值就继续跑没有的话再用随机初始化做几次探索。因为随机初始化覆盖面更广有可能找到K-Means初始化碰不到的更优解。4.3 协方差矩阵奇异问题GMM在实践里最常见的崩溃点是迭代过程中某个组件只获得了极少样本或者某些样本被某个组件以接近1的概率抓住导致这个组件的Σ_k趋近于奇异矩阵。高维数据下尤其频繁因为样本数量不足时协方差矩阵的估计天生就不稳定。出现奇异矩阵后计算multivariate_normal.logpdf会直接报错或者返回无穷大。我的处理方式是在协方差矩阵的对角线上加一个很小的扰动比如1e-6这一步通常叫jitter或regularization。如果需要更平滑还可以用每个维度的全局方差乘以0.01作为下限。另一个做法是设置一个“最小责任度总计数”当某个组件对全部样本的责任度总和小于某个极小阈值时跳过该组件更新并重新初始化它。这些手段是在“拟合精度”和“数值稳定”之间做取舍加了正则后模型往往不容易过拟合在验证集上也更稳。4.4 什么时候能判断EM“没跑好”判断GMM拟合质量不能只看收敛时的对数似然。我习惯再检查每个组件的责任度分布如果一个组件的平均责任度特别低等于说模型里混进了一个几乎不承担职责的组件。这种组件不会贡献有效聚类应该删掉或合并。还可以用轮廓系数之类的聚类指标辅助判断但要注意GMM是软聚类模型直接把硬标签拿出来做轮廓系数会损失一部分概率信息。5. EM算法常见问题排查速查表与我的避坑经验5.1 一张表解决高频问题先整理一张速查表把EM/GMM实战中最容易遇到的问题、可能原因和处理办法列一下问题现象可能原因处理办法迭代中似然值下降E步用错上一轮参数或M步用错更新后均值重新检查责任度与参数更新的先后顺序协方差矩阵奇异某个组件样本太少或数据维度太高加对角正则项限制最小方差多个组件收敛到同一个位置初始值太差EM陷入对称局部最优换K-Means初始化多次重启责任度全是NaN密度计算下溢或协方差非法全程log空间计算用logsumexp归一化收敛非常慢初值差、组件K过多、数据维度高用K-Means初始化适当增加正则收敛模型聚类结果明显不合理局部最优或K大小不合理用BIC/AIC选K多做几次随机重启这张表里的每一行都不是编出来的至少前四类问题我在写自己的EM代码时都见过。尤其是责任度出现NaN那次我查了半天最后发现log密度里有负无穷减去负无穷一归一化就全成了NaN。从那以后我的所有概率计算都强制在log空间完成如果必须回原始空间也一定先做log-sum-exp归一化。5.2 常规文档里不会写的两个操作细节第一个细节协方差矩阵更新后要做对称化。理论上加权外积求和出来的矩阵一定是对称的但在浮点数运算下因为有舍入误差Σ_ij和Σ_ji可能出现微小差异。这个差异平时无所谓但后续做Cholesky分解或逆矩阵时可能会触发告警。所以每次更新完协方差最好执行一行M (M M.T) / 2把它重新变成准确的对称矩阵。第二个细节打印对数似然时不要直接算p(x)而是用logsumexp直接得到每行的log概率和。我见过不少人用概率密度乘权重求和后取对数在高维情况下数组里可能全是0或NaN输出结果就废了。优先在log密度矩阵上做计算最后按行加总。很多成熟的库内部都在使用这类技巧目的只有一个避免中间过程的数值下溢。5.3 关于“K选多少”的一点个人建议GMM需要事先指定组件数量K这个K和聚类数量含义类似。最简单的方法是用信息准则对若干候选K逐个训练模型分别计算BICBIC值越低说明模型在拟合复杂度和样本解释之间越平衡。我自己通常跑K从1到10的序列每个K做几次随机启动取最优然后画一张BIC曲线。曲线会有一个明显的拐点或最低点这个K一般就是比较合理的选择。我踩过一个坑为了追求每个子簇更“精细”把K设得很大结果模型把正常分布切成了很多窄片段聚类结果完全没有可解释性。所以选定K后还要结合业务去看每个组件的主特征如果两个组件的中心几乎一样只是方差被拆开那K大概率设大了。6. 从EM延伸开硬EM、变分EM与我的使用建议6.1 K-Means其实可以看成硬EM理解了GMM之后回头看K-Means会发现两者有很强的联系。如果把GMM中每个高斯的协方差固定为σ²I并且让σ²趋于0那么责任度计算会趋于把每个样本硬分配给离它最近的那个中心M步就等价于把每个簇的中心移到样本均值上。所以K-Means其实是GMM在“协方差各向同性、且极限小”情况下的一种简化也叫硬EM。这个关系对我理解两者的行为很有帮助。K-Means跑的聚类结果通常比GMM“更干脆”但缺少不确定性度量GMM则能给出每个样本属于各个组的概率分布还允许不同簇有不同形状的协方差结构。实际选哪个取决于需求。如果下游业务只在意最终分组标签K-Means更快如果后续流程需要把概率传给另一个模型或者数据各簇的宽度差异很大GMM更合适。6.2 当后验算不出来时就进入变分EM的范畴标准EM有一个前提E步需要计算隐变量的后验分布p(Z | X, θ)在GMM这类共轭结构里它是解析可求的。但很多现代模型没有这么幸运比如主题模型、深度生成模型中的隐变量后验分布形式很复杂无法直接算。那怎么办退一步用变分推断限定隐变量后验在一个参数化的简单分布族q(Z)里通过最小化q与真实后验的KL散度来近似E步。这就变成了变分EM。很多人在看“变分自编码器”或其他生成模型时觉得很难换个角度想它也就是在EM框架下用神经网络去拟合那个无法写出闭合解的E步。先理解EM再理解“E步被换成了近似计算”后面接触这些复杂模型会顺畅得多。6.3 我建议的“要不要用EM”判断标准做了几年实际项目后我形成了一个简单的判断标准如果某个模型在补全隐变量后参数的极大似然解非常简单那就适合用EM反之如果连补全隐变量后的参数解都很麻烦可能更适合直接用梯度优化整体目标函数。EM的优势在于可以复用那些成熟、稳定的加权更新公式而且不需要手动调学习率。但当目标函数复杂、隐变量结构又很灵活时EM的E步和M步很容易卡在推导和数值计算上这时用自动求导框架对整个ELBO做优化反而省力。听上去像是偷懒但工程上确实是“哪个好用用哪个不必对某个算法有执念”。我见过不少团队在模型里死磕EM推导最后发现用基于梯度的优化方法两小时就调通了效果还更好。了解一种算法的边界和替代方案往往比把它的每一步推导都刻进脑子里更重要。EM算法本身的价值在于它提供了一个经典思路未知量太多时先猜再估反复逼近。很多模型即使最终不用它也仍然继承了这个循环迭代的思想。