""" ============================================================================= 第 3 讲 · 代码 4：正则化 —— ridge = MAP、lasso 与稀疏、early stopping 与数据增强 ============================================================================= 对应讲义：lectures/lecture03_generalization_regularization.md 第 4 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- 演示 1：第 1 讲的共线性房价数据（面积和卧室数相关 0.88）上加 ridge。用和第 1 讲完全相同的 2000 份数据： 卧室系数的标准差能降多少？代价（偏差）是什么方向的？包括第 1 讲那份"不走运"的样本。 演示 2：Ridge = 高斯先验下的 MAP。但 MAP 只是后验的一个点；完整的后验还给出**预测的不确定性**： 对"和训练数据相似"的房子很确定，对"面积大、卧室少"这种训练集里没有的组合很不确定。 演示 3：Lasso vs ridge：30 个相关特征里只有 5 个有用。lasso 的路径上系数一个个"恰好变成 0"，ridge 永远不会。 演示 4：两种"看起来不像正则化"的正则化： (a) early stopping —— 从 0 出发的 GD 跑 t 步 ≈ λ ≈ 1/(αt) 的 ridge； (b) 给输入加高斯噪声的数据增强 —— 期望上恰好等于 λ = s² 的 ridge。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 04_regularization.py # 约 3 秒 依赖本目录的 gen_tools.py。 """ import numpy as np from gen_tools import (SPARSE_TRUE_THETA, Standardizer, add_intercept, fit_lasso_cd, fit_least_squares, fit_ridge, lasso_kkt_violation, lasso_lambda_max, make_housing, make_sine, make_sparse_regression, poly_features, print_header, sine_target, to_unit_interval) HOUSING_TRUE = np.array([50.0, 0.12, 10.0]) # 第 1 讲 datasets.py 的真值 [截距, 面积 k$/ft², 卧室 k$/间] HOUSING_SIGMA = 40.0 def to_original_units(theta_std, sc): """ 把"标准化特征上的 θ"换算回原始单位（和第 1 讲 01_linear_regression.py 里的同名函数相同）。 【思路】x̃_j = (x_j − m_j)/s_j，所以 θ_j = θ̃_j / s_j，截距 θ₀ = θ̃₀ − Σ θ̃_j m_j / s_j。 """ slopes = theta_std[1:] / sc.std intercept = theta_std[0] - np.sum(theta_std[1:] * sc.mean / sc.std) return np.concatenate([[intercept], slopes]) def fit_housing(X, y, lam): """在标准化后的房价特征上拟合 ridge（不惩罚截距），返回 (原始单位的 θ, Standardizer)。""" sc = Standardizer().fit(X) th = fit_ridge(add_intercept(sc.transform(X)), y, lam) return to_original_units(th, sc), sc # ============================================================================= # 演示 1：ridge 遇上共线性 —— 先验对不对很重要 # ============================================================================= def ridge_repeated(true_bed, lams, reps=2000): """ 用第 1 讲的 2000 份房价数据（default_rng(10000 + rep)）重复拟合 ridge，统计卧室系数和预测误差。 【参数】 true_bed：卧室的真实系数（k$/间）。第 1 讲的真值是 10；为了做对照实验，可以换成别的值 —— 做法是把 make_housing 生成的 y 减去 10·bedrooms 再加上 true_bed·bedrooms（X 和噪声都不变）。 lams ：λ 列表；reps：重复次数 【返回】dict：λ → (卧室系数均值, 标准差, RMSE, 面积系数 RMSE, 平均 test MSE) test MSE 在 5000 套固定的测试房子（default_rng(4100)）上算，它们的 y 也按 true_bed 生成。 """ X_te, y_te = make_housing(5000, np.random.default_rng(4100)) y_te = y_te + (true_bed - HOUSING_TRUE[2]) * X_te[:, 1] ests = {lam: [] for lam in lams} tmse = {lam: [] for lam in lams} for rep in range(reps): X, y = make_housing(200, np.random.default_rng(10000 + rep)) y = y + (true_bed - HOUSING_TRUE[2]) * X[:, 1] for lam in lams: th, _ = fit_housing(X, y, lam) ests[lam].append(th) tmse[lam].append(np.mean((add_intercept(X_te) @ th - y_te) ** 2)) out = {} for lam in lams: E = np.array(ests[lam]) out[lam] = (float(E[:, 2].mean()), float(E[:, 2].std()), float(np.sqrt(np.mean((E[:, 2] - true_bed) ** 2))), float(np.sqrt(np.mean((E[:, 1] - HOUSING_TRUE[1]) ** 2))), float(np.mean(tmse[lam]))) return out def demo_ridge_collinear(): print_header("演示 1：共线性房价数据上的 ridge —— 和第 1 讲相同的 2000 份数据") lams = [0.0, 0.003, 0.01, 0.03, 0.1] X_ref, _ = make_housing(5000, np.random.default_rng(4100)) sd = X_ref.std(axis=0) print("每份数据 200 套房子；特征先标准化，λ 作用在标准化后的系数上（J = (1/2n)RSS + (λ/2)‖θ_{1:}‖²）。") print(f"面积的标准差 ≈ {sd[0]:.0f} ft²，卧室数的标准差 ≈ {sd[1]:.2f} 间。test MSE 在 5000 套新房子上算（噪声下限 σ² = 1600）。") print(f"(i) 第 1 讲的真值：面积 0.12、卧室 10 → 标准化后的真系数 ≈ [{0.12 * sd[0]:.0f}, {10 * sd[1]:.1f}]（很不平衡）") print(f"(ii) 对照实验：把卧室的真系数换成 80 → 标准化后 ≈ [{0.12 * sd[0]:.0f}, {80 * sd[1]:.0f}]（几乎相等）\n") print(f"{'truth':>6}{'lambda':>8}{'bed mean':>10}{'bed std':>9}{'bed RMSE':>10}{'area RMSE':>11}{'test MSE':>10}") res = {} for tb in [10.0, 80.0]: res[tb] = ridge_repeated(tb, lams) for lam in lams: r = res[tb][lam] print(f"{tb:>6g}{lam:>8g}{r[0]:>10.2f}{r[1]:>9.2f}{r[2]:>10.2f}{r[3]:>11.4f}{r[4]:>10.1f}") print() # 第 1 讲那份"不走运"的样本：default_rng(1) 生成的 200 套房子，OLS 把卧室系数估成 −3.3。 X1, y1 = make_housing(200, np.random.default_rng(1)) unlucky = {lam: fit_housing(X1, y1, lam)[0] for lam in lams} print(" 第 1 讲演示 1 那份不走运的样本（default_rng(1)，真值 10）的卧室系数：") print(" " + " ".join(f"λ={lam:g}: {unlucky[lam][2]:.2f}" for lam in lams)) a, b = res[10.0], res[80.0] best_b = min(lams, key=lambda lam: b[lam][2]) print(f""" 观察 Observations： · (i) λ = 0 一行复现了第 1 讲的重复实验（卧室平均 {a[0.0][0]:.2f}，标准差 {a[0.0][1]:.2f}）。加 ridge 后标准差确实下降 （λ = 0.1：{a[0.1][1]:.2f}），但卧室系数的**平均值往上跑**（{a[0.1][0]:.2f}），RMSE 从 {a[0.0][2]:.2f} 变成 {a[0.1][2]:.2f} —— 越正则越差。 只有很小的 λ（0.003：{a[0.003][2]:.2f}）勉强持平。预测误差也只升不降（{a[0.0][4]:.1f} → {a[0.1][4]:.1f}）。 · 为什么？ridge 不是把每个系数各自往 0 拉，而是沿着**数据方差小的方向**收缩。共线时方差最小的方向是 "两个标准化系数之差"，ridge 把两者往一起拉：大的（面积 ≈ 111）变小，小的（卧室 ≈ 15）变大。 用贝叶斯的话说，先验 θ_j ~ N(0, τ²) 独立同分布，相当于假设"两个特征差不多重要"—— 在这份数据上这个假设是错的。 · (ii) 对照实验里两个标准化系数几乎相等，先验"对了"：ridge 的偏差很小，RMSE 从 {b[0.0][2]:.2f} 降到 {b[best_b][2]:.2f}（λ = {best_b:g}）， 面积系数的 RMSE 也从 {b[0.0][3]:.4f} 降到 {b[best_b][3]:.4f}。λ 再大，整体往 0 收缩的偏差又占上风（λ = 0.1：{b[0.1][2]:.2f}）。 · 两组的改善都不大：n = 200、只有 2 个特征，方差本来就不是主要矛盾。ridge 真正大显身手的场景是 特征多、样本少（代码 01 的 k = 14、n = 15，代码 02 的 p ≈ n，演示 3 的 30 个特征 / 50 个样本）。 · 第 1 讲那份不走运的样本：ridge 把 −3.32 拉到 {unlucky[0.1][2]:.2f}（λ = 0.1）。看起来"修好了"，其实只是收缩的方向恰好 和这份样本的误差方向相反 —— 平均而言（上表）ridge 在这个问题上是往错的方向拉的。单个样本上的"成功"不能说明方法好。""") # ============================================================================= # 演示 2：MAP vs 完整后验 —— 预测的不确定性 # ============================================================================= def demo_posterior(): print_header("演示 2：Ridge = MAP；完整的后验还告诉你\"哪里不确定\"") rng = np.random.default_rng(42) X, y = make_housing(200, rng) n = len(y) sc = Standardizer().fit(X) A = add_intercept(sc.transform(X)) corr = float(np.corrcoef(X[:, 0], X[:, 1])[0, 1]) print(f"一份新的房价数据（n = {n}，面积和卧室数的相关系数 {corr:.2f}）。") print("模型：y = θᵀx̃ + N(0, σ²)，σ = 40 已知；先验：两个标准化系数 θ_j ~ N(0, τ²)，截距用平坦先验（不惩罚）。") print("后验是高斯： Σ_post = σ²(ÃᵀÃ + (σ²/τ²)D)⁻¹， μ_post = Σ_post Ãᵀy / σ² = ridge 解，λ = σ²/(nτ²)。") print("高斯后验的众数 = 均值，所以 MAP（ridge）就是后验均值；但后验还有协方差，它给出每个预测的不确定性。\n") houses = { "typical (2000 ft2, 3 bed)": [2000.0, 3.0], "big + few bedrooms (3800, 1)": [3800.0, 1.0], "small + many bedrooms (1000, 6)": [1000.0, 6.0], "extrapolate (6000 ft2, 8 bed)": [6000.0, 8.0], } D = np.diag([0.0, 1.0, 1.0]) print(f"{'house':<34}{'tau':>6}{'lambda':>9}{'mean':>9}{'sd of mean':>12}{'true f(x)':>11}") rows = {} for tau in [np.inf, 30.0, 10.0]: lam = 0.0 if np.isinf(tau) else HOUSING_SIGMA ** 2 / (n * tau ** 2) prec = A.T @ A / HOUSING_SIGMA ** 2 + (0.0 if np.isinf(tau) else 1.0 / tau ** 2) * D # 后验精度矩阵 cov = np.linalg.inv(prec) mu = cov @ A.T @ y / HOUSING_SIGMA ** 2 assert np.allclose(mu, fit_ridge(A, y, lam)) # 后验均值 = ridge 解（数值检查） for name, h in houses.items(): a = add_intercept(sc.transform(np.array([h])))[0] mean, sd = float(a @ mu), float(np.sqrt(a @ cov @ a)) truth = HOUSING_TRUE[0] + HOUSING_TRUE[1] * h[0] + HOUSING_TRUE[2] * h[1] rows[(name, tau)] = (mean, sd, truth) print(f"{name:<34}{tau:>6g}{lam:>9.3g}{mean:>9.1f}{sd:>12.1f}{truth:>11.1f}") print() t, b = rows[("typical (2000 ft2, 3 bed)", np.inf)], rows[("big + few bedrooms (3800, 1)", np.inf)] bt = rows[("big + few bedrooms (3800, 1)", 10.0)] e = rows[("extrapolate (6000 ft2, 8 bed)", np.inf)] print(f"""观察 Observations： · 平坦先验（τ = ∞，λ = 0，就是 OLS）下，典型房子的预测标准差只有 {t[1]:.1f}（k$），"面积大、卧室少"的房子是 {b[1]:.1f}， 外推到 6000 ft²、8 间是 {e[1]:.1f}。同一个模型，对不同输入的把握差好几倍 —— 训练集里几乎没有"大而卧室少"的房子 （面积和卧室相关 {corr:.2f}），系数在那个方向上没被数据约束住。MAP / 点估计把这个信息丢掉了。 · 注意外推到训练范围之外的房子（沿着"大房子卧室多"的主方向）反而比"3800 ft²、1 间"更确定： 决定不确定性的不是"离数据有多远"，而是"在数据没有约束的方向上走了多远"。 · τ 越小（先验越强，λ 越大），后验越窄（"大而卧室少"：{b[1]:.1f} → {bt[1]:.1f}），预测也被拉向先验。 τ = 10（λ ≈ {HOUSING_SIGMA ** 2 / (n * 100):.3g}）时典型房子几乎不受影响（{rows[("typical (2000 ft2, 3 bed)", 10.0)][0]:.1f}），"大而卧室少"的房子却被拉到 {bt[0]:.1f}（真值 {bt[2]:.1f}）， 偏了 {(bt[2] - bt[0]) / bt[1]:.1f} 个后验标准差：先验把两个系数往一起拉（演示 1），偏差恰好出现在训练集没覆盖的方向上。 **先验错了，后验会"自信地错"**：不确定性估计只在模型假设（包括先验）合理时才可信。 · 贝叶斯视角的价值不在于点估计（那就是 ridge），而在于**不确定性**：主动学习、Bayesian optimization、 以及"模型什么时候该说'我不知道'"都要用到它。代价是要对整个后验积分 —— 线性高斯模型能写出封闭解， 神经网络只能近似（Laplace 近似、变分推断、ensembles；第 4 讲的 bagging 是最朴素的版本）。""") # ============================================================================= # 演示 3：lasso vs ridge # ============================================================================= def demo_lasso(): print_header("演示 3：Lasso vs ridge —— 30 个相关特征，只有 5 个有用") rng = np.random.default_rng(43) X, y = make_sparse_regression(50, rng) X_te, y_te = make_sparse_regression(5000, rng) sc = Standardizer().fit(X) Xs, Xs_te = sc.transform(X), sc.transform(X_te) y_mean = y.mean() yc = y - y_mean # lasso 不拟合截距：先把 y 中心化 sup = [int(j) for j in np.flatnonzero(SPARSE_TRUE_THETA)] # 【Python 小知识】int(...) 把 np.int64 变成普通整数，打印更干净 true_support = set(sup) print(f"n = 50，d = 30，相邻特征相关系数 0.5；真实的非零系数在位置 {sup}，值 {SPARSE_TRUE_THETA[sup].tolist()}。") print("噪声 σ = 1，所以 test MSE 的下限是 1。λ 从 λ_max（解恰好全为 0）开始往下扫，每一步用上一步的解做 warm start。\n") lam_max = lasso_lambda_max(Xs, yc) lam_grid = lam_max * np.logspace(0, -3, 13) print("relaxed = 只用 lasso 选出的特征、重新做不加正则的 OLS（“用 lasso 选特征，用 OLS 估系数”）。\n") print(f"{'lambda':>9}{'nonzeros':>10}{'true pos':>10}{'false pos':>11}{'test MSE':>10}{'relaxed':>9}{'sweeps':>8}" f"{'KKT viol':>10}") theta = np.zeros(Xs.shape[1]) path = {} for lam in lam_grid: theta, sweeps = fit_lasso_cd(Xs, yc, lam, theta0=theta) nz = set(int(j) for j in np.flatnonzero(theta)) te = float(np.mean((y_mean + Xs_te @ theta - y_te) ** 2)) cols = sorted(nz) th_rx = fit_least_squares(add_intercept(Xs[:, cols]), y) # 没有选中任何特征时就是只有截距 te_rx = float(np.mean((add_intercept(Xs_te[:, cols]) @ th_rx - y_te) ** 2)) path[lam] = (theta.copy(), te, te_rx, len(nz), len(nz & true_support), len(nz - true_support)) print(f"{lam:>9.4f}{len(nz):>10}{len(nz & true_support):>10}{len(nz - true_support):>11}{te:>10.3f}{te_rx:>9.3f}" f"{sweeps:>8}{lasso_kkt_violation(Xs, yc, theta, lam):>10.1e}") best_lasso = min(path, key=lambda lam: path[lam][1]) best_rx = min(path, key=lambda lam: path[lam][2]) print("\n Ridge 路径（同样的数据；系数从不恰好为 0，所以数一下 |θ_j| > 0.05 的个数）：") print(f"{'lambda':>9}{'#|theta|>0.05':>15}{'exact zeros':>13}{'test MSE':>10}") ridge = {} Xi, Xi_te = add_intercept(Xs), add_intercept(Xs_te) for lam in [1e-3, 1e-2, 0.03, 0.1, 0.3, 1.0, 3.0]: th = fit_ridge(Xi, y, lam) te = float(np.mean((Xi_te @ th - y_te) ** 2)) ridge[lam] = (th, te) print(f"{lam:>9g}{np.sum(np.abs(th[1:]) > 0.05):>15}{np.sum(th[1:] == 0):>13}{te:>10.3f}") best_ridge = min(ridge, key=lambda lam: ridge[lam][1]) ols_te = float(np.mean((Xi_te @ fit_least_squares(Xi, y) - y_te) ** 2)) oracle = fit_least_squares(add_intercept(Xs[:, sup]), y) oracle_te = float(np.mean((add_intercept(Xs_te[:, sup]) @ oracle - y_te) ** 2)) th_l = path[best_lasso][0] print(f"\n test MSE 汇总（λ 都按测试集挑最好的，只为比较潜力）：OLS 全部 30 个特征 {ols_te:.3f}；" f"最好的 ridge {ridge[best_ridge][1]:.3f}（λ = {best_ridge:g}）；") print(f" 最好的 lasso {path[best_lasso][1]:.3f}（λ = {best_lasso:.4f}，{path[best_lasso][3]} 个非零）；" f"最好的 relaxed lasso {path[best_rx][2]:.3f}（λ = {best_rx:.4f}）；知道真实支撑集的 oracle OLS {oracle_te:.3f}。") th_true_std = SPARSE_TRUE_THETA * sc.std # 真值换到标准化坐标，方便和 lasso 系数对比 print("\n 最好的 lasso 在真实支撑集上的系数（标准化坐标）vs 真值：") print(" " + " ".join(f"j={j}: {th_l[j]:+.2f} (true {th_true_std[j]:+.2f})" for j in sup)) n_shrunk = sum(abs(th_l[j]) < abs(th_true_std[j]) for j in sup) print(f""" 观察 Observations： · Lasso 的系数一个个**恰好**变成 0（"exact zeros"），ridge 一个都没有（只是变小）：这是 L1 的尖角造成的（讲义 §4.3）。 λ 很大时 lasso 只保留最强的几个特征，而且先进来的都是真特征；λ 变小，假特征也陆续进来。 · 在这份数据上，最好的 lasso（{path[best_lasso][1]:.3f}）{'优于' if path[best_lasso][1] < ridge[best_ridge][1] else '不如'}最好的 ridge（{ridge[best_ridge][1]:.3f}），两者都远好于 OLS（{ols_te:.3f}）。 真相稀疏时 L1 的先验（Laplace）更对路；真相是"很多小系数"时 ridge 往往更好（练习 4）。 · 预测最好的 lasso 并不是"选得最准"的 lasso：它有 {path[best_lasso][3]} 个非零，其中 {path[best_lasso][5]} 个是假的。 为什么还要留着假特征？因为 lasso 用同一个 λ 同时做两件事 —— 选特征、收缩系数。λ 大到足以踢掉假特征时， 真特征的系数也被收缩得太多（软阈值把每个系数都往 0 挪 λ）；最好的 λ 是两者的折中。 5 个真特征里有 {n_shrunk} 个的系数比真值小（下方对比）。 · relaxed lasso 把两件事拆开：λ = {best_rx:.4f} 时 lasso 选出 {path[best_rx][3]} 个特征（{path[best_rx][4]} 真 {path[best_rx][5]} 假），再用 OLS 无收缩地估系数， test MSE {path[best_rx][2]:.3f}，{'好于' if path[best_rx][2] < path[best_lasso][1] else '不如'}最好的 lasso。但在 λ 小、选进很多假特征时 relaxed 反而更差 （λ = {best_lasso:.4f}：{path[best_lasso][2]:.3f} vs lasso {path[best_lasso][1]:.3f}）—— 没有收缩，假特征的 variance 全额进账。 · 知道真实支撑集的 oracle（{oracle_te:.3f}）是任何方法的天花板：只估 5 个参数，variance 最小。 "选对特征"的价值比"在错误的特征上调好正则化"大得多。 · 坐标下降的 KKT 违反量都在 1e−10 量级：每个 λ 都收敛到了真正的最优解。""") # ============================================================================= # 演示 4：early stopping 与噪声增强 —— 隐藏的 ridge # ============================================================================= def demo_hidden_ridge(): print_header("演示 4：(a) early stopping ≈ ridge；(b) 输入加噪声的数据增强 → ridge") rng = np.random.default_rng(44) # (a) Early stopping：k = 14、n = 15 的多项式（和代码 01 演示 2 相同的设置），R 个训练集同时跑 GD。 k, n, R = 14, 15, 200 x_te = rng.uniform(0.0, 1.0, size=1000) Phi_te = poly_features(to_unit_interval(x_te), k) f_te = sine_target(x_te) data = [make_sine(n, rng) for _ in range(R)] Phi = np.stack([poly_features(to_unit_interval(x), k) for x, _ in data]) # (R, n, p) Y = np.stack([y for _, y in data]) # (R, n) # 【numpy 小知识】np.stack 把 R 个 (n, p) 矩阵叠成 (R, n, p) 的三维数组； # Phi @ theta[:, :, None] 是"批量矩阵乘法"：对 R 个数据集各做一次 (n, p) @ (p, 1)，一行代码、没有 Python 循环。 G = np.einsum("rnp,rnq->rpq", Phi, Phi) / n # (R, p, p)，每个数据集的 ΦᵀΦ/n lr = 1.0 / np.linalg.eigvalsh(G)[:, -1] # (R,) 每个数据集自己的步长 1/λ_max theta = np.zeros((R, k + 1)) checkpoints = [10, 30, 100, 300, 1000, 3000, 10000, 30000] rows = {} for t in range(1, checkpoints[-1] + 1): resid = (Phi @ theta[:, :, None])[:, :, 0] - Y # (R, n) theta -= lr[:, None] * np.einsum("rnp,rn->rp", Phi, resid) / n if t in checkpoints: gd_te = np.mean((theta @ Phi_te.T - f_te) ** 2, axis=1) + 0.09 # (R,) gd_tr = np.mean(((Phi @ theta[:, :, None])[:, :, 0] - Y) ** 2, axis=1) lam = 1.0 / (lr * t) # 经验对应：λ ≈ 1/(α t) rd_te = [np.mean((Phi_te @ fit_ridge(Phi[r], Y[r], lam[r], intercept=False) - f_te) ** 2) + 0.09 for r in range(R)] rows[t] = (float(np.mean(gd_tr)), float(np.mean(gd_te)), float(np.median(lam)), float(np.mean(rd_te))) print(f"(a) k = {k}、n = {n} 的多项式回归，GD 从 θ = 0 出发、步长 1/λ_max，{R} 个训练集的平均。") print(" 和它比较的是 λ = 1/(α t) 的 ridge（这里连截距一起惩罚，因为 early stopping 对所有方向一视同仁）。\n") print(f"{'GD steps t':>11}{'train MSE':>11}{'GD test':>10}{'lambda=1/(a t)':>16}{'ridge test':>12}") for t in checkpoints: print(f"{t:>11}{rows[t][0]:>11.4f}{rows[t][1]:>10.4f}{rows[t][2]:>16.2e}{rows[t][3]:>12.4f}") best_t = min(checkpoints, key=lambda t: rows[t][1]) best_r = min(checkpoints, key=lambda t: rows[t][3]) # (b) 数据增强：在标准化的房价特征上加 N(0, s²) 噪声，复制 K 份，做普通最小二乘。 X, y = make_housing(200, np.random.default_rng(1)) # 第 1 讲那份不走运的样本 sc = Standardizer().fit(X) Xs = sc.transform(X) s = 0.3 th_ridge = fit_ridge(add_intercept(Xs), y, s ** 2) th_ols = fit_least_squares(add_intercept(Xs), y) print(f"\n(b) 第 1 讲那份房价数据（标准化特征），每个样本复制 K 份、每份加 N(0, {s}²) 的噪声，在 200K 行上做 OLS；") print(f" 对照：λ = s² = {s ** 2:g} 的 ridge。每个 K 重复 20 次（20 次独立的增强），系数是标准化坐标 [截距, 面积, 卧室]。\n") print(f"{'K copies':>9}{'mean theta over 20 draws':>30}{'sd of bed coef':>16}{'RMS |theta - ridge|':>21}") aug_rng = np.random.default_rng(45) aug = {} for K in [1, 10, 100, 1000]: ths = [] for _ in range(20): Xa = np.repeat(Xs, K, axis=0) + aug_rng.normal(0.0, s, size=(len(y) * K, 2)) ths.append(fit_least_squares(add_intercept(Xa), np.repeat(y, K))) ths = np.array(ths) aug[K] = (ths.mean(axis=0), float(ths[:, 2].std()), float(np.sqrt(np.mean(np.sum((ths - th_ridge) ** 2, axis=1))))) txt = "[" + ", ".join(f"{v:.2f}" for v in aug[K][0]) + "]" print(f"{K:>9}{txt:>30}{aug[K][1]:>16.3f}{aug[K][2]:>21.3f}") print(f"{'ridge':>9}{'[' + ', '.join(f'{v:.2f}' for v in th_ridge) + ']':>30}") print(f"{'OLS':>9}{'[' + ', '.join(f'{v:.2f}' for v in th_ols) + ']':>30}") print(f""" 观察 Observations： · (a) GD 的测试误差随步数 t 先降后升（t = {best_t} 时最好，{rows[best_t][1]:.4f}；30000 步时 {rows[30000][1]:.3g}）—— 和 ridge 随 λ 的 U 形是同一条曲线倒过来走：t 小 ≈ λ 大（欠拟合），t 大 ≈ λ 小（过拟合）。 对应关系 λ ≈ 1/(αt) 只是粗略的：ridge 一列的最好值 {rows[best_r][3]:.4f} 出现在 t = {best_r} 对应的 λ，比 GD 的最好值略好， 两条曲线的谷底错开了一点。原因：从 0 出发的 GD 在第 i 个奇异方向上拟合了 1 − (1 − α s_i²/n)^t 的比例， ridge 是 s_i²/(s_i² + nλ)；两者都是"大奇异值方向几乎全拟合、小奇异值方向几乎不拟合"的滤波器，形状相似但不相同。 · 这就是为什么"按验证误差提前停止训练"是深度学习里最常用、最便宜的正则化：它不用改损失函数，只要多看一眼验证集。 · (b) 增强后的 OLS 在 20 次抽样上的平均值和 λ = s² 的 ridge 几乎一样（卧室系数 {aug[1000][0][2]:.2f} vs {th_ridge[2]:.2f}；OLS 是 {th_ols[2]:.2f}）； 单次增强的随机波动随 K 按 1/√K 缩小（离 ridge 的 RMS 距离：K = 1 时 {aug[1][2]:.3f}，K = 1000 时 {aug[1000][2]:.3f}）。 原因是 E_ε[(y − (x + ε)ᵀθ)²] = (y − xᵀθ)² + s²‖θ_{{1:}}‖²：对输入加噪声，期望上就是给损失加了 L2 惩罚。 "让模型对输入的小扰动不敏感"和"惩罚 ‖θ‖"是一回事。Dropout、随机裁剪、mixup 都是这个思路的推广： 它们对应的正则项更复杂（和数据相关），但本质相同。 · 反过来也成立：gen_tools.fit_ridge 就是用"加 p 行假数据"实现 ridge 的 —— 正则化 = 往训练集里加入你的先验知识。""") def main(): demo_ridge_collinear() demo_posterior() demo_lasso() demo_hidden_ridge() if __name__ == "__main__": main()