""" ============================================================================= 第 5 讲 · 代码 2：高斯混合模型与 EM —— ELBO、单调性、k-means 是"硬"EM、似然的塌缩 ============================================================================= 对应讲义：lectures/lecture05_unsupervised.md 第 3 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- EM 的公式你大概背得出来。这里用数值实验把"为什么它能工作、什么时候会出问题"逐条验证： 演示 1：EM 拟合 3 个成分的 GMM。对数似然逐轮单调不减（assert 检查）；学到的参数和真值对比； 多次重启的结果。 演示 2：EM 为什么单调？ELBO 的链条 ℓ(θ_t) = ELBO(Q_t, θ_t) ≤ ELBO(Q_t, θ_{t+1}) ≤ ℓ(θ_{t+1})，逐轮数值验证； 对任意的 Q，ℓ − ELBO 恰好等于 Σᵢ KL(Q_i ‖ p(z | x⁽ⁱ⁾))。 演示 3：soft vs hard。(a) 固定 Σ = σ²I 的"软 k-means"在 σ → 0 时收敛到 k-means 的解（k-means = hard EM）； (b) 在代码 01 的两份"k-means 失败"的数据上，完整的 GMM 能做到多好。 演示 4：似然无界：把一个成分放在一个数据点上、让它的方差 → 0，ℓ → +∞。多次随机重启中真的会发生； 协方差正则化（MAP）阻止它。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 02_gmm_em.py # 约 4 秒 依赖本目录的 unsupervised.py。 """ import itertools import numpy as np from unsupervised import (best_perm_accuracy, elbo, fit_gmm_em, fmt, gmm_e_step, gmm_init, gmm_log_joint, gmm_m_step, init_kmeanspp, kmeans, kmeans_best_of, make_spec, print_header, soft_kmeans, spec_params) def gmm3_data(): """代码 02 演示 1~3 共用的数据：CLUSTER_SPECS["gmm3"]，600 个点，固定种子。""" return make_spec("gmm3", np.random.default_rng(21)) def align(mu_true, mu_est): """ 找一个排列 perm，使 mu_est[perm[k]] 和 mu_true[k] 最接近（成分的编号是任意的 —— label switching）。 【返回】perm：长度 K 的元组。K 很小（≤ 6），直接枚举所有排列。 """ K = len(mu_true) return min(itertools.permutations(range(K)), key=lambda p: sum(np.sum((mu_true[k] - mu_est[p[k]]) ** 2) for k in range(K))) # ============================================================================= # 演示 1：EM 拟合 GMM # ============================================================================= def demo_em_fit(): print_header("演示 1：EM for a 3-component GMM —— log-likelihood never decreases") X, z = gmm3_data() truth = spec_params("gmm3") hist = [] params, Q, ll, n_it, status = fit_gmm_em(X, gmm_init(X, 3, np.random.default_rng(22)), history=hist) monotone = all(b >= a - 1e-12 * abs(a) for a, b in zip(hist, hist[1:])) assert monotone, "EM decreased the log-likelihood" ll_truth = gmm_e_step(X, truth)[1] print(f"数据：n = {len(X)}，3 个成分（真实权重 {fmt(truth['phi'], '.2f')}）；随机初始化，{status}，{n_it} 轮\n") show = [0, 1, 2, 3, 5, 10, 15, 20, n_it] print(f"{'iter':>6}{'log-likelihood':>17}{'increase':>12}") for t in show: inc = "" if t == 0 else f"{hist[t] - hist[t - 1]:.3g}" print(f"{t:>6}{hist[t]:>17.3f}{inc:>12}") print(f"\n 对数似然单调不减：{monotone}（assert 检查了全部 {len(hist) - 1} 步）") print(f" ℓ(EM 的解) = {ll:.3f}；ℓ(生成数据的真参数) = {ll_truth:.3f}") perm = align(truth["mu"], params["mu"]) print(f"\n{'comp':>5}{'phi true':>10}{'phi EM':>9} {'mu true':<16}{'mu EM':<18}{'Sigma true [11,12,22]':<24}" f"{'Sigma EM [11,12,22]'}") for k in range(3): j = perm[k] St, Se = truth["Sigma"][k], params["Sigma"][j] print(f"{k:>5}{truth['phi'][k]:>10.3f}{params['phi'][j]:>9.3f} {fmt(truth['mu'][k], '.2f'):<16}" f"{fmt(params['mu'][j], '.2f'):<18}{fmt([St[0, 0], St[0, 1], St[1, 1]], '.2f'):<24}" f"{fmt([Se[0, 0], Se[0, 1], Se[1, 1]], '.2f')}") acc = best_perm_accuracy(z, Q.argmax(axis=1), 3) acc_bayes = best_perm_accuracy(z, gmm_e_step(X, truth)[0].argmax(axis=1), 3) print("\n 另外 5 次随机重启：") print(f" {'restart':>8}{'iters':>7}{'final log-lik':>15}") finals = [] for r in range(5): _, _, llr, itr, st = fit_gmm_em(X, gmm_init(X, 3, np.random.default_rng(230 + r))) finals.append(llr) print(f" {r:>8}{itr:>7}{llr:>15.3f} {st}") same = max(finals) - min(finals) < 1e-3 print(f""" 观察 Observations： · 每一轮对数似然都上升或不变（演示 2 解释为什么）。第 1 轮 +{hist[1] - hist[0]:.0f}；之后有一段很长的"爬坡" （第 10 轮、第 15 轮每轮仍有 +{hist[10] - hist[9]:.0f}、+{hist[15] - hist[14]:.0f}，成分在慢慢调整形状和归属）， 最后一轮只有 +{hist[n_it] - hist[n_it - 1]:.1e}。EM 的收敛是线性的（一阶方法），成分重叠越多越慢。 · 学到的参数接近真值（成分的编号是任意的，表里已经按均值对齐）。用后验的 argmax 当聚类标签， 准确率 {acc:.3f}；用真参数算后验（Bayes 最优）是 {acc_bayes:.3f} —— 成分有重叠，再好的方法也到不了 1。 · ℓ(EM) = {ll:.1f} 比 ℓ(真参数) = {ll_truth:.1f} 还高：MLE 拟合的是这一份样本，自然比生成它的参数更"贴合"它。 这是过拟合的最温和的形式，n 越大差距越小。 · 5 次重启{'都到达了同一个对数似然（差别 < 0.001）' if same else '到达了不同的对数似然 —— 局部最优'}。 EM 和 k-means 一样只保证局部最优；这份数据的成分分得比较开，所以这里没有遇到。""") # ============================================================================= # 演示 2：ELBO 的链条 # ============================================================================= def demo_elbo(): print_header("演示 2：why EM is monotone —— ℓ(θ_t) = ELBO(Q_t,θ_t) ≤ ELBO(Q_t,θ_{t+1}) ≤ ℓ(θ_{t+1})") X, _ = gmm3_data() params = gmm_init(X, 3, np.random.default_rng(22)) # 和演示 1 相同的初始化 print("Q_t = E 步得到的后验（用 θ_t 算）；θ_{t+1} = M 步的结果。gap = ℓ(θ_{t+1}) − ELBO(Q_t, θ_{t+1})" " = Σᵢ KL(Q_t,i ‖ p(z|x⁽ⁱ⁾; θ_{t+1})) ≥ 0\n") print(f"{'t':>3}{'l(theta_t)':>14}{'ELBO(Q_t,th_t)-l':>18}{'ELBO(Q_t,th_t+1)':>18}{'l(theta_t+1)':>14}{'gap (KL)':>11}") worst_tight, worst_order = 0.0, 0.0 for t in range(200): Q, ll, LJ = gmm_e_step(X, params) # E 步：Q_t，同时得到 ℓ(θ_t) e_tight = elbo(LJ, Q) # ELBO(Q_t, θ_t)：应该等于 ℓ(θ_t) new = gmm_m_step(X, Q) # M 步：θ_{t+1} e_new = elbo(gmm_log_joint(X, new), Q) # ELBO(Q_t, θ_{t+1})：Q 不变，θ 换成新的 ll_new = gmm_e_step(X, new)[1] # ℓ(θ_{t+1}) worst_tight = max(worst_tight, abs(e_tight - ll)) worst_order = min(worst_order, e_new - e_tight, ll_new - e_new) if t < 6: print(f"{t:>3}{ll:>14.3f}{e_tight - ll:>18.1e}{e_new:>18.3f}{ll_new:>14.3f}{ll_new - e_new:>11.3f}") params = new assert worst_tight < 1e-6 and worst_order > -1e-6 print(f"\n 前 200 轮：max |ELBO(Q_t,θ_t) − ℓ(θ_t)| = {worst_tight:.1e}；" f"两个 ≤ 的最小余量 = {worst_order:.1e}（在舍入误差 ~1e-12 以内 ≥ 0：链条成立）") # 任意的 Q：均匀分布、k-means 的硬分配。ℓ − ELBO 应该恰好等于 Σᵢ KL(Q_i ‖ 后验)。 post, ll, LJ = gmm_e_step(X, params) _, km_labels, _, _ = kmeans(X, params["mu"]) Qs = {"posterior (E-step)": post, "uniform 1/K": np.full_like(post, 1.0 / 3), "hard (one-hot)": np.eye(3)[km_labels]} print(f"\n 在收敛后的 θ 处，换不同的 Q（ℓ(θ) = {ll:.3f}）：\n") print(f" {'Q':<20}{'ELBO(Q,theta)':>15}{'l - ELBO':>12}{'sum_i KL(Q_i||post_i)':>23}") rows = {} for name, Qx in Qs.items(): e = elbo(LJ, Qx) logQ = np.log(np.where(Qx > 0, Qx, 1.0)) kl = float(np.sum(Qx * (logQ - np.log(post)))) # KL(Q‖P) = Σ_k Q log(Q/P) rows[name] = (e, ll - e, kl) print(f" {name:<20}{e:>15.3f}{ll - e:>12.3f}{kl:>23.3f}") print(f""" 观察 Observations： · 第 3 列是 0 或 ±1e-13（舍入误差）：E 步之后 ELBO 和对数似然**相等** —— Jensen 不等式的等号成立（Q 取后验时， log 里面的 p(x, z)/Q(z) = p(x) 是常数，"期望的 log" = "log 的期望"）。 · M 步把 ELBO 推高（第 4 列 > 第 2 列对应的 ℓ），之后重新做 E 步，ℓ 又至少涨了 gap = KL 那么多。 所以 ℓ 单调不减：EM 是在 ELBO(Q, θ) 上做**坐标上升**，和 k-means 在 J 上做坐标下降是同一个结构。 · 换别的 Q 时，ℓ − ELBO 恰好等于 Σᵢ KL(Q_i ‖ p(z | x⁽ⁱ⁾))（两列逐位相同）。均匀的 Q 差了 {rows['uniform 1/K'][1]:.0f}； k-means 的硬分配只差 {rows['hard (one-hot)'][1]:.1f} —— 成分分得开时，后验本来就接近 one-hot。 · 这个分解 ℓ(θ) = ELBO(Q, θ) + KL(Q ‖ p(z|x; θ)) 是变分推断（variational inference）和 VAE 的出发点： 后验算不出来时，就在一个容易算的分布族里找 Q，让 KL 尽量小（讲义 7.3 节）。""") # ============================================================================= # 演示 3：soft vs hard # ============================================================================= def demo_soft_vs_hard(): print_header("演示 3：soft vs hard —— k-means is the σ → 0 limit of EM; full GMMs fix k-means' failures") X, _ = gmm3_data() C0 = init_kmeanspp(X, 3, np.random.default_rng(31)) Ck, lk, _, n_km = kmeans(X, C0) print("(a) \"软 k-means\"：GMM 的权重固定为 1/K、协方差固定为 σ²I，只用 EM 更新均值。和 k-means 从同一组初始中心出发。\n") print(f"{'sigma':>7}{'iters':>7}{'mean max resp':>15}{'same label as k-means':>23}{'max |mu - mu_kmeans|':>22}") rows = [] for s in [3.0, 2.0, 1.0, 0.5, 0.25, 0.1, 0.05, 0.02, 0.01]: C, R, it = soft_kmeans(X, C0, s) maxr = float(R.max(axis=1).mean()) agree = float(np.mean(R.argmax(axis=1) == lk)) diff = float(np.max(np.abs(C - Ck))) rows.append((s, maxr, agree, diff)) print(f"{s:>7.2f}{it:>7}{maxr:>15.4f}{agree:>23.4f}{diff:>22.2e}") print(f" （k-means 从同一组初始中心出发，{n_km} 轮收敛）") print("\n(b) 代码 01 演示 4 的两份数据（同样的种子）：k-means vs 完整的 GMM（各重启 10 次取最好）vs Bayes 最优\n") print(f"{'data':<12}{'k-means acc':>13}{'GMM acc':>10}{'Bayes acc':>11}{'GMM log-lik':>13}{'truth log-lik':>15}") res = {} for name, seed in [("elongated", 41), ("unequal", 42)]: Xs, ys = make_spec(name, np.random.default_rng(seed)) _, lab_km, _ = kmeans_best_of(Xs, 2, np.random.default_rng(seed + 100), n_init=10) best = None for r in range(10): _, Qg, llg, _, st = fit_gmm_em(Xs, gmm_init(Xs, 2, np.random.default_rng(seed + 200 + r))) if st in ("converged", "max_iters") and (best is None or llg > best[1]): best = (Qg, llg) truth = spec_params(name) Qt, ll_t, _ = gmm_e_step(Xs, truth) a_km = best_perm_accuracy(ys, lab_km, 2) a_g = best_perm_accuracy(ys, best[0].argmax(axis=1), 2) a_b = best_perm_accuracy(ys, Qt.argmax(axis=1), 2) res[name] = (a_km, a_g, a_b) print(f"{name:<12}{a_km:>13.3f}{a_g:>10.3f}{a_b:>11.3f}{best[1]:>13.1f}{ll_t:>15.1f}") big = rows[0] exact_from = max(r[0] for r in rows if r[3] < 1e-8) exact_diff = max(r[3] for r in rows if r[0] <= exact_from) print(f""" 观察 Observations： · (a) σ 越小，责任矩阵越接近 one-hot（平均最大责任 → 1）。σ ≤ {exact_from:g} 时软 k-means 和 k-means 是**同一个解**（中心的差别 ≤ {exact_diff:.0e}，σ = {rows[-1][0]:g} 时恰好是 {rows[-1][3]:g}）。这就是"k-means = hard EM"的精确含义：E 步的软分配在 σ → 0 时 变成"分给最近的中心"，M 步的加权平均变成"簇内平均"。 · 中间的 σ（1 到 0.05）给出的解和 k-means 差 1e-2 量级、有几个点归属不同 —— 而且不随 σ 单调变小。 原因是两者走了不同的路径：迭代早期中心还没到位时，边界附近的点被"软"地分配，结果停在了 k-means 的另一个 相邻的不动点上（σ = 0.05 时最终的责任已经全是 one-hot，但差别仍然存在）。k-means 本身就有很多不动点。 · σ 很大（{big[0]:g}）时发生了另一件事：平均最大责任 = {big[1]:.4f} ≈ 1/3，三个中心合并成了一个（都在数据的均值处）。 σ 是一个"分辨率"：比簇间距离还大时，模型看不出有几个簇。（deterministic annealing 就是从大 σ 慢慢降下来。） · (b) 同样两份数据，k-means 的准确率是 {res['elongated'][0]:.3f} / {res['unequal'][0]:.3f}；GMM 是 {res['elongated'][1]:.3f} / {res['unequal'][1]:.3f}， Bayes 最优是 {res['elongated'][2]:.3f} / {res['unequal'][2]:.3f}。每个成分有自己的协方差（能表示"细长"）和权重 （能表示"大小悬殊"），代码 01 里 k-means 的两个失败就都消失了。 · 代价：GMM 的参数更多（每个成分 d + d(d+1)/2 + 1 个），需要更多数据；协方差可能塌缩（演示 4）。""") # ============================================================================= # 演示 4：似然无界与塌缩 # ============================================================================= def demo_collapse(): print_header("演示 4：the GMM likelihood is unbounded —— collapse onto a data point, and how a prior stops it") X, _ = gmm3_data() base, _, ll_base, _, _ = fit_gmm_em(X, gmm_init(X, 3, np.random.default_rng(22))) lam = 0.01 print(f"(a) 在演示 1 学到的 3 成分 GMM 上再加第 4 个成分：均值 = 第 0 个数据点 {fmt(X[0], '.2f')}，" f"协方差 σ²I，权重 0.01（其余权重 × 0.99）。") print(f" 原来的 ℓ = {ll_base:.2f}。penalized = ℓ − (λ/2) Σ_k tr(Σ_k⁻¹)，λ = {lam}\n") print(f"{'sigma':>9}{'log-likelihood':>17}{'penalized (MAP)':>18}") lls = [] for s in [1.0, 1e-1, 1e-2, 1e-4, 1e-6, 1e-8]: p = {"phi": np.append(base["phi"] * 0.99, 0.01), "mu": np.vstack([base["mu"], X[0]]), "Sigma": np.concatenate([base["Sigma"], [s ** 2 * np.eye(2)]])} ll = gmm_e_step(X, p)[1] pen = sum(float(np.trace(np.linalg.inv(S))) for S in p["Sigma"]) lls.append(ll) print(f"{s:>9.0e}{ll:>17.2f}{ll - 0.5 * lam * pen:>18.6g}") per_decade = (lls[-1] - lls[-2]) / 2 print("\n(b) 真的会发生吗？一团高斯噪声（没有任何簇结构），成分数偏多，从 100 个随机初始化跑 EM：\n") print(f"{'setting':<16}{'reg':>6}{'converged':>11}{'collapsed':>11}{'empty':>7}{'max_iters':>11}") stats, example = {}, None for (n, K) in [(30, 4), (40, 6)]: Xn = np.random.default_rng(1000 + n).normal(size=(n, 2)) for reg in [0.0, lam]: counts = {"converged": 0, "collapsed": 0, "empty": 0, "max_iters": 0} for r in range(100): hist = [] p_end, _, _, it, st = fit_gmm_em(Xn, gmm_init(Xn, K, np.random.default_rng(5000 + r)), reg=reg, max_iters=1000, history=hist) counts[st] += 1 if st == "collapsed" and example is None: example = (n, K, r, it, hist, p_end) # 记下第一个塌缩的例子，后面详细看 stats[(n, K, reg)] = counts print(f"{'n=' + str(n) + ', K=' + str(K):<16}{reg:>6g}{counts['converged']:>11}{counts['collapsed']:>11}" f"{counts['empty']:>7}{counts['max_iters']:>11}") c0, c1 = stats[(30, 4, 0.0)], stats[(40, 6, 0.0)] n_ex, K_ex, r_ex, it_ex, hist_ex, p_ex = example eigs = [float(np.linalg.eigvalsh(S)[0]) for S in p_ex["Sigma"]] k_bad = int(np.argmin(eigs)) print(f"\n 第一个塌缩的例子（n={n_ex}, K={K_ex}, 第 {r_ex} 次重启）：在第 {it_ex} 轮被检测到。最后 6 轮的对数似然：") print(" " + " ".join(f"{v:.2f}" for v in hist_ex[-6:])) print(f" 塌缩的成分 {k_bad}：N_k = Σᵢ Q_ik = {p_ex['phi'][k_bad] * n_ex:.3f} 个点，" f"协方差的特征值 = {fmt(np.linalg.eigvalsh(p_ex['Sigma'][k_bad]), '.1e')}") empties = sum(stats[(n, K, lam)]["empty"] for (n, K) in [(30, 4), (40, 6)]) collapses_reg = sum(stats[(n, K, lam)]["collapsed"] for (n, K) in [(30, 4), (40, 6)]) if empties > 0: empty_note = (f"代价是出现了 {empties} 次另一种退化 “empty”：\n 某个成分的总责任 → 0（模型“放弃”了一个多余的成分）。" "惩罚项给小成分的下限是 λ/N_k，N_k 越小它越宽、越抢不到点。") else: empty_note = ("这里也没有出现另一种退化 “empty”（某个成分的总责任 → 0）。\n 惩罚项给小成分的下限是 λ/N_k，" "N_k 越小它越宽、越抢不到点，所以多余的成分有可能慢慢“饿死”（编写时在另一份数据上见过）。") print(f""" 观察 Observations： · (a) σ 每缩小 10 倍，ℓ 大约增加 {per_decade:.2f}（= d·ln 10 = 2 × 2.30）：那个数据点的密度 ∝ 1/σ^d → ∞， 而其他点几乎不受影响（它们由另外 3 个成分解释）。**ℓ 没有上界**，"最大似然估计"在字面意义上不存在； 我们平时说的"GMM 的 MLE"其实是指一个"好的"局部最大值（每个成分都覆盖一群点），而不是全局的 sup。 · σ 从 1 降到 0.1 时 ℓ 先**下降**：多出来的成分拿走了 1% 的权重，却只解释了一小块区域；σ 再小下去， 那一个点的密度 1/σ^d 才压倒一切。 · 惩罚项 −(λ/2) tr(Σ⁻¹) 在 σ → 0 时以 1/σ² 的速度趋于 −∞，远快于 log(1/σ) 的增长，于是 MAP 目标函数有界。 · (b) 没有正则化时，{c0['collapsed']}/100（n=30, K=4）和 {c1['collapsed']}/100（n=40, K=6）的随机初始化最终让某个成分塌缩。 例子里塌缩的成分只"拥有" {p_ex['phi'][k_bad] * n_ex:.1f} 个点，协方差有一个特征值掉到了 ~1e-10 以下： 在 d 维空间里，一个成分只要抓住 ≤ d 个点，就能塌缩到它们张成的点 / 直线上，不一定非得是 1 个点。 检测到之后代码就停下（再算下去，Mahalanobis 距离会数值溢出）。 · 加上 λ = {lam} 后塌缩次数是 {collapses_reg}。{empty_note} sklearn 的 reg_covar 是更简单的常数下限 εI。 · 实践中的对策：协方差加下限 / 先验（MAP）、检测塌缩后重新初始化该成分、多次重启、选更少的成分（BIC）。 更一般的教训：**MLE 本身可能是病态的**。第 1 讲的可分数据（‖θ‖ → ∞）是同一类问题，修法也一样：加先验。""") def main(): demo_em_fit() demo_elbo() demo_soft_vs_hard() demo_collapse() if __name__ == "__main__": main()