如果你翻到课程目录的第九章看到EM算法四个字多半已经做了点心理准备。我在读研那会儿第一次看教材里那几页公式脑子里的画面是一个没见过数据的老师非要让一个蒙着眼睛的学生在教室里找黑板擦。后来工作里真正处理带缺失值的数据才意识到EM其实是个非常朴素的策略先拿一个不靠谱的初始参数去猜隐藏信息再用猜出来的信息反过来优化参数重复到满意为止。它不是什么高深魔法而是一套把缺信息这个问题切成两步的通用框架。这一篇我把整章脉络拆开讲清楚它为什么这么设计、公式每一步在算什么再用一个高斯混合模型的小实验把过程跑通顺便把我在实际建模时踩过的坑也写出来。1. EM算法先搞清楚它是来解决什么问题的1.1 一个看似很简单的估计场景假设你收到一万条身高数据想用两个高斯分布去拟合代表男性和女性的身高分布。理想情况下如果每条数据都标记了性别求两个分布的均值和方差就是小菜一碟对男性样本求均值和方差对女性样本求均值和方差收工。这个操作谁都会。难点在于性别字段没记录或者你研究的东西本身就不能直接观测。这时候你面对的数据是X想推测的类别Z却不可见。如果直接对总似然函数做极大化会得到一个带 log 内嵌累加的式子。以两个高斯分量为例目标函数长这样log ℓ(θ) Σᵢ log [ π₁ N(xᵢ|μ₁,Σ₁) π₂ N(xᵢ|μ₂,Σ₂) ]问题在于 log 无法穿进加法号两个子分布的参数缠在一起求导之后没有漂亮的闭式解。如果Z是可见的情况完全不同样本按Z分组后每组内部的似然可以独立分解梯度也自然拆开。一旦Z不可见所有子分布的项都搅合在一个和式里。EM算法本质上就是给这个“纠缠”解套的思路。1.2 把“一步到位”改成“两步循环”既然直接求最优θ很难那就换一种玩法。先随便猜一组初始参数θ₀这组参数可以很糟糕但只要有一点点区分能力就能回答“每条样本更像哪个分布”这个问题。基于θ₀给每条数据算一个软标签它属于男生分布的概率是0.7属于女生分布的概率是0.3。有了软标签之后再回头更新参数。这里有个关键点更新参数时不是按硬分类结果去分组而是按概率加权。第i个样本对男生组的贡献是0.7对女生组的贡献是0.3两边都用到这条数据只是权重不同。E步和M步交替循环每轮都让完整数据对数似然的某个下界抬升。用大白话说E步是在猜身份M步是在依据身份更新群体画像。这个过程很像认识新朋友先凭第一印象判断对方大概属于哪类人再根据一群人的共性修正印象然后再重新判断循环几轮之后画像会越来越准。1.3 EM不是无中生有它是一套通用的“补洞”框架不少同学学完EM会以为这算法就是为高斯混合模型服务的这种理解太窄了。任何参数模型只要似然能够写成“观测变量隐藏变量”的联合概率积分或求和并且隐藏变量一旦补齐模型就变成容易优化的问题就可以套用EM。隐马尔可夫模型里的前向-后向训练算法本质就是EMLDA主题模型里的变分推断走的也是类似路径甚至一些缺失数据填补、推荐系统里的隐因子模型都能落到完全相同的大框架里。所以学EM时我建议大家先吸收建模习惯判断哪个变量在“偷懒”不肯露面再把推断分解成E步和M步。只要把模型化成这种完整数据形式后半段往往只是套用标准求解器的问题。把这一层想明白第九章的学习曲线会突然变平。2. 核心公式拆解E步和M步到底在算什么2.1 从“对数似然拆不开”说起公式还是要看一点但可以看得轻松一些。设观测数据为X参数为θ隐变量为Z。我们希望最大化的是ℓ(θ) log p(X|θ) log Σ_Z p(X,Z|θ)麻烦就麻烦在log不能直接钻进求和号里。但如果拥有完整数据(X,Z)要优化的目标变成 log p(X,Z|θ)。因为联合分布通常会按变量关系分解成若干项的乘积取对数之后就会变成一堆可独立处理的项相加极大化难度显著降低。于是问题被分到两层。log 外与求和号之间有一堵墙但如果我们引入一个分布q(Z)做加权再利用Jensen不等式可以得到log p(X|θ) ≥ Σ_Z q(Z) log [ p(X,Z|θ) / q(Z) ]右侧是一个下界。如果把q选成基于旧参数θ_old的后验分布 p(Z|X,θ_old)那么在θ_old 这一点下界与原似然函数相等。接下来固定q不动去寻找能让下界变得更大的新参数θ。E步就是在做这件事拿着旧参数θ_old把隐藏信息的后验分布算出来。M步则是让q保持不变去最大化这个下界。因为每一步都让原似然的下界上升而原似然本身有限最终会收敛到一个局部极大值。初次接触的话别强迫自己从下界公式直接跳到底层代码。先记住“E步给权、M步做加权极大似然”这个宏观流程公式更重要的作用是解释这个流程为什么不跑偏。2.2 E步里的“期望”到底是对谁求的E步中经常出现这样的Q函数Q(θ; θ_old) E_{Z|X, θ_old} [ log p(X,Z|θ) ]这行式子很容易把心理素质差的人劝退。翻译成人话就是把完整数据的对数似然看成Z的函数这里Z服从后验分布 p(Z|X, θ_old)然后对Z求期望。如果Z是离散的所谓的期望就是一个按后验概率加权的和。每个样本的Z不是一个确定值而是一组概率也就是软标签。所以E步真正输出的通常不是把一个样本硬塞给某个簇而是一张N行K列的权重表。例如第i个样本属于第k个簇的权重是0.82属于另一个簇的权重是0.18。这一点是和K-means最大的差异。K-means在分配样本时只保留唯一的归属EM则保留不确定性。正是这种“软”属性让GMM在类别边界不清晰的数据上通常会比K-means稳定得多。2.3 M步为什么往往“白给”把Q函数写出来以后M步反而显得很平凡。因为log p(X,Z|θ)展开后许多项只与某个子参数相关不相关的部分求导后会被直接消掉。拿高斯混合模型来说Q函数按类别展开后均值参数μ_k的更新会变成一个典型的带权平均问题μ_k_new (Σᵢ rᵢₖ xᵢ) / (Σᵢ rᵢₖ)其中rᵢₖ就是E步算出的第i个样本属于第k个高斯的责任权重。分母是第k个聚类总权重分子是样本加权求和这正是加权版本的极大似然闭式解。方差矩阵和混合系数也类似全部有显式公式。所以完全可以把各种现成极大似然求解器装进M步口袋里使用。如果M步没有闭式解还可以用一次梯度上升去近似最大化这类做法叫广义EM实际项目中用得非常普遍。不要被M步的符号吓到它本质上就是做一次加权拟合。2.4 收敛性只需要一个直觉严格数学证明容易劝退人但EM的收敛逻辑可以靠直觉掌握。E步和M步每一步都在优化同一个下界函数。当参数位于θ_old时E步把下界函数挪到与原似然相切的位置M步沿着下界往上爬。由于下界在旧点是贴着原函数的所以爬升高下界后原似然只能跟着往上走不会掉下去。这个循环就像登山时先架一座云梯再顺势往上走一段反复架桥、前进总能到达某个山顶。不过如果山体有很多峰最终到的是不是最高峰取决于初始位置。所以我在实际项目中几乎默认使用多组随机初始化最后挑似然最大的结果。3. 实操手写一个GMM聚类把流程跑通3.1 先选一个看得见摸得着的例子教科书一般用二维高斯混合模型来演示EM因为能画图、能肉眼判断聚类效果比一上来就在高维数据上折腾直观得多。我这里构造一个包含三个分量的二维数据每个分量各300个样本分量之间有一定重叠再用EM估计三个高斯的均值、协方差和混合系数。这类问题在无监督聚类、异常检测、图像分割里有大量真实投影。3.2 代码前的三个准备工作手写EM之前先做三件事确定分量数K、初始化参数、定义好数据结构。K需要提前指定在真实场景里可以用BIC或轮廓系数去选在小实验里我直接设K3。初始化方式我推荐先用K-means跑几轮把样本粗略分成K组再以每组均值作为高斯均值μ以组内样本的协方差作为初始协方差Σ以组内样本占比作为混合系数α。这样初始化比完全随机初始化稳定得多。随机初始化如果不配合多次启动很容易让两个高斯分量在初始阶段抢同一批数据最终收敛到很差的结果。3.3 E步实现先算权重别忘归一化E步的核心是计算责任矩阵 r[i,k]也就是第i个样本属于第k个高斯的概率。代码上可以用 SciPy 的 multivariate_normal 计算多元高斯密度也可以用 NumPy 手写公式。先算score[i,k] αₖ × pdfₖ(xᵢ)再把每一行除以该行所有score之和。这里必须提醒一个新手的经典错误直接用样本到均值向量的欧氏距离作为归属依据忽略了协方差矩阵的形态和混合系数的先验。那样做等于把GMM降级成了K-means会丢掉软聚类的核心信息。3.4 M步实现三个更新公式各司其职责任矩阵算好之后M步只需要做三个更新。更新混合系数αₖ时把所有样本对第k个高斯的责任相加再除以总样本数N得到的就是新的混合权重。更新第k个高斯的均值向量μₖ时每个样本用责任r[i,k]作为权重计算加权平均再从权重总和归一化。更新协方差矩阵Σₖ时用带责任的样本构造差异向量外积并做同样的加权平均。这三个步骤本质上就是带权重的高斯极大似然估计只是把传统属于某个类的0/1权重换成了[0,1]之间的连续责任值。循环体写完以后我要在每一轮迭代里打印对数似然肉眼确认数值在单调上升。如果发现似然下降那几乎可以确定是代码里某个归一化或取对数的地方写岔了。3.5 完整示例代码与运行观察下面是一个最小可运行版本只依赖 NumPy 和 SciPyimport numpy as np from scipy.stats import multivariate_normal # 生成三簇模拟数据 np.random.seed(42) true_means [np.array([-2, 0]), np.array([2, 0]), np.array([0, 3])] true_covs [np.array([[0.6, 0.2], [0.2, 0.5]]), np.array([[0.6, -0.2], [-0.2, 0.5]]), np.array([[0.5, 0.1], [0.1, 0.4]])] X [] for m, c in zip(true_means, true_covs): X.append(np.random.multivariate_normal(m, c, 300)) X np.vstack(X) N, D X.shape K 3 # 用K-means初始化参数 from sklearn.cluster import KMeans km KMeans(n_clustersK, random_state0).fit(X) labels km.labels_ alpha np.array([np.mean(labels k) for k in range(K)]) mu np.array([X[labels k].mean(axis0) for k in range(K)]) cov np.array([np.cov(X[labels k], rowvarFalse) 1e-6 * np.eye(D) for k in range(K)]) def log_likelihood(X, alpha, mu, cov): n X.shape[0] log_ll 0.0 for i in range(n): row np.array([ alpha[k] * multivariate_normal.pdf(X[i], meanmu[k], covcov[k]) for k in range(K) ]) log_ll np.log(row.sum()) return log_ll for it in range(100): # E-step: 计算责任矩阵 r[i, k] r np.zeros((N, K)) for k in range(K): r[:, k] alpha[k] * multivariate_normal.pdf(X, meanmu[k], covcov[k]) r r / r.sum(axis1, keepdimsTrue) # M-step: 更新参数 Nk r.sum(axis0) alpha_new Nk / N mu_new (r.T X) / Nk[:, None] cov_new np.zeros((K, D, D)) for k in range(K): diff X - mu_new[k] cov_new[k] (r[:, k][:, None] * diff).T diff / Nk[k] cov_new[k] 1e-6 * np.eye(D) alpha, mu, cov alpha_new, mu_new, cov_new if it % 10 0: print(iter, it, log-likelihood:, log_likelihood(X, alpha, mu, cov))运行的时候你会看到似然函数在前面几次迭代快速爬升随后增量变小最终趋于平缓。由于这个实验的数据是从已知参数生成的收敛得到的中心会和真实中心非常接近。我实测时还发现如果完全换成随机初始化偶尔会让两个高斯分量在数据中互相纠缠最后收敛到局部最优导致某个簇被切分。这个问题也是下一节要聊的重头戏。代码里写逐样本循环是为了让公式对应关系直白数据量大了以后需要改成矩阵化计算不然会很慢。但作为学习版本清晰度优先。4. 实战中必踩的坑与排查手段4.1 局部最优每次结果都不一样EM对初始参数非常敏感尤其当K偏大、数据分布重叠度高时同一份数据跑十次很可能得到几种截然不同的聚类结果。我做项目的标准操作是先用K-means初始化这是降低随机性的第一步。第二步再做多次启动比如10次甚至20次每次固定不同随机种子最后取对数似然最大的那一组参数。第三步如果业务对稳定性有硬性要求还可以先做标准化或PCA降维让数据分布更接近各项同性EM跑起来会稳一些。如果发现不同启动的结果差异极大优先怀疑初始化太差而不是代码写错了。这是EM的固有特性不需要慌。4.2 协方差矩阵奇异大多是某个分量样本太少运行GMM时最常见的报错是LinAlgError: Singular matrix源头往往是某一个高斯分量的责任权重总和接近0或为0导致协方差矩阵无法求逆。当K设得比数据真实分布多或者初始化把某一簇里样本扔得太少时就会出现这种问题。最简单的应对是给每次M步更新后的协方差矩阵加一个很小的正则项例如1e-6乘以单位矩阵。这个操作在代码里已经列出来了虽然增加了极小的偏差但能稳定矩阵求逆。如果加了正则仍然报错就要检查Nk是不是掉到了0附近。若是这样减少K或者用更保守的初始化会更有效。4.3 收敛慢或震荡先检查是不是用了硬标签收敛过程如果出现似然震荡而不是稳定上升通常是在E步做了多余操作。比如有人觉得软标签不如硬标签清晰手动把责任矩阵里最大值置1、其余置0拿硬EM去跑。那样做不是完全不行但它已经不是标准EM更容易出现震荡。还有一种常见问题是没有对责任矩阵做行归一化。如果责任矩阵每一行的和不等于1M步的若干加权平均就会失去概率意义参数更新方向会变得奇怪。排查时打印几行责任矩阵确认每个元素都在0到1之间并且行和严格等于1。4.4 分量标签翻转和可视化陷阱高斯混合模型本身不区分哪个分量排第一、哪个排第二。同一组聚类结果把分量1和分量2互换似然不会变化。这会给代码调试带来困惑两次运行聚类效果一样但中心点看起来对不上甚至标签名字都对不上。遇到这种情形不要用分量序号直接与真实簇类别做比较而是使用调整兰德指数、归一化互信息等不受标签顺序影响的指标。如果必须对比参数就给分量做一个按位置排列的后处理比如按均值向量在空间中的位置排序否则对比结果很容易自我怀疑。4.5 决定K值不只看图还要算账GMM的K值是需要预先指定的。如果数据本身有四类你偏用K3算法会把两类糊成一团。如果你用K5多余的分量又会去抢数据甚至某个高斯形状会被拉成细长条来拟合离群点。一套比较可复现的决策方式是画BIC或AIC曲线看指标下降的拐点。BIC的计算并不麻烦在最优参数时给对数似然加一个与自由参数数量的惩罚项。比起纯粹人眼看散点图BIC更客观也更容易写进流程。不过我要多提醒一句无监督场景里统计上的最优K不一定等于业务想要的分类数。之前做个客服文本聚类项目统计结果偏好K8但业务侧只需要5个标签那就拿K8的结果做归类合并把细节留到后续人工打标而不是机械地追求统计最优。5. 走出课本EM在真实项目里的变体与体会5.1 不止是聚类缺失值填补也能用日常数据处理时我们经常遇到部分字段缺失的情况。设备传感器可能在某段时间断传用户画像里某个字段可能成片缺失。直接删掉缺失行浪费信息用均值填充又太粗糙。如果能将缺失值视作隐变量将已观测字段视作已知部分完全可以用EM框架估计出不带主观先验的填补结果。对于待填补的连续变量E步要计算缺失变量的条件期望M步基于期望后的补全结果重新拟合回归或分类参数。这套逻辑比独立插值更能捕捉变量间的关联结构。不过要记牢普通EM默认缺失机制可以被忽略如果字段是否缺失与字段值本身有关例如高收入用户更不愿意填收入那么直接用EM会有偏差需要引入更复杂的模型才能处理。5.2 从GMM到HMM、LDA方法异曲同工EM骨架可以顺滑地搬到许多经典模型上。隐马尔可夫模型的Baum-Welch算法本质上就是EME步用前向-后向算法估计每个时刻落在某个隐状态的后验概率M步再重估状态转移矩阵和发射概率。LDA主题模型里的训练过程也与EM同源唯一的区别在于它需要用变分分布去近似真正的后验分布不然E步可能不可解。在第九章掌握的那套思路能帮你快速理解其他高级模型的训练步骤。如果真的只是为了用库多数时候我们不需要手写EM库内部已经完成了相关计算。但理解原理仍然能帮我们判断参数设置是否合理以及当模型结果异常时排查问题出在E步还是M步。5.3 手推小例子和调包之间的分寸我自己的习惯是学习任何迭代式算法先在小数据集上手动迭代三到五轮把每轮参数的变化写下来。EM尤其适合这种练法用两个高斯分量、几十个样本把责任矩阵逐轮打印出来你会非常直观地看到软标签怎样从先验均衡走向分化均值向量怎样一步一步移动到聚类中心。等这个流程做完再看网上的高速矩阵版本或者看框架源码里的加速技巧都会有底气。工作后我用现成高斯混合库的机会很多但每当效果异常例如某个类始终不饱满我第一反应仍然是调出中间责任矩阵看几眼看看能量是不是被某个分量吸走了。这种“打开黑盒”的直觉来自当初那几页手推笔记。所以我的建议很朴素第九章别急着跳到sklearn里的GaussianMixture先花半小时手写一个不加加速的最小版本收益比刷十遍教学视频都高。
