""" ============================================================================= 第 3 讲 · 代码 3：学习理论 —— Hoeffding、union bound、一致收敛与 VC 维 ============================================================================= 对应讲义：lectures/lecture03_generalization_regularization.md 第 3 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题 What this file answers】 ----------------------------------------------------------------------------- "训练误差小 ⇒ 测试误差小"什么时候成立？学习理论的回答是：当训练误差对假设类里的**所有**假设 同时是测试误差的好估计（uniform convergence）时。这里用能精确计算真实误差的玩具问题检验每一步： 演示 1：Hoeffding 不等式 vs 精确的二项分布尾概率：上界对不对？松多少？为什么松？ 演示 2：有限假设类（k 个网格阈值）：sup_h |训练误差 − 真实误差| 的实际分布 vs 一致收敛的上界 γ = sqrt(log(2k/δ)/(2m))；ERM 的超额误差 vs 2γ。k 从 2 加到 10⁴，上界一直涨，实际却不涨 —— 因为阈值函数的"有效大小"由 VC 维（= 1）决定，不由网格点数决定。 演示 3：VC 维与 shattering：平面上的直线（半平面）能实现 n 个点的多少种标注？ 3 个点：8/8（shattered）；4 个点：14/16（XOR 做不到）；n 个点：多项式增长（Sauer 引理），不是 2ⁿ。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 03_learning_theory.py # 约 1 秒 依赖本目录的 gen_tools.py。只用到精确公式和简单的模拟，没有训练任何模型。 """ from math import comb, exp, lgamma, log import numpy as np from gen_tools import (hoeffding_bound, make_noisy_threshold, print_header, sample_complexity, threshold_true_error, uniform_gamma) # ============================================================================= # 演示 1：Hoeffding vs 精确的二项分布尾概率 # ============================================================================= def binomial_two_sided_tail(m, phi, gamma): """ 精确计算 P(|φ̂ − φ| > γ)，φ̂ = (1/m)·Binomial(m, φ)。 【思路】把所有满足 |k/m − φ| > γ 的 k 的二项概率加起来。C(m, k) 在 m 大时是天文数字， 所以在对数空间里算：log P(k) = lgamma(m+1) − lgamma(k+1) − lgamma(m−k+1) + k log φ + (m−k) log(1−φ)。 【参数】m：样本数；phi ∈ (0, 1)；gamma > 0 【返回】float 【边界】|k − φm| 恰好等于 γm 的 k 不算（严格大于）；用 1e−9 的容差，避免 0.6 − 0.5 = 0.0999999… 这类浮点误差 把边界上的 k 错误地算进去或漏掉。 【Python 小知识】math.lgamma(n + 1) = log(n!)，对很大的 n 也不会溢出。 """ total = 0.0 for k in range(m + 1): if abs(k - phi * m) > gamma * m + 1e-9: total += exp(lgamma(m + 1) - lgamma(k + 1) - lgamma(m - k + 1) + k * log(phi) + (m - k) * log(1 - phi)) return total def demo_hoeffding(): print_header("演示 1：Hoeffding 上界 2exp(−2γ²m) vs 精确的尾概率（γ = 0.1）") gamma = 0.1 ms = [10, 25, 50, 100, 200, 500, 1000] print("φ̂ 是 m 个 Bernoulli(φ) 的平均（比如：一个固定的分类器在 m 个测试样本上的错误率）。") print("表中是 P(|φ̂ − φ| > 0.1) 的精确值（二项分布直接求和）和 Hoeffding 上界。\n") print(f"{'m':>6}{'exact phi=0.5':>15}{'exact phi=0.1':>15}{'Hoeffding':>12}{'bound/exact(0.5)':>18}") rows = {} for m in ms: e5 = binomial_two_sided_tail(m, 0.5, gamma) e1 = binomial_two_sided_tail(m, 0.1, gamma) hb = float(hoeffding_bound(m, gamma)) rows[m] = (e5, e1, hb) print(f"{m:>6}{e5:>15.3g}{e1:>15.3g}{hb:>12.3g}{hb / e5:>18.3g}") print("\n 样本复杂度：要以 95% 的把握让 |φ̂ − φ| ≤ γ，Hoeffding 说需要 m ≥ log(2/δ)/(2γ²)：") for g in [0.1, 0.05, 0.01]: print(f" γ = {g:<5g} → m ≥ {sample_complexity(1, g, 0.05):>6}") print(f""" 观察 Observations： · 上界永远成立（Hoeffding 列 ≥ 两个精确列），而且和精确值一样随 m **指数**下降 —— 这是它有用的原因。 · 但它很松：m = 10、25 时上界大于 1（{rows[10][2]:.2f}、{rows[25][2]:.2f}），没有任何信息；m = 1000 时上界是精确值的 {rows[1000][2] / rows[1000][0]:.0f} 倍。 原因是 Hoeffding 只用了"每个变量在 [0, 1] 之间"，必须对最坏的分布也成立；它没有用到方差。 · φ = 0.1 时精确尾概率小得多（m = 200：{rows[200][1]:.1e} vs φ = 0.5 的 {rows[200][0]:.1e}），因为方差 φ(1−φ) 只有 0.09 而不是 0.25。 用上方差的不等式（Bernstein）会紧得多；学习理论里用 Hoeffding 是为了简单和通用。 · 样本复杂度按 1/γ² 增长：精度提高 10 倍，样本要多 100 倍。这也是"测试集要多大"的经验法则： 想把错误率估到 ±1%（95% 置信），需要上万个测试样本。""") # ============================================================================= # 演示 2：有限假设类的一致收敛；ERM 的保证 # ============================================================================= def threshold_train_errors(x, y, grid): """ 对网格上的每个阈值 t，计算 h_t(x) = 1{x > t} 在训练集上的错误率 —— 一次算完所有 t。 【思路】先把 x 排序。对阈值 t，令 j = #{x_i ≤ t}（np.searchsorted 给出）。h_t 把前 j 个点预测为 0、其余预测为 1： 错误数 = （前 j 个点里 y = 1 的个数） + （后 m − j 个点里 y = 0 的个数） 用前缀和 c[j] = 前 j 个点里 y = 1 的个数，两项都是 O(1)：c[j] + (m − j) − (c[m] − c[j])。 【参数】x, y：训练集（y 是 0/1）；grid：(k,) 阈值 【返回】(k,) 每个阈值的训练错误率 【numpy 小知识】np.searchsorted(xs, grid, side="right") 对每个 t 返回"xs 里 ≤ t 的元素个数"（xs 必须已排序）。 np.concatenate([[0], np.cumsum(ys)]) 是带一个前导 0 的前缀和，c[j] = ys[:j].sum()。 """ order = np.argsort(x) xs, ys = x[order], y[order] m = len(xs) c = np.concatenate([[0], np.cumsum(ys)]) j = np.searchsorted(xs, grid, side="right") return (c[j] + (m - j) - (c[m] - c[j])) / m def uniform_deviation_experiment(m, k, R, rng, delta=0.05): """ 对网格阈值类 H_k（k 个等间距的阈值），重复 R 次"采 m 个训练样本 → 看所有假设的训练误差和真实误差"。 【返回】dict： sup_dev：(R,) 每次的 max_{h ∈ H_k} |ε̂(h) − ε(h)| excess ：(R,) 每次 ERM 选出的假设的真实误差 − H_k 里最好的真实误差 gamma ：一致收敛上界 γ = sqrt(log(2k/δ)/(2m)) """ grid = np.linspace(0.0, 1.0, k) true_err = threshold_true_error(grid) sup_dev = np.empty(R) excess = np.empty(R) for r in range(R): x, y = make_noisy_threshold(m, rng) tr = threshold_train_errors(x, y, grid) sup_dev[r] = np.max(np.abs(tr - true_err)) excess[r] = true_err[np.argmin(tr)] - true_err.min() # argmin 遇到并列时取第一个 return {"sup_dev": sup_dev, "excess": excess, "gamma": float(uniform_gamma(m, k, delta))} def demo_uniform_convergence(): print_header("演示 2：有限假设类 —— 一致收敛的上界 vs 实际（阈值分类器，δ = 0.05）") rng = np.random.default_rng(52) R = 2000 print("数据：x ~ U(0,1)，y = 1{x > 0.3}，再以 10% 的概率翻转（Bayes error = 0.1）。") print("H_k = { h_t(x) = 1{x > t} : t 取 [0, 1] 上 k 个等间距的点 }。真实误差有精确公式 ε(t) = 0.1 + 0.8|t − 0.3|。") print(f"每个设置重复 R = {R} 次。sup dev = max_h |ε̂(h) − ε(h)|；excess = ε(ERM 选的 h) − min_h ε(h)。\n") m = 100 print(f"(a) m = {m}，k 变大：") print(f"{'k':>7}{'gamma bound':>13}{'P(sup dev>gamma)':>18}{'95% sup dev':>13}{'mean excess':>13}" f"{'95% excess':>12}{'2*gamma':>9}") rows_k = {} for k in [2, 11, 101, 1001, 10001]: res = uniform_deviation_experiment(m, k, R, rng) sd, ex, g = res["sup_dev"], res["excess"], res["gamma"] rows_k[k] = (g, float(np.mean(sd > g)), float(np.quantile(sd, 0.95)), float(np.mean(ex))) print(f"{k:>7}{g:>13.3f}{np.mean(sd > g):>18.4f}{np.quantile(sd, 0.95):>13.3f}{np.mean(ex):>13.4f}" f"{np.quantile(ex, 0.95):>12.3f}{2 * g:>9.3f}") print(f"\n(b) k = 1001，m 变大：") print(f"{'m':>7}{'gamma bound':>13}{'95% sup dev':>13}{'ratio':>8}{'mean excess':>13}") rows_m = {} for m2 in [25, 100, 400, 1600]: res = uniform_deviation_experiment(m2, 1001, R, rng) q = float(np.quantile(res["sup_dev"], 0.95)) rows_m[m2] = (res["gamma"], q, float(np.mean(res["excess"]))) print(f"{m2:>7}{res['gamma']:>13.3f}{q:>13.3f}{res['gamma'] / q:>8.2f}{np.mean(res['excess']):>13.4f}") print(f""" 观察 Observations： · 上界成立：在所有设置下，sup dev 超过 γ 的频率都远小于 δ = 0.05（最大的一个是 {max(r[1] for r in rows_k.values()):.4f}）， ERM 的超额误差的 95% 分位也都小于 2γ。 · (a) γ 随 log k 一直增长：k 从 2 到 10001，γ 从 {rows_k[2][0]:.3f} 涨到 {rows_k[10001][0]:.3f}。实际的 95% sup dev 从 {rows_k[2][2]:.3f} 涨到 {rows_k[101][2]:.3f}（k = 101）后 就基本停住了（k = 1001：{rows_k[1001][2]:.3f}，k = 10001：{rows_k[10001][2]:.3f}）。union bound 把 k 个假设当成互不相关的坏事件来加； 但相邻的阈值几乎做同样的预测，它们的偏差高度相关，"有效"的假设个数远小于 k。 这正是 VC 维要解决的问题：阈值类不管取多少个网格点、甚至取遍所有实数，VC 维都是 1（演示 3 讲 VC 维）。 · (b) m 每乘以 4，γ 和实际的 sup dev 都大约减半（1/√m）：上界的**速率**是对的，只是常数松了约 {min(r[0] / r[1] for r in rows_m.values()):.1f}～{max(r[0] / r[1] for r in rows_m.values()):.1f} 倍。学习理论的界通常就是这样用的：看它怎样随 m、模型大小变化，而不是看具体数值。 · ERM 的平均超额误差（m = 100 时约 {rows_k[1001][3]:.3f}）远小于 2γ：上界是"最坏情况 + 最坏的 5%"，平均情况好得多。""") # ============================================================================= # 演示 3：VC 维 —— 半平面能实现 n 个点的多少种标注？ # ============================================================================= def halfplane_labelings(P): """ 精确枚举平面上 n 个点能被半平面 {x : wᵀx > c} 实现的所有 0/1 标注（labelings / dichotomies）。 【思路】固定方向 w，按投影 wᵀx 从大到小排序，半平面能实现的标注恰好是"排名前 j 的点标 1"（j = 0..n）。 排序只在 w 垂直于某个 x_i − x_j 时才改变（两个点的投影相等），这样的"临界方向"最多 n(n−1) 个。 所以只需要在相邻临界方向之间各取一个方向（取中点），就遍历了所有可能的排序，从而得到全部标注。 这是精确枚举，不是随机搜索。 【参数】P：(n, 2) 点的坐标（不要求一般位置；重合的点会被当作同一个点） 【返回】set，每个元素是一个整数 bitmask：第 i 位是 1 表示第 i 个点标为 1 【Python 小知识】bitmask：sum(1 << i for i in S) 把集合 S 编码成一个整数，方便放进 set 去重。 """ n = len(P) angles = [] for i in range(n): for j in range(i + 1, n): dx, dy = P[j] - P[i] a = np.arctan2(dy, dx) + np.pi / 2 # 垂直于 x_j − x_i 的方向 angles += [a % (2 * np.pi), (a + np.pi) % (2 * np.pi)] angles = np.unique(angles) if angles else np.array([0.0]) mids = (angles + np.roll(angles, -1)) / 2.0 # 相邻临界角的中点 mids[-1] = (angles[-1] + angles[0] + 2 * np.pi) / 2.0 # 最后一个扇区跨过 2π，单独处理 labelings = set() for a in mids: order = np.argsort(-(P @ np.array([np.cos(a), np.sin(a)]))) mask = 0 labelings.add(mask) # j = 0：全标 0（c 取很大） for i in order: # 依次把排名最前的点加进"标 1"的集合 mask |= 1 << int(i) labelings.add(mask) return labelings def describe_missing(labelings, n): """把没能实现的标注写成 '1010' 这样的字符串（第 i 个字符是第 i 个点的标注）。""" missing = sorted(set(range(2 ** n)) - labelings) return [''.join(str((mk >> i) & 1) for i in range(n)) for mk in missing] def demo_vc(): print_header("演示 3：VC 维与 shattering —— 平面上的半平面（线性分类器 + 截距）") rng = np.random.default_rng(53) print("对给定的 n 个点，枚举半平面能实现的全部 0/1 标注。全部 2ⁿ 种都能实现，就说这组点被 shattered。\n") cases = { "3 points, general position": np.array([[0.0, 0.0], [1.0, 0.0], [0.3, 1.0]]), "3 points, collinear": np.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]]), "4 points, square": np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]), "4 points, one inside": np.array([[0.0, 0.0], [2.0, 0.0], [1.0, 2.0], [1.0, 0.6]]), } print(f"{'point set':<30}{'achieved':>10}{'2^n':>6} missing labelings") found = {} for name, P in cases.items(): L = halfplane_labelings(P) found[name] = (len(L), describe_missing(L, len(P))) print(f"{name:<30}{len(L):>10}{2 ** len(P):>6} {', '.join(found[name][1]) if found[name][1] else '(none)'}") print("\n n 个随机点（一般位置）：实现的标注数 vs 2ⁿ vs Cover 公式 vs Sauer 引理的上界（VC 维 d = 3）") print(f"{'n':>4}{'achieved':>10}{'2^n':>8}{'Cover':>8}{'Sauer bound':>13}") counts = {} for n in [3, 4, 5, 6, 8, 10, 15, 20]: L = halfplane_labelings(rng.normal(size=(n, 2))) cover = 2 * sum(comb(n - 1, i) for i in range(3)) # 一般位置的点：恰好这么多（Cover 1965） sauer = sum(comb(n, i) for i in range(4)) # Σ_{i ≤ d} C(n, i)，d = 3 counts[n] = (len(L), cover, sauer) print(f"{n:>4}{len(L):>10}{2 ** n:>8}{cover:>8}{sauer:>13}") print(f""" 观察 Observations： · 一般位置的 3 个点：{found['3 points, general position'][0]}/8 种标注都能实现 —— shattered。所以半平面类的 VC 维 ≥ 3。 · 共线的 3 个点只有 {found['3 points, collinear'][0]}/8（中间的点和两端不同的 {', '.join(found['3 points, collinear'][1])} 做不到）。 这不影响 VC 维 ≥ 3：定义只要求**存在**一组能被 shatter 的点，不要求所有点集都行。 · 4 个点：正方形缺的恰好是两个 XOR 标注（{', '.join(found['4 points, square'][1])}，对角线上的点同类）；一个点在三角形内部时， 缺的是"内部点和外面三个点都不同"的两个标注（{', '.join(found['4 points, one inside'][1])}）。任何 4 个点都有做不到的标注 （Radon 定理），所以 VC 维恰好是 3 = 参数个数（w₁, w₂, c）。一般地，Rᵈ 里的线性分类器 VC 维是 d + 1。 · n 个点上实现的标注数和 Cover 公式完全一致，并且**多项式**增长（n = 20：{counts[20][0]} 种，而 2²⁰ = {2 ** 20}）。 Sauer 引理：VC 维为 d 的类在 n 个点上最多实现 Σ_{{i≤d}} C(n, i) = O(nᵈ) 种标注（上界 {counts[20][2]} ≥ {counts[20][0]}）。 · 这就是 VC 理论把无限假设类变成"有限"的方法：在 n 个点上，无限多条直线只能做出 O(n³) 种不同的事情， union bound 里的 log k 换成 log(O(nᵈ)) ≈ d·log n，于是样本复杂度 ∝ VC 维 d（讲义 §3.3）。""") def main(): demo_hoeffding() demo_uniform_convergence() demo_vc() if __name__ == "__main__": main()