""" ============================================================================= 第 5 讲 · 代码 3：因子分析 —— d > n 时为什么需要"低秩 + 对角"的协方差 ============================================================================= 对应讲义：lectures/lecture05_unsupervised.md 第 4 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- CS229 讲 factor analysis 的动机是：特征数 d 比样本数 n 还多时，连"拟合一个高斯分布"都做不到。 演示 1：d = 30、n = 20。样本协方差 S 的秩只有 n − 1 = 19，是奇异的：高斯 MLE 的训练似然 = +∞， 测试似然 = −∞。加 εI 或只用对角线能"修好"，但各有代价。 演示 2：用 EM 拟合 factor analysis（k = 1~6 个因子）：对数似然单调不减；训练 / 测试似然； 因子数太多时出现 Heywood case（某个特征的噪声方差 Ψ_j → 0）。 演示 3：FA 学到了什么：Λ 只能确定到一个旋转（换一个旋转，似然逐位相同）；协方差矩阵的估计误差； 后验 E[z | x] 能多好地恢复真实的隐因子。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 03_factor_analysis.py # 约 1 秒 依赖本目录的 unsupervised.py。 """ import numpy as np from unsupervised import chol_logdet, fa_loglik, fit_fa_em, gaussian_logpdf, make_lowrank, print_header D, K_TRUE, N_TRAIN, N_TEST = 30, 3, 20, 5000 def fa_data(): """ 生成 FA 的数据：x = Λz + ε，d = 30 个特征由 k = 3 个因子驱动，每个特征有自己的噪声标准差 ∈ [0.5, 1.5]。 【返回】字典：Lam, psi（真噪声方差）, Sigma（真协方差 ΛΛᵀ + diag ψ）, X, Xte（训练 / 测试）, Zte（测试集的真因子） 【边界】所有演示共用同一份数据（固定种子 31），这样三个演示的数字可以互相对照。 """ rng = np.random.default_rng(31) Lam = rng.normal(size=(D, K_TRUE)) noise_std = rng.uniform(0.5, 1.5, size=D) X, _, _ = make_lowrank(N_TRAIN, Lam, noise_std, rng) Xte, Zte, _ = make_lowrank(N_TEST, Lam, noise_std, rng) psi = noise_std ** 2 return {"Lam": Lam, "psi": psi, "Sigma": Lam @ Lam.T + np.diag(psi), "X": X, "Xte": Xte, "Zte": Zte} def avg_ll(X, mu, Sigma): """每个样本的平均对数似然 (1/n) Σᵢ log N(x⁽ⁱ⁾; μ, Σ)。比总和好比较（训练集和测试集大小不同）。""" return float(np.mean(gaussian_logpdf(X, mu, Sigma))) # ============================================================================= # 演示 1：d > n 时高斯 MLE 失效 # ============================================================================= def demo_singular(): print_header("演示 1：d = 30 > n = 20 —— the Gaussian MLE has a singular covariance") data = fa_data() X, Xte = data["X"], data["Xte"] mu = X.mean(axis=0) S = (X - mu).T @ (X - mu) / len(X) # 高斯协方差的 MLE（除以 n） eig = np.linalg.eigvalsh(S) rank = int(np.sum(eig > 1e-10 * eig[-1])) try: chol_logdet(S) chol_msg = "Cholesky 成功（意外！）" except np.linalg.LinAlgError as err: chol_msg = f"Cholesky 失败：LinAlgError: {err}" print(f"样本协方差 S：{D}×{D}，秩 = {rank}（= n − 1：中心化用掉了一个自由度），" f"最小的 3 个特征值 = [{', '.join(f'{v:.1e}' for v in eig[:3])}]") print(f" {chol_msg}") print(f" → 高斯 MLE 的 log det S = −∞：训练似然 = +∞，而训练数据张成的 19 维子空间之外，密度为 0（测试似然 = −∞）\n") ll_truth = avg_ll(Xte, np.zeros(D), data["Sigma"]) print(f"{'model':<26}{'#params':>9}{'train ll/sample':>17}{'test ll/sample':>16}") rows = {} for eps in [1e-6, 1e-3, 1e-1, 0.3, 1.0, 3.0]: Se = S + eps * np.eye(D) rows[eps] = (avg_ll(X, mu, Se), avg_ll(Xte, mu, Se)) print(f"{'full Gaussian, S + ' + format(eps, 'g') + ' I':<26}{D + D * (D + 1) // 2:>9}" f"{rows[eps][0]:>17.2f}{rows[eps][1]:>16.2f}") Sd = np.diag(np.diag(S)) diag_tr, diag_te = avg_ll(X, mu, Sd), avg_ll(Xte, mu, Sd) print(f"{'diagonal Gaussian':<26}{2 * D:>9}{diag_tr:>17.2f}{diag_te:>16.2f}") print(f"{'true parameters':<26}{'-':>9}{avg_ll(X, np.zeros(D), data['Sigma']):>17.2f}{ll_truth:>16.2f}") best_eps = max(rows, key=lambda e: rows[e][1]) print(f""" 观察 Observations： · 20 个点（中心化后）只张成 19 维，30 维空间里剩下的 11 个方向上样本方差恰好是 0。 MLE 想让这些方向的方差 → 0，于是 log det → −∞、训练似然 → +∞：ε 从 1e-3 降到 1e-6， 训练似然从 {rows[1e-3][0]:.1f} 涨到 {rows[1e-6][0]:.1f}，测试似然从 {rows[1e-3][1]:.0f} 跌到 {rows[1e-6][1]:.0f}。 这是第 1 讲"XᵀX 奇异、d > n"在密度估计里的翻版。 · 加 εI（shrinkage / ridge）能修好，但 ε 很关键：这几个值里测试似然最好的是 ε = {best_eps:g}（{rows[best_eps][1]:.2f}）， 而 ε = 0.1 时只有 {rows[0.1][1]:.2f}。注意这里是**看着测试集**挑的 ε，是乐观的；真正使用时要交叉验证（第 3 讲）。 · 只用对角线（假设特征互相独立）永远可逆，参数也少（{2 * D} 个），但完全丢掉了特征之间的相关性： 测试似然 {diag_te:.2f}，比真参数的 {ll_truth:.2f} 差很多。 · 满协方差有 {D * (D + 1) // 2} 个参数，而训练集只有 {N_TRAIN}×{D} = {N_TRAIN * D} 个数。需要一个介于"对角"和"满"之间的结构 —— 演示 2。""") # ============================================================================= # 演示 2：EM 拟合 factor analysis # ============================================================================= def demo_fa_em(): print_header("演示 2：factor analysis by EM —— Σ = ΛΛᵀ + Ψ (low rank + diagonal)") data = fa_data() X, Xte = data["X"], data["Xte"] S_diag = X.var(axis=0) ll_truth = avg_ll(Xte, np.zeros(D), data["Sigma"]) mu = X.mean(axis=0) ridge_te = avg_ll(Xte, mu, (X - mu).T @ (X - mu) / len(X) + np.eye(D)) print(f"真实因子数 k = {K_TRUE}。每个 k 从同一个随机种子初始化，最多 5000 轮 EM。min Ψ_j/S_jj 到达下限 1e-4 = Heywood case\n") print(f"{'k':>3}{'#params':>9}{'EM iters':>10}{'monotone':>10}{'train ll/sample':>17}{'test ll/sample':>16}" f"{'min Psi_j/S_jj':>16}") res = {} for k in [1, 2, 3, 4, 6]: hist = [] mu_k, Lam, Psi, it = fit_fa_em(X, k, np.random.default_rng(32), max_iters=5000, history=hist) mono = all(b >= a - 1e-9 * abs(a) for a, b in zip(hist, hist[1:])) assert mono, f"EM decreased the FA log-likelihood (k={k})" tr = hist[-1] / len(X) te = float(np.mean(gaussian_logpdf(Xte, mu_k, Lam @ Lam.T + np.diag(Psi)))) ratio = float(np.min(Psi / S_diag)) n_par = D + D * k + D - k * (k - 1) // 2 res[k] = (tr, te, ratio, it) print(f"{k:>3}{n_par:>9}{it:>10}{str(mono):>10}{tr:>17.2f}{te:>16.2f}{ratio:>16.1e}") print(f"\n 对照（演示 1）：S + 1·I 的测试似然 {ridge_te:.2f}；真参数 {ll_truth:.2f}") best_k = max(res, key=lambda k: res[k][1]) heywood = [k for k in res if res[k][2] <= 1.0001e-4] print(f""" 观察 Observations： · 每个 k 的 EM 都让对数似然单调不减（assert 检查），和 GMM 的 EM 是同一个 ELBO 论证：E 步算 z | x 的高斯后验， M 步对 Λ 做一个"把 E[z] 当输入"的线性回归。 · 训练似然随 k 单调上升（模型越来越大），测试似然在 k = {best_k}（真实的因子数是 {K_TRUE}）最好：{res[best_k][1]:.2f}。 和第 1 讲的多项式次数、代码 01 的 k 一样，训练集上的指标不能用来选模型。 · k = {K_TRUE} 的 FA 只有 {D + D * K_TRUE + D - K_TRUE * (K_TRUE - 1) // 2} 个参数，测试似然和"看着测试集挑了 ε"的 S + 1·I（{ridge_te:.2f}）差不多， 远好于对角高斯。FA 的好处在于结构是**从数据里学出来的**：哪些特征共享哪些因子、每个特征自己的噪声有多大。 · k ∈ {heywood} 时出现了 Heywood case：某个 Ψ_j 被推到下限（这个特征被因子"完全"解释，自己的噪声方差 → 0）。 参数跑到了参数空间的边界上，EM 在边界附近收敛得极慢（这几个 k 都用满了 5000 轮）。k 太大时这是过拟合的信号； k 太小时也会发生 —— 一个常见的退化是模型把某个因子"专门"留给一个特征。只有 k = {K_TRUE} 干净地收敛（{res[K_TRUE][3]} 轮）。 Heywood case 通常提示模型设定或样本量有问题；它和 GMM 的塌缩都是"MLE 跑向边界"，只是这里有下限挡着。""") # ============================================================================= # 演示 3：FA 学到了什么 # ============================================================================= def demo_what_fa_learns(): print_header("演示 3：what FA learns —— rotation ambiguity, covariance error, recovering z") data = fa_data() X, Xte, Zte = data["X"], data["Xte"], data["Zte"] mu, Lam, Psi, _ = fit_fa_em(X, K_TRUE, np.random.default_rng(32), max_iters=5000) # (a) 旋转不确定性：Λ → ΛR（R 正交）不改变 ΛΛᵀ，所以似然完全相同 R, _ = np.linalg.qr(np.random.default_rng(33).normal(size=(K_TRUE, K_TRUE))) ll1, ll2 = fa_loglik(X, mu, Lam, Psi), fa_loglik(X, mu, Lam @ R, Psi) print(f"(a) 随机正交矩阵 R：ℓ(Λ) = {ll1:.6f}，ℓ(ΛR) = {ll2:.6f}，差 = {abs(ll1 - ll2):.1e}；" f"但 ‖Λ − ΛR‖_F = {np.linalg.norm(Lam - Lam @ R):.2f}") # (b) 协方差的估计误差（相对 Frobenius 范数） S = (X - X.mean(axis=0)).T @ (X - X.mean(axis=0)) / len(X) Sig = data["Sigma"] ests = {"sample covariance S": S, "S + 1·I": S + np.eye(D), "diagonal of S": np.diag(np.diag(S)), f"FA, k = {K_TRUE}": Lam @ Lam.T + np.diag(Psi)} print("\n(b) 协方差矩阵的估计误差 ‖Σ̂ − Σ‖_F / ‖Σ‖_F：") errs = {} for name, E in ests.items(): errs[name] = np.linalg.norm(E - Sig) / np.linalg.norm(Sig) print(f" {name:<22}{errs[name]:>8.3f}") off = Sig - np.diag(np.diag(Sig)) off_frac = np.linalg.norm(off) / np.linalg.norm(Sig) print(f" （真 Σ 的非对角部分占 ‖Σ‖_F 的比例 = {off_frac:.2f}）") psi_corr = np.corrcoef(Psi, data["psi"])[0, 1] print(f" FA 估计的 Ψ 和真实噪声方差的相关系数 = {psi_corr:.2f}") # (c) 后验均值 E[z | x] 作为"降维后的表示"：能多好地恢复真实因子（只能恢复到一个旋转，所以用回归的 R²） Sx = Lam @ Lam.T + np.diag(Psi) Ez = (Xte - mu) @ np.linalg.solve(Sx, Lam) # (n, k)：E[z | x] = Λᵀ Sx⁻¹ (x − μ) A = np.column_stack([np.ones(len(Ez)), Ez]) r2 = [] for j in range(K_TRUE): coef = np.linalg.lstsq(A, Zte[:, j], rcond=None)[0] r2.append(1 - np.var(Zte[:, j] - A @ coef) / np.var(Zte[:, j])) # 理论上限：用真参数算 E[z | x]，它和 z 的 R² 就是"给定 x，z 还剩多少不确定性" Ez_true = Xte @ np.linalg.solve(Sig, data["Lam"]) r2_true = [1 - np.var(Zte[:, j] - Ez_true[:, j]) / np.var(Zte[:, j]) for j in range(K_TRUE)] print(f"\n(c) 测试集上，用 E[z | x]（3 维）线性回归真实因子 z 的 R² = [{', '.join(f'{v:.3f}' for v in r2)}]" f"；用真参数算的 E[z | x] = [{', '.join(f'{v:.3f}' for v in r2_true)}]") print(f""" 观察 Observations： · (a) Λ 和 ΛR 给出相同的似然（差 = {abs(ll1 - ll2):.1g}），但 Λ 本身完全不同。所以 FA 的因子载荷 "只确定到一个旋转"：心理测量里常用 varimax 等方法再挑一个"好解释"的旋转。ICA（代码 05）要解决的正是 这个旋转不确定性 —— 前提是隐变量**不是**高斯的。 · (b) 只用 {N_TRAIN} 个样本，三种满协方差估计的 Frobenius 误差都很大：FA {errs[f'FA, k = {K_TRUE}']:.3f}、 样本协方差 {errs['sample covariance S']:.3f}、S + I {errs['S + 1·I']:.3f}，FA 只是略好。对角估计最差（{errs['diagonal of S']:.3f}）： 真 Σ 的非对角部分占了 {off_frac:.0%}，它全丢了。 注意误差度量会改变排序：S 的 Frobenius 误差不比 S + I 大，可它是奇异的、测试似然是 −∞。 Frobenius 范数只看"大的元素"，似然却对"小特征值的方向"极其敏感 —— 用哪个度量取决于你拿估计去做什么。 · (c) 后验均值 E[z | x] 就是 FA 给出的"低维表示"：它保留了真实因子的大部分信息（R² ≈ {np.mean(r2):.2f}， 用真参数也只有 ≈ {np.mean(r2_true):.2f} —— 噪声让 z 不能从 x 中完全确定）。 把 Ψ 限制成 σ²I 就得到 probabilistic PCA；再令 σ → 0，E[z | x] 就变成"投影到主方向上"—— 这就是 PCA（代码 04）。""") def main(): demo_singular() demo_fa_em() demo_what_fa_learns() if __name__ == "__main__": main()