""" ============================================================================= 第 5 讲 · 代码 4：PCA —— 最大方差 = 最小重构误差、怎样选 k、标准化什么时候重要 ============================================================================= 对应讲义：lectures/lecture05_unsupervised.md 第 5 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- 演示 1：低秩 + 噪声的数据（d = 10，真实秩 3）。特征值谱：前 3 个是"信号 + 噪声"，其余是平坦的"噪声地板" σ²。 SVD 和协方差特征分解给出同样的结果；按"解释方差比例"选 k 的几个常见阈值分别给出什么。 演示 2：重构误差恰好等于丢掉的特征值之和（逐个 k 数值验证）。训练 / 测试的重构误差随 k 单调下降， 所以**不能**用它选 k；和"干净信号"的误差才是 U 形的（PCA 用于去噪）。 演示 3：PCA 不是尺度不变的。(A) 一个单位很大的无关特征霸占了第一主成分；(B) 同一单位下的 "几乎是常数"的噪声特征被标准化放大。标准化不是永远正确的默认选项。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 04_pca.py # 不到 1 秒 依赖本目录的 unsupervised.py。 """ import numpy as np from unsupervised import fit_pca, make_lowrank, pca_reconstruct, pca_transform, print_header SIGNAL_SCALES = np.array([3.0, 2.0, 1.0]) # 三个隐因子方向上的标准差（信号强度） NOISE_STD = 0.5 def lowrank_data(): """ d = 10、真实秩 3 的数据：Λ 的三列是正交的单位向量乘以 [3, 2, 1]，再加各向同性噪声 N(0, 0.5²)。 【思路】这样真实协方差 ΛΛᵀ + σ²I 的特征值是精确已知的：[9.25, 4.25, 1.25, 0.25 × 7]。 【返回】(X, Xte, clean_te, Lam)：训练 500 个、测试 2000 个，以及测试集的"干净信号" """ rng = np.random.default_rng(41) Q, _ = np.linalg.qr(rng.normal(size=(10, 3))) # 10×3，列是正交的单位向量 Lam = Q * SIGNAL_SCALES X, _, _ = make_lowrank(500, Lam, NOISE_STD, rng) Xte, _, clean_te = make_lowrank(2000, Lam, NOISE_STD, rng) return X, Xte, clean_te, Lam # ============================================================================= # 演示 1：特征值谱与选 k # ============================================================================= def demo_spectrum(): print_header("演示 1：the eigenvalue spectrum of low-rank + noise data (d = 10, true rank 3)") X, _, _, Lam = lowrank_data() model = fit_pca(X) lam = model["eigvals"] Xc = X - X.mean(axis=0) lam_eig = np.linalg.eigvalsh(Xc.T @ Xc / len(X))[::-1] # 另一种算法：协方差矩阵的特征值（从大到小） theory = np.concatenate([SIGNAL_SCALES ** 2 + NOISE_STD ** 2, np.full(7, NOISE_STD ** 2)]) evr = lam / lam.sum() cum = np.cumsum(evr) print(f"n = {len(X)}。theory = 真实协方差的特征值 [3², 2², 1²] + σ²（σ = {NOISE_STD}）\n") print(f"{'j':>3}{'eigval (SVD)':>14}{'eigval (eigh)':>15}{'theory':>9}{'explained':>11}{'cumulative':>12}") for j in range(10): print(f"{j + 1:>3}{lam[j]:>14.4f}{lam_eig[j]:>15.4f}{theory[j]:>9.2f}{evr[j]:>11.3f}{cum[j]:>12.3f}") # 前 3 个主方向张成的子空间 vs Λ 的列空间：最大主角（principal angle）的余弦 U = model["components"][:3].T # 10×3，正交的列 Qt, _ = np.linalg.qr(Lam) cos_min = np.linalg.svd(U.T @ Qt, compute_uv=False).min() picks = {thr: int(np.searchsorted(cum, thr) + 1) for thr in [0.8, 0.95, 0.99]} ratios = lam[:-1] / lam[1:] gap_k = int(np.argmax(ratios) + 1) print(f"\n SVD vs eigh 的特征值最大差 = {np.max(np.abs(lam - lam_eig)):.1e}") print(f" 前 3 个主方向的子空间 vs 真实的 span(Λ)：最大主角的余弦 = {cos_min:.4f}（1 = 完全重合）") print(" 按累计解释方差选 k：" + "，".join(f"≥{thr:.0%} → k = {k}" for thr, k in picks.items()) + f"；相邻特征值之比最大的位置 → k = {gap_k}（λ_{gap_k}/λ_{gap_k + 1} = {ratios[gap_k - 1]:.2f}）") print(f""" 观察 Observations： · 前 3 个特征值 ≈ 信号方差 + σ²，后 7 个都挤在 σ² = {NOISE_STD ** 2} 附近（{lam[3:].min():.3f}~{lam[3:].max():.3f}）： 这就是"噪声地板"。地板不是完全平的 —— 有限样本让本来相等的特征值散开（随机矩阵理论的 Marchenko–Pastur 现象）， d/n 越大散得越开。 · SVD 和协方差特征分解给出同样的特征值（差 {np.max(np.abs(lam - lam_eig)):.0e}）。实践中用 SVD：不用构造 XᵀX（条件数是 X 的平方）， 而且 d ≫ n 时可以只算前几个奇异向量（randomized SVD）。 · 前 3 个主方向和真实的信号子空间几乎重合（余弦 {cos_min:.4f}）。 · "解释 X% 的方差"是最常用的选 k 规则，但阈值是任意的：80% 给 k = {picks[0.8]}（丢了一个真实的因子）， 95% 给 k = {picks[0.95]}、99% 给 k = {picks[0.99]}（把噪声方向也收进来了）—— 7 个噪声方向合起来占了 {lam[3:].sum() / lam.sum():.0%} 的方差。 · 看谱的"断崖"更可靠：相邻特征值之比 λ₁/λ₂ = {ratios[0]:.2f}、λ₂/λ₃ = {ratios[1]:.2f}、λ₃/λ₄ = {ratios[2]:.2f}、 之后都 ≤ {ratios[3:].max():.2f}。最大的比值在 j = {gap_k}，断崖之后是平坦的地板（scree plot 的"肘部"）。 但它也不是万能的：如果最弱的信号因子只比噪声稍强，它和地板之间的比值可能小于强因子之间的比值。 更有原则的做法是用概率模型（PPCA / FA）的 held-out 似然选 k（代码 03 演示 2）。""") # ============================================================================= # 演示 2：重构误差 = 丢掉的特征值之和 # ============================================================================= def demo_reconstruction(): print_header("演示 2：max variance = min reconstruction error —— error equals the sum of discarded eigenvalues") X, Xte, clean_te, _ = lowrank_data() model = fit_pca(X) lam = model["eigvals"] total = lam.sum() print("误差都是每个样本的平均平方误差 (1/n) Σᵢ ‖x⁽ⁱ⁾ − x̂⁽ⁱ⁾‖²；PCA 只用训练集拟合\n") print(f"{'k':>3}{'train error':>13}{'sum discarded':>15}{'difference':>12}{'captured+error':>16}" f"{'test error':>12}{'error vs clean':>16}") rows = [] for k in range(11): err = float(np.mean(np.sum((X - pca_reconstruct(X, model, k)) ** 2, axis=1))) disc = float(lam[k:].sum()) captured = float(np.mean(np.sum((pca_transform(X, model, k)) ** 2, axis=1))) R = pca_reconstruct(Xte, model, k) err_te = float(np.mean(np.sum((Xte - R) ** 2, axis=1))) err_clean = float(np.mean(np.sum((clean_te - R) ** 2, axis=1))) rows.append((k, err, disc, err_te, err_clean)) print(f"{k:>3}{err:>13.6f}{disc:>15.6f}{err - disc:>12.1e}{captured + err:>16.6f}{err_te:>12.4f}{err_clean:>16.4f}") noisy_vs_clean = float(np.mean(np.sum((Xte - clean_te) ** 2, axis=1))) best = min(rows, key=lambda r: r[4]) print(f"\n 总方差 Σλ_j = {total:.6f}；测试集上“什么都不做”（直接用带噪声的 x）和干净信号的误差 = {noisy_vs_clean:.4f}") print(f""" 观察 Observations： · 第 2、3 列相同（差 ≤ {max(abs(r[1] - r[2]) for r in rows):.0e}，舍入误差）：用前 k 个主成分重构，平均平方误差**恰好**等于丢掉的特征值之和； "保留的方差 + 重构误差"恒等于总方差 {total:.4f}（第 5 列）。所以"最大化投影方差"和"最小化重构误差"是同一个问题 —— 这是勾股定理：‖x‖² = ‖投影‖² + ‖残差‖²。 · 训练误差和测试误差都随 k **单调下降**，k = 10 时为 0：更大的子空间总能把点投影得更近。 所以"在验证集上算重构误差"选不出 k（和代码 01 的 distortion 一样）。需要一个概率模型（PPCA / FA 的 held-out 似然） 或者别的准则。 · 和**干净信号**的误差是 U 形的：k = {best[0]} 时最小（{best[4]:.4f}），比直接用带噪声的数据（{noisy_vs_clean:.4f}）好得多。 丢掉的方向主要是噪声，PCA 起到了去噪（denoising）的作用；k 再大，就把噪声也"重构"回来了。 这一列在真实问题里算不出来（不知道干净信号），这里是为了说明"PCA 保留了什么"。""") # ============================================================================= # 演示 3：标准化 # ============================================================================= def factor_r2(F, Z): """ 每个真实因子 Z[:, j] 能被 F 的线性组合解释多少（带截距的线性回归 R²）。 【思路】PCA 的主成分和真实因子之间差一个旋转（以及可能的符号、顺序），不能逐个比较； 用回归的 R² 衡量"F 张成的子空间里包含了多少 z_j 的信息"，这个度量和坐标、单位无关。 【返回】长度 Z.shape[1] 的数组 """ A = np.column_stack([np.ones(len(F)), F]) out = [] for j in range(Z.shape[1]): coef = np.linalg.lstsq(A, Z[:, j], rcond=None)[0] out.append(1.0 - np.var(Z[:, j] - A @ coef) / np.var(Z[:, j])) return np.array(out) def demo_standardization(): print_header("演示 3：PCA is not scale-invariant —— when to standardize") rng = np.random.default_rng(43) L = rng.normal(size=(5, 2)) # 5 个有信息的特征，由 2 个隐因子驱动，噪声标准差 0.5 X5, Z, _ = make_lowrank(500, L, 0.5, rng) junk_big = rng.normal(size=(500, 1)) * 50.0 # (A) 一个和因子无关、但单位很大的特征（标准差 50） junk_tiny = rng.normal(size=(500, 5)) * 0.01 # (B) 5 个同一单位下几乎是常数的噪声特征（标准差 0.01） scenarios = { "oracle: 5 informative only": X5, "(A) + 1 junk feature, std 50": np.column_stack([X5, junk_big]), "(B) + 5 junk features, std 0.01": np.column_stack([X5, junk_tiny]), } print("5 个有信息的特征由 2 个隐因子 z 驱动；R² = 用前 2 个主成分线性回归 z₁、z₂ 的 R²（越接近 1 越好）\n") print(f"{'data':<34}{'standardize':>12}{'PC1 explains':>14}{'PC1 top feature':>17}{'R2 z1':>8}{'R2 z2':>8}") res = {} for name, Xs in scenarios.items(): for std in [False, True]: m = fit_pca(Xs, standardize=std) r2 = factor_r2(pca_transform(Xs, m, 2), Z) top = int(np.argmax(np.abs(m["components"][0]))) res[(name, std)] = (m["eigvals"][0] / m["eigvals"].sum(), top, r2) if name.startswith("(A)") and std: share = m["components"][:, 5] ** 2 # 无关特征 x6 的"载荷平方"在各个主成分上的分布（每列平方和 = 1） junk_pc, junk_share = int(np.argmax(share)) + 1, float(share.max()) print(f"{name:<34}{str(std):>12}{res[(name, std)][0]:>14.3f}{'x' + str(top + 1):>17}" f"{r2[0]:>8.3f}{r2[1]:>8.3f}") # 换单位：把第 1 个特征从"米"换成"毫米"（×1000）。标准化的 PCA 应该完全不受影响。 X_mm = X5 * np.array([1000.0, 1, 1, 1, 1]) e1, e2 = fit_pca(X5, standardize=True)["eigvals"], fit_pca(X_mm, standardize=True)["eigvals"] u1, u2 = fit_pca(X5)["eigvals"], fit_pca(X_mm)["eigvals"] print(f"\n 换单位（x1 × 1000）：标准化 PCA 的特征值变化 = {np.max(np.abs(e1 - e2)):.1e}；" f"不标准化时 PC1 解释的比例从 {u1[0] / u1.sum():.3f} 变成 {u2[0] / u2.sum():.6f}") A0, A1 = res[("(A) + 1 junk feature, std 50", False)], res[("(A) + 1 junk feature, std 50", True)] B0, B1 = res[("(B) + 5 junk features, std 0.01", False)], res[("(B) + 5 junk features, std 0.01", True)] O0 = res[("oracle: 5 informative only", False)] O1 = res[("oracle: 5 informative only", True)] snr = np.sum(L ** 2, axis=1) # 每个有信息特征的信号方差 ‖Λ_j‖²（噪声方差都是 0.25） print(f""" 观察 Observations： · PCA 找的是"方差最大"的方向，而方差取决于单位。(A) 那个无关特征的标准差是 50，不标准化时 PC1 几乎就是它 （解释了 {A0[0]:.1%} 的方差，最大载荷在 x{A0[1] + 1}），前 2 个主成分对真实因子的 R² 只有 {A0[2][0]:.2f} / {A0[2][1]:.2f}。 标准化后 R² 回到 {A1[2][0]:.2f} / {A1[2][1]:.2f}，几乎和"只用 5 个有信息的特征、同样标准化"（{O1[2][0]:.2f} / {O1[2][1]:.2f}）一样： 标准化后无关特征的方差也是 1，但它和别的特征都不相关，几乎单独构成了第 {junk_pc} 个主成分（{junk_share:.1%} 的载荷在那里）。 · 换单位的检验更直接：把 x1 乘以 1000，不标准化的 PCA 完全变了（PC1 几乎全是 x1），标准化的 PCA 一点不变。 特征的单位不同（米 vs 元 vs 人数）时，方差之间根本不可比，**应该**标准化（用相关矩阵做 PCA）。 · 但 (B) 是反例：同一单位下，几个"几乎是常数"的特征本来对 PCA 毫无影响（不标准化：R² = {B0[2][0]:.2f} / {B0[2][1]:.2f}）， 标准化把它们的方差从 0.0001 放大到 1，和有信息的特征平起平坐，R² 掉到 {B1[2][0]:.2f} / {B1[2][1]:.2f}。 图像的像素、同一种传感器的读数、词频 —— 单位相同、方差大小本身有意义时，标准化反而有害。 · 连没有无关特征的 oracle 数据上，标准化也让 R² 从 {O0[2][0]:.2f} / {O0[2][1]:.2f} 降到 {O1[2][0]:.2f} / {O1[2][1]:.2f}：5 个特征的噪声相同（0.5²）， 信号方差却差很多（{', '.join(f'{v:.2f}' for v in snr)}），方差大 = 信噪比高。不标准化时这些特征权重大，恰好是对的； 标准化把它们拉平了，几乎没有信号的 x{int(np.argmin(snr)) + 1} 被放大。 · 规则：单位不同 → 标准化；单位相同且方差有意义 → 不标准化（只中心化）。拿不准时两种都试，看下游任务。""") def main(): demo_spectrum() demo_reconstruction() demo_standardization() if __name__ == "__main__": main()