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-Maximization），EM 是它在含隐变量概率模型上的一个特例。理解了这一层，你就明白 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_{j=1}^{K} \pi_j \mathcal{N}(x_i|\mu_j,\Sigma_j)}$$

这个式子是贝叶斯公式的直接应用，分母是样本 $x_i$ 的边缘密度，分子是第 $k$ 个成分对它的贡献。$\gamma_{ik}$ 可以理解为"第 $k$ 个成分对解释样本 $i$ 负多大责任"。

2.3 M 步：对 Q 函数求极值

M 步要解的是 $\theta^{(t+1)} = \arg\max_\theta Q(\theta|\theta^{(t)})$。对高斯混合模型，把 Q 函数写开：

$$Q = \sum_{i=1}^{N}\sum_{k=1}^{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^{(t+1)} = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}$$

这就是责任度加权的样本均值。对 $\Sigma_k$ 求导得到加权协方差：

$$\Sigma_k^{(t+1)} = \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^{(t+1)} = \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^{(t+1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t)})$。而真实似然总是不小于下界：$\log p(X|\theta^{(t+1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t+1)})$。三式串起来：

$$\log p(X|\theta^{(t+1)}) \geq \mathcal{L}(q^{(t)},\theta^{(t+1)}) \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_comp=200): 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, replace=False) 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, axis=1)) 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(axis=1, keepdims=True) + 1e-300) return resp

这里有个数值细节值得说：np.sum(diff @ inv * diff, axis=1)算的是马氏距离的平方，利用了 $(x-\mu)^T\Sigma^{-1}(x-\mu)$ 可以写成逐元素乘再求和。比显式写双重循环快一个数量级。归一化时加1e-300是防止某些样本在所有成分下密度都下溢到零导致除零，虽然概率极低，但生产代码里必须防。

M 步直接用前面推的公式，注意用resp.sum(axis=0)得到每个成分的总责任度。

def m_step(X, resp): n, d = X.shape K = resp.shape[1] Nk = resp.sum(axis=0) + 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_iter=200, tol=1e-6, n_init=5): 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_iter=200，对中等规模数据够用。n_init=5是多次重启取最优，能显著降低陷入坏局部最优的概率。

跑一下看看结果：

ll, means, covs, weights, iters = train_em(X, K=3) 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_components=3, covariance_type='full', n_init=5, random_state=42).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_covar=1e-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(d+1)/2$ 个参数，加上 $K-1$ 个权重参数。BIC 对复杂模型惩罚更重，倾向选小 K；AIC 惩罚轻，倾向选大 K。实践中 BIC 更常用，因为它的一致性性质在样本量增大时能选到真实 K。

bics = [] for k in range(1, 8): gm = GaussianMixture(n_components=k, n_init=5, random_state=42).fit(X) bics.append(gm.bic(X)) print(bics)

画出来看拐点，或者直接取最小值。注意 BIC 在成分重叠严重时也会失灵，这时候得结合业务判断——比如你明确知道数据来自三个群体，那就别让 BIC 把你带偏。

4.3 收敛慢与标签切换

EM 在最优附近收敛慢是出了名的，因为它的收敛速率由缺失信息比例决定，缺失信息越多越慢。如果跑了几百轮似然还在缓慢爬升，可以考虑：

• 用更好的初始化（K-means++ 而不是随机）
• 加动量项或者改用 ECM（Expectation Conditional Maximization）变体
• 对参数做缩放，改善条件数

另一个隐蔽的坑是标签切换（label switching）：EM 对成分的编号没有约束，两次运行可能得到完全一样的划分但成分顺序不同。做结果对比时不能直接比均值数组，要先做匹配（比如按均值排序后比）。这个问题在单次训练里看不出来，做多次重启或者交叉验证时才会暴露。

4.4 常见问题速查表

提示：调试 EM 时，把每一轮的对数似然、各成分权重、最小特征值都打印出来，异常往往在数值上早有征兆，等到报错才查就晚了。

5. 从高斯混合往外看：EM 的延展与替代方案

5.1 隐马尔可夫模型里的 Baum-Welch

EM 最著名的应用之一就是隐马尔可夫模型的参数学习，那个算法叫Baum-Welch，本质就是 EM 在序列数据上的特化。E 步用前向-后向算法算状态的后验概率，M 步用这些后验更新转移矩阵和发射矩阵。公式形式和高斯混合高度相似，都是"后验加权 + 重新估计"的套路。理解了高斯混合的 EM，Baum-Welch 基本就是换个概率模型的事。

5.2 变分推断：当后验算不出来时

标准 EM 的 E 步要求后验 $p(Z|X,\theta)$ 能精确计算。有些模型里这个后验不可解，这时候就得用变分推断：不再要求 $q$ 等于后验，而是在一个受限分布族里找最接近后验的 $q$，最大化下界。这催生了变分自编码器那一系的方法。从 EM 到变分推断，核心思想没变——还是构造下界、优化下界，只是下界更松了。

5.3 什么时候该放弃 EM

EM 不是唯一选择。如果模型可微且你愿意调学习率，随机梯度变分贝叶斯（SGVB）那套在深度模型里更主流；如果数据量极大，小批量 EM 或者在线 EM 更合适；如果对全局最优有执念，可以试多次重启加模拟退火。我的经验是：中小规模、模型结构清晰、需要可解释的参数估计时，EM 依然是首选；一旦上了深度网络或者超大数据，就该考虑别的工具了。

6. 一些不那么教科书的心得

写到这里，公式和代码都摆完了，说几个文档里不会写但实际会遇到的点。第一，初始化比算法本身更值得花时间。我见过太多人把 EM 实现得完美无缺，结果因为随机初始化踩了局部最优，回头怀疑公式推错了。多花十分钟设计初始化策略，能省下几小时的调试。第二，对数似然的绝对值不重要，趋势才重要。不同初始化下似然的绝对值可能差很多，但只要单调上升且收敛到相近的划分，就说明结果可信。第三，协方差的正则化不是可选项而是必选项，哪怕你觉得数据很干净，数值误差也会在某个角落埋雷。

还有一个我踩过的坑：早期实现时忘了在 E 步归一化后检查是否有整行为零，结果某个离群点在所有成分下密度都下溢，责任度全零，M 步分母为零直接 NaN。后来养成习惯，每次 E 步后断言resp.sum(axis=1)全为 1，问题立刻暴露。这种防御性检查看着啰嗦，但在调参阶段能救命。

最后分享一个理解 EM 的小视角：把它想成"用当前信念去重新分配信任，再用信任去修正信念"的循环。E 步是分配信任，M 步是修正信念，两者互为因果，迭代到自洽。这个视角比任何公式都更容易记住，也更容易迁移到别的含隐变量问题上。