EM 算法这个名字第一次听到的人多半会愣一下——期望最大化听起来像是某种管理学课程里的名词跟机器学习有什么关系我当初也是这么想的直到真正动手推导了一遍高斯混合模型的参数估计才明白这套东西为什么能在无监督学习里占据这么稳固的位置。这篇内容我打算把 EM 算法从直觉到公式、从推导到代码完整拆一遍面向的是刚接触机器学习不久、想搞明白为什么这么算的读者尤其是那些看了一遍公式仍然觉得云里雾里的朋友。我会尽量用生活化的例子把每一步的动机讲清楚公式该推的地方不省略代码该跑的地方给完整实现读完你应该能自己手写一个能收敛的高斯混合模型训练流程并且知道它什么时候会坑你。1. 先搞清楚 EM 算法到底在解决什么问题1.1 一个数据有隐变量的典型场景假设你面前有一堆身高数据来自两个不同群体的人混在一起但你拿到的记录里只有身高数值没有属于哪个群体这个标签。现在你想估计两个群体各自的身高均值和方差。如果标签已知这就是个再简单不过的统计题——分组求均值方差就完事了。问题恰恰在于标签缺失而缺失的标签又依赖于我们想估计的参数参数估计又依赖于标签形成了一个死循环。这个缺失的标签就是所谓的隐变量latent variable。EM 算法的全部精妙之处就在于它给出了一套打破这个死循环的迭代框架先猜一组参数用这组参数去推断每个样本属于各群体的概率这就是 E 步期望步再拿这些概率作为权重重新估计参数这就是 M 步最大化步如此往复直到收敛。听起来朴素得不像话但它背后有严格的数学保证——每一次迭代都不会让似然函数变小。我用一个更贴近实际的类比来说明。想象你在调配两种不同浓度的糖水但杯子混在一起分不清哪杯是哪种。你先随便假设两个浓度值然后根据这个假设去判断每一杯更可能是哪种浓度软判断给概率而不是硬分类接着按这些概率加权重新算两个浓度的估计值再回头修正判断。反复几轮假设会越来越接近真实浓度。这就是 EM 的骨架没有任何神秘成分。1.2 为什么不能直接用梯度下降有人会问既然目标是最大化似然函数直接对参数求导做梯度上升不行吗技术上可以但实际会遇到两个麻烦。第一对数似然里出现对数里面套求和的结构比如高斯混合模型的似然是 $\log \sum_k \pi_k \mathcal{N}(x|\mu_k,\Sigma_k)$求导之后各项耦合在一起解析解很难直接写出来。第二隐变量的存在让似然曲面变得复杂梯度法对初值敏感、容易停在糟糕的局部区域。EM 的思路完全不同它不去直接啃那个难解的对数似然而是构造一个下界函数Q 函数这个下界在当前位置与真实对数似然相等并且容易优化。优化下界就等于间接抬高了真实似然。这个构造下界—优化下界的套路在数学优化里叫MM 算法Minorize-MaximizationEM 是它在含隐变量概率模型上的一个特例。理解了这一层你就明白 EM 的收敛性不是巧合而是下界性质直接推出来的。1.3 EM 的适用边界在哪里EM 不是什么都能干的万能钥匙。它适用的前提是模型的对数似然在给定隐变量后变得容易优化。高斯混合模型、隐马尔可夫模型、混合专家模型、含缺失数据的回归都满足这个条件。反过来如果你的模型没有隐变量结构或者给定隐变量后依然难解EM 就帮不上忙。另外要清楚它的几个固有特点收敛是单调的似然不会下降这点很让人安心但收敛到的是局部最优不是全局最优初值选得好不好直接决定结果质量收敛速度在接近最优时会变慢因为它在最优点附近的表现接近线性收敛不如牛顿法那种二次收敛来得快。这些特点我在后面的实操部分会结合具体现象再展开。2. 数学推导从对数似然到下界函数2.1 用 Jensen 不等式构造下界设观测数据为 $X$隐变量为 $Z$待估参数为 $\theta$。我们要最大化的对数似然是$$\log p(X|\theta) \log \sum_Z p(X,Z|\theta)$$直接优化右边很难因为对数里面是求和。引入一个关于隐变量的分布 $q(Z)$满足 $\sum_Z q(Z)1$做恒等变形$$\log p(X|\theta) \log \sum_Z q(Z) \frac{p(X,Z|\theta)}{q(Z)}$$这一步只是乘了一个 $q(Z)/q(Z)$没有改变任何东西。接下来用Jensen 不等式对数函数是凹函数所以 $\log \mathbb{E}[Y] \geq \mathbb{E}[\log Y]$。把求和看成对 $q(Z)$ 的期望得到$$\log p(X|\theta) \geq \sum_Z q(Z) \log \frac{p(X,Z|\theta)}{q(Z)}$$右边这个量记作 $\mathcal{L}(q,\theta)$它就是对数似然的下界。Jensen 不等式的取等条件是随机变量为常数也就是 $\frac{p(X,Z|\theta)}{q(Z)}$ 不随 $Z$ 变化这等价于 $q(Z) \propto p(X,Z|\theta)$归一化后就是 $q(Z) p(Z|X,\theta)$即隐变量的后验分布。这个结论极其关键当 $q$ 取成后验分布时下界与真实对数似然相等。所以 EM 的每一轮都是先让下界贴紧真实似然E 步再抬高下界M 步真实似然自然跟着被抬高。这就是单调收敛性的来源。2.2 E 步计算后验并写出 Q 函数E 步做的事就是把 $q(Z)$ 固定为当前参数 $\theta^{(t)}$ 下的后验 $p(Z|X,\theta^{(t)})$然后写出下界作为 $\theta$ 的函数$$Q(\theta|\theta^{(t)}) \sum_Z p(Z|X,\theta^{(t)}) \log p(X,Z|\theta)$$注意这里省略了与 $\theta$ 无关的项 $\sum_Z p(Z|X,\theta^{(t)}) \log p(Z|X,\theta^{(t)})$因为它不影响对 $\theta$ 的优化。所以 Q 函数本质上就是完全数据对数似然在后验分布下的期望。这个名字很直白完全数据指 $(X,Z)$ 都已知时的似然期望是对隐变量的不确定性求的。对于高斯混合模型E 步具体算的是每个样本属于每个成分的责任度responsibility$$\gamma_{ik} \frac{\pi_k \mathcal{N}(x_i|\mu_k,\Sigma_k)}{\sum_{j1}^{K} \pi_j \mathcal{N}(x_i|\mu_j,\Sigma_j)}$$这个式子是贝叶斯公式的直接应用分母是样本 $x_i$ 的边缘密度分子是第 $k$ 个成分对它的贡献。$\gamma_{ik}$ 可以理解为第 $k$ 个成分对解释样本 $i$ 负多大责任。2.3 M 步对 Q 函数求极值M 步要解的是 $\theta^{(t1)} \arg\max_\theta Q(\theta|\theta^{(t)})$。对高斯混合模型把 Q 函数写开$$Q \sum_{i1}^{N}\sum_{k1}^{K} \gamma_{ik} \left[ \log \pi_k - \frac{1}{2}\log|\Sigma_k| - \frac{1}{2}(x_i-\mu_k)^T\Sigma_k^{-1}(x_i-\mu_k) - \frac{d}{2}\log(2\pi) \right]$$对 $\mu_k$ 求导并令其为零得到$$\mu_k^{(t1)} \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}$$这就是责任度加权的样本均值。对 $\Sigma_k$ 求导得到加权协方差$$\Sigma_k^{(t1)} \frac{\sum_i \gamma_{ik}(x_i-\mu_k)(x_i-\mu_k)^T}{\sum_i \gamma_{ik}}$$对 $\pi_k$ 的优化要带约束 $\sum_k \pi_k 1$用拉格朗日乘子法结果是$$\pi_k^{(t1)} \frac{\sum_i \gamma_{ik}}{N}$$即第 $k$ 个成分的平均责任度也等于它软分配到的样本比例。三个更新公式形式高度统一都是用责任度加权的思路这也是 EM 在混合模型上如此优雅的原因。注意协方差更新公式里的分母是 $\sum_i \gamma_{ik}$ 而不是 $N_k - 1$这是最大似然估计的有偏形式。样本量小的时候有偏估计会让协方差偏小实践中如果成分样本很少可以考虑加正则项或者改用有偏/无偏的权衡。2.4 收敛性证明的直觉版本严格证明需要用到下界性质我这里给一个能说服自己的直觉版本。设第 $t$ 轮参数为 $\theta^{(t)}$。E 步让下界在该点等于真实似然$\mathcal{L}(q^{(t)},\theta^{(t)}) \log p(X|\theta^{(t)})$。M 步抬高下界$\mathcal{L}(q^{(t)},\theta^{(t1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t)})$。而真实似然总是不小于下界$\log p(X|\theta^{(t1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t1)})$。三式串起来$$\log p(X|\theta^{(t1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t)}) \log p(X|\theta^{(t)})$$单调性得证。这个证明链条干净利落也解释了为什么 EM 不会出现似然震荡——每一步都在往上走或者持平。3. 高斯混合模型完整实现从公式到可运行代码3.1 数据生成与初始化策略理论讲完动手才有感觉。我用 Python 和 NumPy 从零实现一个二维高斯混合模型的 EM 训练不依赖 sklearn这样每一步都看得见。先造数据三个成分均值分别放在不同位置协方差各不相同每个成分抽若干样本。import numpy as np np.random.seed(42) def generate_data(n_per_comp200): means [np.array([0.0, 0.0]), np.array([5.0, 5.0]), np.array([0.0, 6.0])] covs [np.array([[1.0, 0.3], [0.3, 1.0]]), np.array([[1.5, -0.4], [-0.4, 1.0]]), np.array([[0.8, 0.0], [0.0, 0.8]])] X [] for m, c in zip(means, covs): X.append(np.random.multivariate_normal(m, c, n_per_comp)) return np.vstack(X) X generate_data() print(X.shape) # (600, 2)初始化是 EM 最容易被忽视却最影响结果的环节。常见做法有三种随机选 K 个样本点当初始均值用 K-means 跑一遍拿聚类中心当初始均值多次随机初始化取似然最高的那次。我实测下来K-means 初始化 多次重启最稳因为 K-means 本身对初值有一定鲁棒性而且它给出的硬划分能避免 EM 一开始就陷入退化解。def init_params(X, K): n, d X.shape # 简单版随机选 K 个样本作为初始均值 idx np.random.choice(n, K, replaceFalse) means X[idx].copy() covs np.array([np.eye(d) for _ in range(K)]) weights np.ones(K) / K return means, covs, weights协方差初始化为单位阵是个保守选择避免一开始就给出过于自信的形状假设。权重均匀分配同理。3.2 E 步与 M 步的向量化实现E 步要算责任度矩阵形状是 $N \times K$。逐样本循环在数据量大时会很慢用广播写成矩阵运算。def gaussian_pdf(X, mean, cov): d X.shape[1] diff X - mean inv np.linalg.inv(cov) det np.linalg.det(cov) norm 1.0 / np.sqrt((2 * np.pi) ** d * det) exp_term np.exp(-0.5 * np.sum(diff inv * diff, axis1)) return norm * exp_term def e_step(X, means, covs, weights): n, K X.shape[0], len(means) resp np.zeros((n, K)) for k in range(K): resp[:, k] weights[k] * gaussian_pdf(X, means[k], covs[k]) # 归一化防止除零加个极小值 resp resp / (resp.sum(axis1, keepdimsTrue) 1e-300) return resp这里有个数值细节值得说np.sum(diff inv * diff, axis1)算的是马氏距离的平方利用了 $(x-\mu)^T\Sigma^{-1}(x-\mu)$ 可以写成逐元素乘再求和。比显式写双重循环快一个数量级。归一化时加1e-300是防止某些样本在所有成分下密度都下溢到零导致除零虽然概率极低但生产代码里必须防。M 步直接用前面推的公式注意用resp.sum(axis0)得到每个成分的总责任度。def m_step(X, resp): n, d X.shape K resp.shape[1] Nk resp.sum(axis0) 1e-10 means (resp.T X) / Nk[:, None] covs np.zeros((K, d, d)) for k in range(K): diff X - means[k] covs[k] (resp[:, k][:, None] * diff).T diff / Nk[k] weights Nk / n return means, covs, weights(resp.T X) / Nk[:, None]这一行就是加权均值的矩阵形式resp.T是 $K \times N$乘 $X$ 是 $N \times d$得到 $K \times d$再按成分样本数归一化。协方差那行用resp[:, k][:, None] * diff给每个样本的偏差加权再转置相乘得到 $d \times d$ 的加权协方差。3.3 对数似然监控与收敛判定训练循环里要监控对数似然一方面判断收敛一方面验证单调性。def log_likelihood(X, means, covs, weights): n X.shape[0] K len(means) ll np.zeros(n) for k in range(K): ll weights[k] * gaussian_pdf(X, means[k], covs[k]) return np.sum(np.log(ll 1e-300)) def train_em(X, K, max_iter200, tol1e-6, n_init5): best None for _ in range(n_init): means, covs, weights init_params(X, K) prev_ll -np.inf for it in range(max_iter): resp e_step(X, means, covs, weights) means, covs, weights m_step(X, resp) ll log_likelihood(X, means, covs, weights) if ll - prev_ll tol: break prev_ll ll if best is None or ll best[0]: best (ll, means, covs, weights, it) return best收敛判据用相邻两轮对数似然的增量小于tol。这里有个坑tol设太小会跑满max_iter设太大会提前停。我一般用1e-6配合max_iter200对中等规模数据够用。n_init5是多次重启取最优能显著降低陷入坏局部最优的概率。跑一下看看结果ll, means, covs, weights, iters train_em(X, K3) print(log-likelihood:, ll) print(means:\n, means) print(weights:, weights) print(iterations:, iters)实测下来三个成分的均值会收敛到接近真实值的位置权重接近 1/3。如果某次初始化不好可能出现两个成分挤在一起、第三个成分覆盖一大片的情况这就是局部最优的典型表现多次重启能缓解但不能根除。3.4 用 sklearn 做交叉验证自己写完再拿成熟库对一遍是确认实现正确的好习惯。sklearn 的GaussianMixture接口很简洁from sklearn.mixture import GaussianMixture gm GaussianMixture(n_components3, covariance_typefull, n_init5, random_state42).fit(X) print(gm.means_) print(gm.weights_) print(gm.lower_bound_) # 对数似然下界对比自己实现的结果均值和权重应该在同一量级。如果差很多多半是初始化或者协方差更新写错了。lower_bound_是 sklearn 内部记录的对数似然可以拿来和自己算的对比误差在数值精度范围内就说明实现没问题。4. 实操中踩过的坑与排查手册4.1 协方差矩阵奇异最常见的崩溃点EM 训练高斯混合模型时最经典的失败模式是某个成分坍缩到一个点上协方差矩阵趋于奇异行列式趋近零密度函数爆炸对数似然冲到正无穷。现象是训练日志里似然突然跳到一个巨大的值然后程序报LinAlgError或者结果完全失去意义。根因是最大似然估计在高斯混合里本身是无界的——只要让某个成分的方差趋于零、均值对准某个样本点似然就能无限大。这不是数值 bug是模型设定问题。解决办法有几种我按推荐程度排序给协方差加正则在协方差对角线上加一个小的正数 $\epsilon I$sklearn 里对应reg_covar参数默认1e-6。这是最省事的做法代价是轻微偏差。限制协方差类型改用diagonal或spherical参数少、更稳但表达能力下降。贝叶斯化的先验给协方差加逆 Wishart 先验从根上避免奇异但实现复杂。早停监控似然出现异常跳变就停回退到上一轮参数。我一般默认开reg_covar1e-6数据维度高或者样本少的时候调到1e-4。4.2 成分数 K 怎么选K 是超参数EM 本身不会告诉你该用几个成分。常用手段是跑不同 K 值比较BIC贝叶斯信息准则或AIC$$\text{BIC} -2\log p(X|\hat\theta) p\log N$$其中 $p$ 是自由参数个数。高斯混合模型 full 协方差下每个成分有 $d d(d1)/2$ 个参数加上 $K-1$ 个权重参数。BIC 对复杂模型惩罚更重倾向选小 KAIC 惩罚轻倾向选大 K。实践中 BIC 更常用因为它的一致性性质在样本量增大时能选到真实 K。bics [] for k in range(1, 8): gm GaussianMixture(n_componentsk, n_init5, random_state42).fit(X) bics.append(gm.bic(X)) print(bics)画出来看拐点或者直接取最小值。注意 BIC 在成分重叠严重时也会失灵这时候得结合业务判断——比如你明确知道数据来自三个群体那就别让 BIC 把你带偏。4.3 收敛慢与标签切换EM 在最优附近收敛慢是出了名的因为它的收敛速率由缺失信息比例决定缺失信息越多越慢。如果跑了几百轮似然还在缓慢爬升可以考虑用更好的初始化K-means 而不是随机加动量项或者改用 ECMExpectation Conditional Maximization变体对参数做缩放改善条件数另一个隐蔽的坑是标签切换label switchingEM 对成分的编号没有约束两次运行可能得到完全一样的划分但成分顺序不同。做结果对比时不能直接比均值数组要先做匹配比如按均值排序后比。这个问题在单次训练里看不出来做多次重启或者交叉验证时才会暴露。4.4 常见问题速查表现象可能原因排查与解决似然突然跳到极大值协方差坍缩加reg_covar检查是否有成分只分到极少样本似然震荡不单调实现 bug 或数值误差检查 E 步归一化、M 步分母确认下界性质收敛极慢初值差或缺失信息多换 K-means 初始化增加重启次数结果每次不同随机初始化固定随机种子或增加n_init取最优某成分权重趋近零该成分冗余减少 K或加权重先验高维下内存爆炸协方差矩阵 $K \times d \times d$改用对角协方差或分块计算提示调试 EM 时把每一轮的对数似然、各成分权重、最小特征值都打印出来异常往往在数值上早有征兆等到报错才查就晚了。5. 从高斯混合往外看EM 的延展与替代方案5.1 隐马尔可夫模型里的 Baum-WelchEM 最著名的应用之一就是隐马尔可夫模型的参数学习那个算法叫Baum-Welch本质就是 EM 在序列数据上的特化。E 步用前向-后向算法算状态的后验概率M 步用这些后验更新转移矩阵和发射矩阵。公式形式和高斯混合高度相似都是后验加权 重新估计的套路。理解了高斯混合的 EMBaum-Welch 基本就是换个概率模型的事。5.2 变分推断当后验算不出来时标准 EM 的 E 步要求后验 $p(Z|X,\theta)$ 能精确计算。有些模型里这个后验不可解这时候就得用变分推断不再要求 $q$ 等于后验而是在一个受限分布族里找最接近后验的 $q$最大化下界。这催生了变分自编码器那一系的方法。从 EM 到变分推断核心思想没变——还是构造下界、优化下界只是下界更松了。5.3 什么时候该放弃 EMEM 不是唯一选择。如果模型可微且你愿意调学习率随机梯度变分贝叶斯SGVB那套在深度模型里更主流如果数据量极大小批量 EM 或者在线 EM 更合适如果对全局最优有执念可以试多次重启加模拟退火。我的经验是中小规模、模型结构清晰、需要可解释的参数估计时EM 依然是首选一旦上了深度网络或者超大数据就该考虑别的工具了。6. 一些不那么教科书的心得写到这里公式和代码都摆完了说几个文档里不会写但实际会遇到的点。第一初始化比算法本身更值得花时间。我见过太多人把 EM 实现得完美无缺结果因为随机初始化踩了局部最优回头怀疑公式推错了。多花十分钟设计初始化策略能省下几小时的调试。第二对数似然的绝对值不重要趋势才重要。不同初始化下似然的绝对值可能差很多但只要单调上升且收敛到相近的划分就说明结果可信。第三协方差的正则化不是可选项而是必选项哪怕你觉得数据很干净数值误差也会在某个角落埋雷。还有一个我踩过的坑早期实现时忘了在 E 步归一化后检查是否有整行为零结果某个离群点在所有成分下密度都下溢责任度全零M 步分母为零直接 NaN。后来养成习惯每次 E 步后断言resp.sum(axis1)全为 1问题立刻暴露。这种防御性检查看着啰嗦但在调参阶段能救命。最后分享一个理解 EM 的小视角把它想成用当前信念去重新分配信任再用信任去修正信念的循环。E 步是分配信任M 步是修正信念两者互为因果迭代到自洽。这个视角比任何公式都更容易记住也更容易迁移到别的含隐变量问题上。