""" ============================================================================= 第 4 讲 · 代码 3：Boosting —— AdaBoost 与 gradient boosting ============================================================================= 对应讲义：lectures/lecture04_trees_ensembles_nn.md 第 4 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题】 ----------------------------------------------------------------------------- Bagging 平均的是“一群独立的、强的、高方差的”模型，降方差；boosting 串行地叠加“一串弱的、高偏差的”模型， 每一个都去修补前面留下的错误，降偏差。它在优化什么？为什么不太容易过拟合？什么时候会过拟合？ 演示 1：AdaBoost = 指数损失上的前向分步加法建模（forward stagewise additive modeling）。数值验证三件事： 权重 w ∝ exp(−y·F)；α = ½ log((1 − ε)/ε) 就是精确的 line search；每轮的归一化常数 Z = 2√(ε(1 − ε))。 演示 2：AdaBoost + stumps 的训练曲线。训练误差 ≤ 指数损失 = Π 2√(ε(1 − ε))，指数下降；训练误差到 0 之后， 测试误差和 margin 怎么变？标签有噪声时又怎样？ 演示 3：Gradient boosting = 函数空间里的梯度下降。学习率（shrinkage）ν 的影响：ν 小 → 需要更多棵树，但泛化更好。 和 bagging 对比：boosting 降的是偏差。 演示 4：树的深度 = 交互的阶数。stumps 只能表示加法模型 Σ g_j(x_j)（数值验证：混合差分恰好为 0）。 演示 5：换损失 = 换伪残差。有离群点时，绝对值损失（伪残差 = sign(残差)）比平方损失稳得多 —— 第 1 讲 §2.4 的老故事。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 03_boosting.py # 约 12 秒 依赖本目录的 tree_models.py、nonlinear_data.py，以及第 1 讲的 linear_models.py、datasets.py。 """ import numpy as np from nonlinear_data import make_circle, make_friedman1, make_sine, sine_target from tree_models import AdaBoost, GradientBoostingRegressor, RandomForest from linear_models import mse, print_header def ternary_min(fun, lo, hi, iters=200): """ 在 [lo, hi] 上找一维凸函数 fun 的最小点（三分搜索，ternary search）。 【思路】每次比较两个三等分点的函数值，扔掉较差那一侧的 1/3。200 次后区间缩小到原来的 (2/3)^200 ≈ 1e−35， 远小于 float64 的精度。只适用于单峰（unimodal）函数 —— 指数损失沿一个方向是凸的，满足条件。 【返回】最小点（float） """ for _ in range(iters): m1, m2 = lo + (hi - lo) / 3.0, hi - (hi - lo) / 3.0 if fun(m1) < fun(m2): hi = m2 else: lo = m1 return 0.5 * (lo + hi) # ============================================================================= # 演示 1：AdaBoost = 指数损失上的前向分步加法建模 # ============================================================================= def demo_adaboost_is_fsam(): print_header("演示 1：AdaBoost 就是在指数损失上一步一步地加 stump（forward stagewise additive modeling）") rng = np.random.default_rng(31) X, y01 = make_circle(300, rng) y = 2 * y01 - 1 ab = AdaBoost(6).fit(X, y01) H = np.stack([2 * s.predict(X) - 1 for s in ab.stumps]) # (T, n)：每个 stump 的 ±1 预测 print("圆形边界数据，300 个点。逐轮检查（F_{t−1} 是前 t − 1 个 stump 的加权和，F_0 = 0）：") print(" (1) AdaBoost 的权重 w_t 和 exp(−y·F_{t−1}) 归一化之后是否相同；") print(" (2) 公式 α_t = ½ log((1 − ε_t)/ε_t) 和数值最小化 L(α) = Σ exp(−y·(F_{t−1} + α·h_t)) 得到的 α 是否相同；") print(" (3) 损失的缩小倍数 L(F_t)/L(F_{t−1}) 是否等于 2√(ε_t(1 − ε_t))。\n") print(f"{'t':>3}{'eps_t':>8}{'max|w - w_exp|':>16}{'alpha (formula)':>17}{'alpha (search)':>16}{'L_t/L_t-1':>11}" f"{'2*sqrt(e(1-e))':>16}") F = np.zeros(len(y)) worst_w, worst_a, worst_z = 0.0, 0.0, 0.0 for t in range(len(ab.alphas)): w_exp = np.exp(-y * F) w_exp /= w_exp.sum() dw = float(np.max(np.abs(ab.weights[t] - w_exp))) h = H[t] a_search = ternary_min(lambda a: float(np.sum(np.exp(-y * (F + a * h)))), 0.0, 5.0) L_old = float(np.sum(np.exp(-y * F))) F = F + ab.alphas[t] * h ratio = float(np.sum(np.exp(-y * F))) / L_old z = 2.0 * np.sqrt(ab.eps[t] * (1.0 - ab.eps[t])) worst_w, worst_a, worst_z = max(worst_w, dw), max(worst_a, abs(a_search - ab.alphas[t])), max(worst_z, abs(ratio - z)) print(f"{t + 1:>3}{ab.eps[t]:>8.3f}{dw:>16.1e}{ab.alphas[t]:>17.6f}{a_search:>16.6f}{ratio:>11.6f}{z:>16.6f}") print(f""" 观察 Observations： · (1) 权重的最大差 {worst_w:.0e}：AdaBoost “把分错的样本权重乘以 e^α” 的规则，累积下来恰好是 w_i ∝ exp(−y_i F(x_i))， 也就是“每个样本在当前指数损失里占的份额”。被反复分错的点 y·F 很负，权重指数级变大。 · (2) α 的公式和数值 line search 一致（最大差 {worst_a:.0e}）：令 dL/dα = 0，分对的权重和 (1 − ε) 乘 e^{{−α}}、 分错的 ε 乘 e^{{α}}，相等时 α = ½ log((1 − ε)/ε)。AdaBoost 最初（Freund & Schapire 1995/1997）是从在线学习（Hedge 算法）的角度提出的； Friedman, Hastie & Tibshirani（2000）才指出它等价于这个优化过程 —— 一个很好的“先有算法、后有解释”的例子。 · (3) 每一轮损失恰好乘以 Z_t = 2√(ε_t(1 − ε_t))（最大差 {worst_z:.0e}）。只要 ε_t < 0.5，Z_t < 1：损失严格下降。 stump 的选择也是“最小化 ε_t”= 最小化 Z_t：stump 和 α 都是贪心地让指数损失下降最多。""") # ============================================================================= # 演示 2：训练曲线、上界、margin # ============================================================================= def demo_adaboost_curves(): print_header("演示 2：AdaBoost + stumps 的训练曲线 —— 训练误差上界、测试误差、margin") rng = np.random.default_rng(32) X, y01 = make_circle(300, rng) Xt, yt01 = make_circle(5000, rng) T = 2000 ab = AdaBoost(T).fit(X, y01) y, yt = 2 * y01 - 1, 2 * yt01 - 1 F_tr, F_te = ab.staged_scores(X), ab.staged_scores(Xt) # (T, n) train_err = np.mean(np.where(F_tr >= 0, 1, -1) != y, axis=1) test_err = np.mean(np.where(F_te >= 0, 1, -1) != yt, axis=1) bound = np.cumprod(2.0 * np.sqrt(ab.eps * (1.0 - ab.eps))) # Π Z_t = 平均指数损失 gamma = 0.5 - ab.eps bound2 = np.exp(-2.0 * np.cumsum(gamma ** 2)) margins = y * F_tr / np.cumsum(ab.alphas)[:, None] # 归一化 margin ∈ [−1, 1] t0 = int(np.argmax(train_err == 0)) + 1 if np.any(train_err == 0) else None print(f"圆形边界（无标签噪声），训练 300 个点，测试 5000 个点，{T} 轮。单个 stump 的测试误差 = {test_err[0]:.3f}。") print(f"训练误差第一次等于 0 是在第 {t0} 轮。margin = y·F(x)/Σα ∈ [−1, 1]。\n") print(f"{'round':>6}{'eps_t':>7}{'train err':>11}{'exp-loss bound':>16}{'exp(-2 sum g^2)':>17}{'test err':>10}" f"{'min margin':>12}{'median margin':>15}") rounds = sorted(set([1, 5, 10, 20, 50, 100, 200, 500, 1000, 2000] + ([t0] if t0 else []))) for r in rounds: i = r - 1 print(f"{r:>6}{ab.eps[i]:>7.3f}{train_err[i]:>11.3f}{bound[i]:>16.4f}{bound2[i]:>17.4f}{test_err[i]:>10.4f}" f"{margins[i].min():>12.3f}{np.median(margins[i]):>15.3f}") best = int(np.argmin(test_err)) + 1 # 有标签噪声的版本 rng = np.random.default_rng(33) Xn, yn = make_circle(300, rng, flip=0.1) Xnt, ynt = make_circle(5000, rng, flip=0.1) abn = AdaBoost(T).fit(Xn, yn) tr_n = np.mean((abn.staged_scores(Xn) >= 0).astype(int) != yn, axis=1) te_n = np.mean((abn.staged_scores(Xnt) >= 0).astype(int) != ynt, axis=1) best_n = int(np.argmin(te_n)) + 1 print("\n同样的设置，但随机翻转 10% 的标签（Bayes error = 0.10）：\n") print(f"{'round':>6}{'train err':>11}{'test err':>10}") for r in sorted(set([1, 10, 50, 100, 200, 500, 1000, 2000, best_n])): print(f"{r:>6}{tr_n[r - 1]:>11.3f}{te_n[r - 1]:>10.4f}") print(f""" 观察 Observations： · “exp-loss bound”一列就是 Π_t 2√(ε_t(1 − ε_t))，它等于训练集上的平均指数损失，而 1{{yF ≤ 0}} ≤ exp(−yF)， 所以它是训练误差的上界；再用 1 − x ≤ e^{{−x}} 得到更松的 exp(−2Σγ_t²)（γ_t = 0.5 − ε_t 是“比随机猜好多少”）。 只要每一轮都有 γ_t ≥ γ > 0（弱学习器假设），训练误差就指数级地趋于 0。 · ε_t 越来越接近 0.5（第 {rounds[-1]} 轮 {ab.eps[-1]:.3f}）：后面的 stump 面对的是被前面的轮次反复加权的“难”样本， 越来越难找到好的 stump。但只要不等于 0.5，损失就还在下降。 · 训练误差在第 {t0} 轮第一次到 0（注意它不是单调的：第 100 轮又是 {train_err[99]:.3f}，单调下降的只有上界）。 之后测试误差没有立刻回升：第 {t0} 轮 {test_err[t0 - 1]:.4f}，第 {best} 轮最低 {test_err[best - 1]:.4f}。同时最小 margin 从 {margins[t0 - 1].min():.3f} 增大到 {margins[-1].min():.3f}：训练点都已经分对了，AdaBoost 还在继续把它们推离边界。 Schapire et al.（1998）用 margin 理论解释这个“训练误差为 0 后继续变好”的现象。 但改善很小，而且后期测试误差又回升到 {test_err[-1]:.4f}（比第 {t0} 轮还差）—— “AdaBoost 不会过拟合”并不成立， 只是过拟合得慢；轮数仍然要用验证集选。 · 有 10% 标签噪声时：测试误差在第 {best_n} 轮最低（{te_n[best_n - 1]:.4f}），之后 AdaBoost 为了把翻转的标签也分对， 给它们越来越大的权重（指数损失对错得离谱的点惩罚是指数级的），训练误差在第 {T} 轮降到 {tr_n[-1]:.3f}，测试误差升到 {te_n[-1]:.4f}。 指数损失对标签噪声很敏感（Dietterich 2000）—— 和第 1 讲平方损失怕离群点是同一类问题。实践中用 early stopping、 shrinkage，或者换成 logistic loss（LogitBoost / gradient boosting）。""") # ============================================================================= # 演示 3：Gradient boosting 与 shrinkage # ============================================================================= def demo_shrinkage(): print_header("演示 3：Gradient boosting —— 学习率（shrinkage）ν 的影响；和 bagging 比较") rng = np.random.default_rng(34) X, y = make_friedman1(300, rng) Xt, yt = make_friedman1(5000, rng) M = 1000 print(f"Friedman #1（d = 10），训练 300 个点，测试 5000 个点（噪声方差 1 是测试 MSE 的下限）。深度 3 的树，{M} 轮。\n") print(f"{'nu':>6}{'test@10':>9}{'test@100':>10}{'test@300':>10}{'test@1000':>11}{'best test':>11}{'best round':>12}" f"{'train@1000':>12}") res = {} for nu in [1.0, 0.3, 0.1, 0.03]: gb = GradientBoostingRegressor(M, nu, max_depth=3).fit(X, y) te = np.mean((gb.staged_predict(Xt) - yt) ** 2, axis=1) tr = np.mean((gb.staged_predict(X) - y) ** 2, axis=1) b = int(np.argmin(te)) res[nu] = (te[b], b + 1, te[-1], tr[-1]) print(f"{nu:>6g}{te[9]:>9.2f}{te[99]:>10.2f}{te[299]:>10.2f}{te[-1]:>11.2f}{te[b]:>11.2f}{b + 1:>12}{tr[-1]:>12.4f}") # ν = 0.03 在 1000 轮时可能还没走完：再给它 3 倍的轮数看看 gb_long = GradientBoostingRegressor(3 * M, 0.03, max_depth=3).fit(X, y) te_long = np.mean((gb_long.staged_predict(Xt) - yt) ** 2, axis=1) b_long = int(np.argmin(te_long)) print(f"\n补充：ν = 0.03 跑 {3 * M} 轮，最好的测试 MSE = {te_long[b_long]:.2f}（第 {b_long + 1} 轮），第 {3 * M} 轮 {te_long[-1]:.2f}") bag = mse(yt, RandomForest(200, "regression", seed=35).fit(X, y).predict(Xt)) print(f"对照：同一份数据上 bagging（200 棵不限深度的树）的测试 MSE = {bag:.2f}") print(f""" 观察 Observations： · 每一轮：算伪残差 r = y − F（平方损失的负梯度），用一棵深度 3 的树拟合 r，再 F ← F + ν·tree。 这就是 GD 的 θ ← θ − α∇J，只不过“参数”是函数 F 在训练点上的取值，“梯度”只能用一棵树来近似。 · ν = 1（不收缩）：第 {res[1.0][1]} 轮就达到最好的 {res[1.0][0]:.2f}，之后很快过拟合（{M} 轮时 {res[1.0][2]:.2f}）。每棵树都把残差（包括噪声） 一口吃掉，几轮之后剩下的“残差”基本就是噪声了。 · ν 越小，需要的轮数越多（最好的轮数：ν = 0.3 → {res[0.3][1]}，0.1 → {res[0.1][1]}，0.03 → {res[0.03][1]}）。 ν 从 1 降到 0.3、0.1，最好的测试误差明显下降（{res[1.0][0]:.2f} → {res[0.3][0]:.2f} → {res[0.1][0]:.2f}）；ν = 0.03 在 {M} 轮时还在下降 （{res[0.03][0]:.2f}），给它 {3 * M} 轮，最好是 {te_long[b_long]:.2f} —— 和 ν = 0.1 几乎一样。收益递减：再小的 ν 只是更慢。 Friedman（2001）的经验法则：ν 取小（≤ 0.1），轮数用验证集 / early stopping 选。 直觉：小步长让很多棵树“分摊”每一处修正，相当于一种正则化，类似于 GD 里的小学习率 + early stopping。 · 训练误差到 {M} 轮时几乎为 0（ν = 0.1：{res[0.1][3]:.4f}），但 ν ≤ 0.3 时测试误差几乎不回升（ν = 0.1：最好 {res[0.1][0]:.2f}，{M} 轮 {res[0.1][2]:.2f}）。 过拟合是有的，只是很慢。 · boosting（{res[0.1][0]:.2f}）明显好于 bagging（{bag:.2f}）。代码 2 说过 bagging 后剩下的主要是偏差；boosting 用很多棵小树 逐步修正残差，降的正是偏差。这也是为什么在表格数据上 gradient boosting（XGBoost、LightGBM）通常是最强的基线。""") # ============================================================================= # 演示 4：树的深度 = 交互的阶数 # ============================================================================= def mixed_difference(model, rng, j=0, k=1, n=500, d=10): """ 特征 j、k 之间“交互作用”的大小：E|F(a, b) − F(a, b') − F(a', b) + F(a', b')|（其他特征固定在随机值上）。 【思路】如果 F 是加法的，F = g_j(x_j) + g_k(x_k) + 其他，那么这个混合差分恰好为 0（四项两两抵消）。 只要它不为 0，模型里就有 x_j 和 x_k 的交互项（比如 x_j·x_k）。 【返回】float """ base = rng.uniform(0.0, 1.0, size=(n, d)) a, a2, b, b2 = (rng.uniform(0.0, 1.0, size=n) for _ in range(4)) def F(u, v): Z = base.copy() Z[:, j], Z[:, k] = u, v return model.predict(Z) return float(np.mean(np.abs(F(a, b) - F(a, b2) - F(a2, b) + F(a2, b2)))) def demo_depth_interactions(): print_header("演示 4：每棵树的深度 = 能表示的交互阶数（stumps → 加法模型）") rng = np.random.default_rng(34) X, y = make_friedman1(300, rng) # 和演示 3 同一份数据（同一个种子） Xt, yt = make_friedman1(5000, rng) M = 1000 print(f"同演示 3 的数据，ν = 0.1，{M} 轮。真实函数含交互项 10·sin(π x₁ x₂)。") print("最后一列：x₁、x₂ 的混合差分 E|F(a,b) − F(a,b') − F(a',b) + F(a',b')|，加法模型恰好为 0。\n") print(f"{'max_depth':>10}{'best test':>11}{'best round':>12}{'test@1000':>11}{'x1-x2 interaction':>19}") res = {} for depth in [1, 2, 3, 6]: gb = GradientBoostingRegressor(M, 0.1, max_depth=depth).fit(X, y) te = np.mean((gb.staged_predict(Xt) - yt) ** 2, axis=1) b = int(np.argmin(te)) inter = mixed_difference(gb, np.random.default_rng(36)) res[depth] = (te[b], b + 1, te[-1], inter) print(f"{depth:>10}{te[b]:>11.2f}{b + 1:>12}{te[-1]:>11.2f}{inter:>19.2e}") best_depth = min(res, key=lambda k: res[k][0]) print(f""" 观察 Observations： · stumps（深度 1）的模型是 Σ_m (一个只看某一个特征的阶梯函数) = Σ_j g_j(x_j)：一个加法模型（additive model）。 混合差分是 {res[1][3]:.0e}（浮点误差级别），证明它确实表示不了任何交互；所以学不好 sin(π x₁ x₂)，最好只有 {res[1][0]:.2f}。 · 深度 2 的树可以表示两两交互，混合差分变成 {res[2][3]:.2f}，测试 MSE 降到 {res[2][0]:.2f}。 Friedman #1 最高只有两个特征之间的交互，所以深度 {best_depth} 最好；深度 6 的树太灵活（每棵树都在拟合残差里的噪声），反而变差（{res[6][0]:.2f}）。 · 经验：gradient boosting 用浅树（深度 2~8，或 4~32 个叶子）；random forest 用深树。 前者靠“很多棵弱学习器的和”降偏差，后者靠“很多棵强学习器的平均”降方差 —— 零件的选择正好相反。 · 演示 2 的圆形边界 x₁² + x₂² < r² 是加法的，所以 AdaBoost + stumps 能做好；换成 XOR 型的边界（x₁·x₂ > 0），stumps 就无能为力了 （见 README 的“动手改一改”）。""") # ============================================================================= # 演示 5：换损失 = 换伪残差（离群点） # ============================================================================= def demo_robust_loss(): print_header("演示 5：换损失 = 换伪残差 —— 有离群点时用绝对值损失（第 1 讲 §2.4 的老故事）") rng = np.random.default_rng(35) x, y = make_sine(200, rng) n_out = 20 idx = rng.choice(200, size=n_out, replace=False) y = y.copy() y[idx] += rng.choice([-1.0, 1.0], size=n_out) * rng.uniform(3.0, 6.0, size=n_out) # 20 个离群点，偏 ±3~6 xt = np.linspace(0.0, 1.0, 1001) ft = sine_target(xt) print(f"y = sin(2πx) + N(0, 0.3²)，200 个训练点，其中 {n_out} 个加了 ±3~6 的大误差。深度 2 的树，ν = 0.1。") print("表中是和真实函数 sin(2πx) 的 MSE（不含噪声）。\n") print(f"{'loss':<10}{'pseudo-residual':<18}{'@10':>8}{'@50':>8}{'@100':>8}{'@200':>8}{'@500':>8}{'best':>8}{'best round':>12}") res = {} for loss, pr in [("squared", "y − F"), ("absolute", "sign(y − F)")]: gb = GradientBoostingRegressor(500, 0.1, max_depth=2, loss=loss).fit(x[:, None], y) e = np.mean((gb.staged_predict(xt[:, None]) - ft) ** 2, axis=1) b = int(np.argmin(e)) res[loss] = (e[b], b + 1, e[-1]) print(f"{loss:<10}{pr:<18}" + "".join(f"{e[r - 1]:>8.3f}" for r in [10, 50, 100, 200, 500]) + f"{e[b]:>8.3f}{b + 1:>12}") print(f""" 观察 Observations： · Gradient boosting 的框架和损失无关：每一轮拟合“损失对预测值的负梯度”。平方损失 → 残差 y − F； 绝对值损失 → sign(y − F)（只看方向，不看大小）；logistic loss → y − σ(F)（又是第 1 讲 GLM 的“目标 − 预测”）。 · 平方损失下，离群点的残差很大，后面的树会专门去拟合它们：和真实函数的 MSE 最好只有 {res['squared'][0]:.3f}， {500} 轮时恶化到 {res['squared'][2]:.3f}。绝对值损失下，每个点的伪残差只有 ±1，离群点没有特殊的拉力：最好 {res['absolute'][0]:.3f}，{500} 轮 {res['absolute'][2]:.3f}。 · 和第 1 讲 §2.4 完全相同的道理：平方损失 = 高斯噪声假设，薄尾，怕离群点；绝对值损失 = Laplace 噪声，厚尾，稳健。 选损失仍然是在选噪声模型 —— 不管模型是线性的、树的和，还是神经网络。""") def main(): demo_adaboost_is_fsam() demo_adaboost_curves() demo_shrinkage() demo_depth_interactions() demo_robust_loss() if __name__ == "__main__": main()