""" ============================================================================= 第 5 讲 · 代码 5：ICA —— 鸡尾酒会问题；为什么白化不够、为什么高斯源不行、源分布的假设有多重要 ============================================================================= 对应讲义：lectures/lecture05_unsupervised.md 第 6 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- 演示 1：鸡尾酒会。3 个源（正弦、锯齿波、Laplace 噪声）被随机矩阵 A 混合。混合信号和每个源都相关； PCA 白化之后协方差是 I（不相关了），但每个分量**仍然是混合的**。 演示 2：白化之后只剩一个旋转要找。在 2 个源的情形，把 logistic 先验下的似然画成旋转角的函数： 高斯源 —— 完全平坦（不可识别）；Laplace 源 —— 在正确的角度最大； 正弦 + 锯齿波 —— 在 45°（最混合的角度）最大！ 演示 3：CS229 的 sigmoid ICA（批量梯度上升）。(a) 3 个 super-Gaussian 源：完美分离（只差排列和缩放）； (b) 正弦 / 锯齿波 / Laplace：只分出了 Laplace，另外两个停在 45° 混合； (c) 换成按峰度切换的源分布（extended infomax）：全部分离。 演示 4：1 个 Laplace + 2 个高斯源。Laplace 每次都被找到；两个高斯源的"答案"随初始化、迭代步数和 抽样噪声漂移，而似然几乎不变 —— 高斯子空间里的旋转无法确定。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 05_ica.py # 约 10 秒 依赖本目录的 unsupervised.py。 """ import numpy as np from unsupervised import abs_corr, fit_ica, fmt, ica_logistic_loglik, make_sources, print_header, whiten N = 5000 SCALE = np.pi / np.sqrt(3.0) # logistic 分布（scale 1）的标准差是 π/√3：把单位方差的信号放大到同样的尺度 def mixed(kinds, seed): """ 生成源 S、随机混合矩阵 A、观测 X = S Aᵀ（每行一个时刻，x = A s），以及白化后的 Z。 【返回】(S, A, X, Z, W_white) 【边界】A 是标准正态的随机矩阵，几乎一定可逆；只要可逆，ICA 就能在原理上恢复源（不需要 A 正交）。 """ rng = np.random.default_rng(seed) S = make_sources(N, kinds, rng) A = rng.normal(size=(len(kinds), len(kinds))) X = S @ A.T Z, W_white, _ = whiten(X) return S, A, X, Z, W_white def print_corr(title, C, row_name, col_names): """打印一个 |相关系数| 矩阵：行是恢复出来的分量，列是真实的源。""" print(f" {title}") print(f" {'':<6}" + "".join(f"{c:>10}" for c in col_names)) for i, row in enumerate(C): print(f" {row_name + str(i + 1):<6}" + "".join(f"{v:>10.3f}" for v in row)) # ============================================================================= # 演示 1：鸡尾酒会；白化只去相关，不分离 # ============================================================================= def demo_cocktail_party(): print_header("演示 1：the cocktail party —— whitening decorrelates, but does not unmix") kinds = ["sine", "sawtooth", "laplace"] S, A, X, Z, _ = mixed(kinds, 51) kurt = np.mean(S ** 4, axis=0) print(f"3 个源，n = {N}。峰度 E[s⁴]（单位方差下；高斯 = 3）：" + "，".join(f"{k} {v:.2f}" for k, v in zip(kinds, kurt))) print(f"混合矩阵 A = {fmt(A[0], '.2f')} / {fmt(A[1], '.2f')} / {fmt(A[2], '.2f')}（按行），cond(A) = {np.linalg.cond(A):.1f}\n") print_corr("|corr(mixture x_i, source s_j)|：", abs_corr(X, S), "x", kinds) cov_err = np.max(np.abs(Z.T @ Z / N - np.eye(3))) print(f"\n PCA 白化后 Z 的协方差和 I 的最大差 = {cov_err:.1e}") Cz = abs_corr(Z, S) print_corr("|corr(whitened z_i, source s_j)|：", Cz, "z", kinds) print(f""" 观察 Observations： · 每个麦克风（x_i）都同时听到 3 个源，和每个源都有明显的相关。 · 白化（= PCA + 把每个主方向缩放到方差 1）之后，各分量互不相关、方差为 1 —— 二阶统计量已经"处理干净"了。 但每个 z_i 仍然是源的混合（每行最大的相关只有 {fmt(Cz.max(axis=1), '.2f')}，不是 1）。 · 原因：白化之后再乘以任何旋转矩阵 R，协方差仍然是 R I Rᵀ = I。所有旋转在二阶统计量上一模一样， PCA 只是其中一个"碰巧"的旋转（按方差排序）。剩下的这个旋转只能靠**高阶统计量 / 非高斯性**确定 —— 演示 2。 · 好消息：白化把 d² 个未知数减到了一个 d×d 正交矩阵（d(d−1)/2 个自由度），所以实践中 ICA 几乎总是先白化。""") # ============================================================================= # 演示 2：似然作为旋转角的函数 # ============================================================================= def rot(theta): """2×2 旋转矩阵。""" c, s = np.cos(theta), np.sin(theta) return np.array([[c, -s], [s, c]]) def demo_rotation_likelihood(): print_header("演示 2：after whitening only a rotation is left —— the logistic-prior likelihood vs rotation angle") pairs = {"gauss+gauss": ["gauss", "gauss"], "laplace+laplace": ["laplace", "laplace"], "sine+sawtooth": ["sine", "sawtooth"]} n = 20000 rng = np.random.default_rng(52) srcs = {name: make_sources(n, kinds, rng) for name, kinds in pairs.items()} angles = list(range(0, 91, 15)) print(f"S 是两个单位方差、独立的源（n = {n}）。W = c·R(θ)，c = π/√3（logistic 先验的标准差）。") print("表中是 CS229 ICA 的平均对数似然 (1/n) Σᵢ [Σ_j log g'(w_jᵀ s⁽ⁱ⁾)] + log|det W|。θ = 0° 和 90° 都是正确的分离\n") print(f"{'theta':>6}" + "".join(f"{name:>18}" for name in pairs)) table = {name: [] for name in pairs} for deg in angles: W = SCALE * rot(np.deg2rad(deg)) row = f"{deg:>5}°" for name in pairs: v = ica_logistic_loglik(W, srcs[name]) table[name].append(v) row += f"{v:>18.4f}" print(row) print(f"{'range':>6}" + "".join(f"{max(v) - min(v):>18.4f}" for v in table.values())) best = {name: angles[int(np.argmax(v))] for name, v in table.items()} print(f""" 观察 Observations： · 高斯源：七个角度的似然只差 {max(table['gauss+gauss']) - min(table['gauss+gauss']):.4f}（有限样本的波动）。 独立的标准高斯向量的分布是旋转对称的：R s 和 s 的分布**完全相同**，所以任何旋转都同样"正确"。 这就是为什么 ICA 要求源不是高斯的（最多允许一个高斯源）。也是 FA（代码 03）只能把 Λ 确定到一个旋转的原因。 · Laplace 源（super-Gaussian，尖峰重尾，像语音）：似然在 {best['laplace+laplace']}°（或 90°）最大，45° 最小。 · 正弦 + 锯齿波（sub-Gaussian，"两头堆积"）：似然在 {best['sine+sawtooth']}° 最大！在这个模型下，**错误的答案似然更高** —— 演示 3 的梯度上升会老老实实地找到它。 · 为什么？看 0° 这一行：Laplace {table['laplace+laplace'][0]:.3f} > 高斯 {table['gauss+gauss'][0]:.3f} > 正弦 + 锯齿波 {table['sine+sawtooth'][0]:.3f}。 方差固定之后，logistic 的 log 密度在尾部像 −|y|，所以每个分量越"尖"（越 super-Gaussian），平均 log 似然越高 —— 固定的 logistic 先验其实是在把每个分量推向"最尖"。而混合会让分布更接近高斯（中心极限定理的直觉）： Laplace 混合后变得不那么尖，似然下降；sub-Gaussian 信号本来比高斯"更平"，混合后反而变尖，似然上升。 · 所以"假设的源分布"不是一个无关紧要的细节：它至少要和真实源的峰度方向（super / sub）一致。""") # ============================================================================= # 演示 3：CS229 的 sigmoid ICA # ============================================================================= def demo_fit_ica(): print_header("演示 3：likelihood-based ICA (CS229 sigmoid prior) —— works for super-Gaussian sources only") cases = [("(a) 3 Laplace sources, logistic prior", ["laplace", "laplace", "laplace"], "logistic"), ("(b) sine/sawtooth/laplace, logistic prior", ["sine", "sawtooth", "laplace"], "logistic"), ("(c) sine/sawtooth/laplace, switching prior", ["sine", "sawtooth", "laplace"], "switching")] print("先白化，再从 W = I 出发做 3000 步批量梯度上升（学习率 0.1）。|corr| 矩阵：行 = 恢复的分量，列 = 真实的源\n") out = {} for title, kinds, prior in cases: S, A, X, Z, W_white = mixed(kinds, 53) hist = [] W = fit_ica(Z, prior=prior, history=hist) Y = Z @ W.T C = abs_corr(Y, S) P = W @ W_white @ A # 从源到恢复分量的总变换：y = W W_white A s Pn = P / np.max(np.abs(P), axis=1, keepdims=True) ll = ica_logistic_loglik(W, Z) out[title] = (C, ll, Y.std(axis=0)) print(title) print_corr("|corr(recovered y_i, source s_j)|：", C, "y", kinds) print(f" W·W_white·A（每行除以绝对值最大的元素）= " + " / ".join(fmt(r, '.2f') for r in Pn)) print(f" logistic 先验下的平均对数似然：初始 {hist[0]:.4f} → 最终 {ll:.4f}；恢复分量的标准差 = {fmt(Y.std(axis=0), '.2f')}\n") Ca, Cb, Cc = (out[t][0] for t, _, _ in cases) llb, llc = out[cases[1][0]][1], out[cases[2][0]][1] print(f"""观察 Observations： · (a) 3 个 super-Gaussian 源：|corr| 矩阵是一个排列矩阵（每个源的最大 |corr| = {fmt(Ca.max(axis=0), '.3f')}）。 W·W_white·A 每行只有一个 ±1，其余接近 0 —— 恰好是"排列 × 缩放"。两个不确定性是 ICA 本质的： 交换源的顺序、或者把一个源乘以 c 同时把 A 的对应列除以 c，观测 x 完全不变，所以没有任何方法能确定它们。 恢复分量的标准差也不是 1（{fmt(out[cases[0][0]][2], '.2f')}），而是由 logistic 先验的"偏好尺度"决定。 · (b) 同一个算法用在正弦 / 锯齿波 / Laplace 上：Laplace 被分离出来（{Cb.max(axis=0)[2]:.3f}），正弦和锯齿波的最大 |corr| 只有 {Cb.max(axis=0)[0]:.3f} / {Cb.max(axis=0)[1]:.3f} ≈ 1/√2 = 0.707 —— 正好是演示 2 预言的 45° 混合，W·W_white·A 里对应的两行是 (±1, ±1, 0)。 这不是没收敛：(b) 的最终似然 {llb:.4f} **高于** (c) 那个正确解在同一个 logistic 模型下的似然 {llc:.4f}。模型错了，MLE 就"正确地"找到了错误的答案。 · (c) 把源分布换成按峰度切换的版本（extended infomax：super-Gaussian 分量用尖峰的密度，sub-Gaussian 分量用双峰的密度， 每步根据当前估计自动选择），三个源都分离了（{fmt(Cc.max(axis=0), '.3f')}）。 · CS229 讲义的 sigmoid 在鸡尾酒会上好用，是因为**语音是 super-Gaussian 的**（大部分时间接近 0，偶尔很响）。 正弦这类"规则"的信号是 sub-Gaussian 的，需要别的先验 —— 和第 1 讲"选损失 = 选噪声模型"是同一个教训。""") # ============================================================================= # 演示 4：高斯源不可识别 # ============================================================================= def rotate_rows(W, i, j, phi): """把 W 的第 i、j 两行在它们张成的平面里旋转 φ（|det W| 不变，只改变这两个分量的边缘分布）。""" W2 = W.copy() W2[[i, j]] = rot(phi) @ W[[i, j]] return W2 def demo_gaussian_sources(): print_header("演示 4：Gaussian sources are not identifiable —— 1 Laplace + 2 Gaussian sources") kinds = ["laplace", "gauss", "gauss"] S, A, X, Z, _ = mixed(kinds, 54) print("(a) 同一份数据，从 5 个随机正交矩阵出发跑 sigmoid ICA，记录 3000 步和 10000 步时每个源的最大 |corr|\n") print(f"{'init':>5}{'steps':>7}{'log-lik':>11}" + "".join(f"{k:>10}" for k in kinds)) lls, cors, Ws = [], [], [] for r in range(5): W = np.linalg.qr(np.random.default_rng(540 + r).normal(size=(3, 3)))[0] for steps, extra in [(3000, 3000), (10000, 7000)]: W = fit_ica(Z, W0=W, iters=extra) # 接着上一次的 W 继续跑，总步数 = steps C = abs_corr(Z @ W.T, S) lls.append(ica_logistic_loglik(W, Z)) cors.append(C.max(axis=0)) print(f"{r:>5}{steps:>7}{lls[-1]:>11.5f}" + "".join(f"{v:>10.3f}" for v in cors[-1])) Ws.append(W) cors = np.array(cors) final = cors[1::2] # 每个初始化 10000 步时的结果 # (b) 换 6 份独立的数据（同一个 A，重新抽样 Laplace 和高斯源），每份都从 W = I 跑 10000 步 print("\n(b) 6 份独立抽样的数据（同一个混合矩阵 A），每份从 W = I 跑 10000 步：\n") print(f"{'data':>5}{'kurt g1':>9}{'kurt g2':>9}" + "".join(f"{k:>10}" for k in kinds)) cors_b = [] for r in range(6): Sb = make_sources(N, kinds, np.random.default_rng(600 + r)) Zb, _, _ = whiten(Sb @ A.T) Wb = fit_ica(Zb, iters=10000) cors_b.append(abs_corr(Zb @ Wb.T, Sb).max(axis=0)) kb = np.mean(Sb[:, 1:] ** 4, axis=0) print(f"{r:>5}{kb[0]:>9.3f}{kb[1]:>9.3f}" + "".join(f"{v:>10.3f}" for v in cors_b[-1])) cors_b = np.array(cors_b) # (c) 在 (a) 最后一个解处：把两个"高斯"分量在它们的平面里旋转 vs 把 Laplace 分量和一个高斯分量一起旋转 W = Ws[-1] C = abs_corr(Z @ W.T, S) i_lap = int(np.argmax(C[:, 0])) i_g = [i for i in range(3) if i != i_lap] phis = np.deg2rad(np.arange(0, 91, 15)) flat = [ica_logistic_loglik(rotate_rows(W, i_g[0], i_g[1], p), Z) for p in phis] steep = [ica_logistic_loglik(rotate_rows(W, i_lap, i_g[0], p), Z) for p in phis] print("\n(c) 在 (a) 的解处把两行一起旋转 φ = 0°, 15°, …, 90°，平均对数似然的变化范围：") print(f" 两个高斯分量的平面： {max(flat) - min(flat):.5f}") print(f" Laplace + 一个高斯分量： {max(steep) - min(steep):.5f}") print(f""" 观察 Observations： · Laplace 源在 (a)、(b) 的每一次都被干净地分离（最大 |corr| ≥ {min(cors[:, 0].min(), cors_b[:, 0].min()):.3f}）：只要一个源是非高斯的，它的方向就是确定的。 · (a) 同一份数据上，5 个初始化在 3000 步时给出不同的高斯分量（|corr| {cors[0::2, 1:].min():.3f}~{cors[0::2, 1:].max():.3f}）， 10000 步时都漂到了同一个解（|corr| {final[:, 1:].min():.3f}~{final[:, 1:].max():.3f}）；全部 10 个结果的对数似然只差 {max(lls) - min(lls):.5f}。 在这份样本里，高斯平面中确实有一个"最好"的旋转 —— 但它是从哪里来的？ · (b) 回答了这个问题：换一份独立抽样的数据，那个"最好"的旋转就换了一个地方，和真实高斯源的 |corr| 在 {cors_b[:, 1:].min():.2f} 到 {cors_b[:, 1:].max():.2f} 之间（45° 混合是 0.71，完全对齐是 1）。{N} 个高斯样本本身不是严格的高斯（样本峰度偏离 3）， 有限样本给了一点虚假的"偏好"，方向由抽样噪声决定，和源本身无关。n → ∞ 时这点偏好消失，任何旋转都一样好。 · (c) 在高斯平面里旋转，似然只变 {max(flat) - min(flat):.4f}；把 Laplace 分量拿去混合，似然变 {max(steep) - min(steep):.4f}。 高斯方向上的似然曲面几乎是平的，所以梯度上升在上面漂得很慢（(a) 里 3000 步还没停下）。 · 这就是"最多一个高斯源"的含义：ICA 能分离所有非高斯的源，把高斯部分作为一个整体子空间交出来， 但没法在子空间内部再分 —— 和 FA（代码 03）的 Λ 只能确定到一个旋转是同一回事。""") def main(): demo_cocktail_party() demo_rotation_likelihood() demo_fit_ica() demo_gaussian_sources() if __name__ == "__main__": main()