""" ============================================================================= 第 1 讲 · 代码 3：逻辑回归、Newton 法与感知机（logistic regression, Newton's method, perceptron） ============================================================================= 对应讲义：lectures/lecture01_linear_models.md 第 3 节 ----------------------------------------------------------------------------- 【这个文件要回答的问题】 ----------------------------------------------------------------------------- 演示 1：为什么不用线性回归做分类？几个"毫无争议"的极端样本就能把最小二乘的分类边界拉歪。 用 margin 的视角看四种损失：平方损失会惩罚"太正确"的样本。 演示 2：GD vs Newton。Newton 几步就收敛（quadratic convergence），而且不怕特征尺度（affine invariance）。 演示 3：数据线性可分时，逻辑回归的最大似然估计不存在：‖θ‖ 一直增长，loss 一直趋向 0。L2 正则化修好它。 演示 4：感知机。可分时有限步收敛（犯错次数 ≤ (R/γ)²）；不可分时永不停止，最后的 θ 时好时坏。 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python3 03_logistic_newton.py # 约 2 秒 依赖本目录的 datasets.py 和 linear_models.py。 """ import numpy as np from datasets import make_blobs, make_tumor from linear_models import (accuracy, add_intercept, fit_least_squares, fit_logistic_gd, fit_logistic_newton, logistic_loss, perceptron, print_header) # ============================================================================= # 演示 1：线性回归做分类的问题；四种损失 # ============================================================================= def demo_why_not_linear_regression(): print_header("演示 1：为什么不用 least squares 做分类？—— 极端但正确的样本会拉歪边界") rng = np.random.default_rng(31) x0, y0 = make_tumor(30, rng, n_huge=0) X0 = add_intercept(x0) print("良性 ~ N(2, 0.7²)，恶性 ~ N(4, 0.7²)，各 30 个。再往训练集里加 k 个大小在 12~16 cm 的恶性肿瘤（标签毫无争议）。") print("分类规则：least squares 预测 ≥ 0.5 判为恶性；logistic 预测概率 ≥ 0.5 判为恶性。准确率在 60 个原始点上算。\n") print(f"{'k huge tumors':>14}{'LS boundary':>14}{'LS acc':>10}{'logistic boundary':>20}{'logistic acc':>14}") for k in [0, 5, 10, 20, 40]: x_huge = np.random.default_rng(32).uniform(12.0, 16.0, size=k) Xv = add_intercept(np.concatenate([x0, x_huge])) yv = np.concatenate([y0, np.ones(k, dtype=int)]) th_ls = fit_least_squares(Xv, yv) b_ls = (0.5 - th_ls[0]) / th_ls[1] # θ₀ + θ₁ x = 0.5 的解 th_lr, _ = fit_logistic_newton(Xv, yv) b_lr = -th_lr[0] / th_lr[1] # θ₀ + θ₁ x = 0 的解（σ = 0.5） acc_ls = accuracy(y0, (X0 @ th_ls >= 0.5).astype(int)) acc_lr = accuracy(y0, (X0 @ th_lr >= 0).astype(int)) print(f"{k:>14}{b_ls:>12.3f}cm{acc_ls:>10.3f}{b_lr:>18.3f}cm{acc_lr:>14.3f}") print("\n用 margin 看损失。令 y' = 2y − 1 ∈ {−1, +1}，z = θᵀx，margin t = y'·z（t > 0 表示分对，越大越'自信地分对'）：\n") ts = [-2.0, -1.0, 0.0, 1.0, 2.0, 4.0] print(f"{'loss':<34}" + "".join(f"{'t=' + format(t, 'g'):>8}" for t in ts)) rows = [ ("0-1 loss 1{t ≤ 0}", lambda t: float(t <= 0)), ("logistic log(1 + e^{−t})", lambda t: float(np.logaddexp(0.0, -t))), ("hinge max(0, 1 − t)", lambda t: max(0.0, 1.0 - t)), ("squared (1 − t)²", lambda t: (1.0 - t) ** 2), ] for name, f in rows: print(f"{name:<34}" + "".join(f"{f(t):>8.3f}" for t in ts)) print(""" 观察 Observations： · 没有大肿瘤时，两种方法的边界几乎一样（2.89 cm），准确率都是 0.95。 · 加入大肿瘤后，least squares 的边界一路右移（k = 40 时到 3.49 cm），原本分对的恶性样本被判成良性， 准确率降到 0.867；逻辑回归的边界一动不动（小数点后三位都不变）。 · 原因看损失表：对 y = 1 的样本，最小二乘希望 θᵀx 恰好等于 1。大肿瘤的 θᵀx 远大于 1（"太正确"，t = 4）， 平方损失仍然给它 9 的惩罚，于是直线被压平去迁就它们。logistic 损失在 t 很大时趋于 0：分对了就不管了。 · 第三行的 hinge loss 是 SVM 的损失（第 2 讲）：t ≥ 1 后损失恰好为 0。hinge 是 0-1 损失的凸上界 （convex surrogate）；logistic 用 log₂ 时也是（表里是自然对数，差一个常数因子 1/ln 2）。 0-1 损失本身不连续、梯度处处为 0，没法直接优化，所以我们优化它的凸替身。""") # ============================================================================= # 演示 2：GD vs Newton # ============================================================================= def first_below(history, tol): """返回 history 中第一个 < tol 的下标（即"迭代多少步之后梯度小于 tol"）；没有就返回 None。""" for t, v in enumerate(history): if v < tol: return t return None def demo_gd_vs_newton(): print_header("演示 2：Gradient descent vs Newton's method") rng = np.random.default_rng(33) X, y = make_blobs(200, [[0.0, 0.0], [2.0, 1.0]], rng) Xi = add_intercept(X) hist_nt = [] th_nt, it_nt = fit_logistic_newton(Xi, y, history=hist_nt) print("Newton 每一步的梯度大小 ‖∇J‖_∞（注意指数大约每步翻倍 —— quadratic convergence）：") print(" " + " ".join(f"it{t}: {g:.1e}" for t, g in enumerate(hist_nt))) hist_gd = [] th_gd, it_gd = fit_logistic_gd(Xi, y, tol=1e-10, max_iters=200000, history=hist_gd) print("\nGD（步长 1/L）达到各个精度需要的步数：") for tol in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10]: print(f" ‖∇J‖_∞ < {tol:.0e}: GD {first_below(hist_gd, tol)!s:>6} 步 Newton {first_below(hist_nt, tol)!s:>3} 步") print(f"\n两者的解的最大差 = {np.max(np.abs(th_gd - th_nt)):.1e}；θ = {np.round(th_nt, 4)}") print("\n把第一个特征乘以 1000（比如单位从 km 换成 m）—— 同一个问题，不同的坐标：") Xs = Xi * np.array([1.0, 1000.0, 1.0]) hist_gd2, hist_nt2 = [], [] th_gd2, _ = fit_logistic_gd(Xs, y, tol=1e-6, max_iters=200000, history=hist_gd2) th_nt2, it_nt2 = fit_logistic_newton(Xs, y, history=hist_nt2) print(f" Newton：{it_nt2} 步收敛（和原来一样）；θ₁ 变成了原来的 1/1000：{th_nt2[1]:.6f} vs {th_nt[1] / 1000:.6f}") print(f" GD：200000 步后 ‖∇J‖_∞ = {hist_gd2[-1]:.1e}，训练 loss 比最优多 " f"{logistic_loss(th_gd2, Xs, y) - logistic_loss(th_nt2, Xs, y):.2e}") print(""" 观察 Observations： · GD 每一步都把误差乘以一个固定的比例（linear convergence），精度每提高 100 倍要多走差不多相同的步数； Newton 一旦接近最优点，每一步把误差平方（quadratic convergence）：几步就到机器精度。 · Newton 的代价：每步要算 d × d 的 Hessian 并解方程组，O(nd² + d³)。d = 3 时这几乎免费； d = 10⁶（神经网络）时完全不可行 —— 所以深度学习用一阶方法（SGD / Adam），而 GLM 软件（R 的 glm、 statsmodels）用 Newton / IRLS。 · Affine invariance：换坐标（缩放特征）后，Newton 的迭代路径"同步变换"，步数完全不变； GD 却被条件数拖垮（和代码 01 的演示 3 一样）。Newton 等于自动做了"最好的" feature scaling。""") # ============================================================================= # 演示 3：线性可分时 MLE 不存在 # ============================================================================= def demo_separable(): print_header("演示 3：线性可分的数据 —— MLE 不存在，‖θ‖ → ∞") rng = np.random.default_rng(34) X, y = make_blobs(50, [[0.0, 0.0], [5.0, 5.0]], rng) Xi = add_intercept(X) th_p, mistakes = perceptron(Xi, y, max_epochs=1000) print(f"先确认数据线性可分：感知机在第 {len(mistakes)} 个 epoch 犯错 0 次 → 存在一个把两类完全分开的超平面。\n") print(f"{'GD 步数':>10}{'‖θ‖ (λ=0)':>14}{'loss (λ=0)':>14}{'方向变化 (度)':>16}{'‖θ‖ (λ=0.01)':>16}") prev_dir = None for T in [10, 100, 1000, 10000, 100000]: th, _ = fit_logistic_gd(Xi, y, tol=0.0, max_iters=T) # tol = 0：强制跑满 T 步 th_r, _ = fit_logistic_gd(Xi, y, l2=0.01, tol=0.0, max_iters=T) d = th / np.linalg.norm(th) angle = "" if prev_dir is None else f"{np.degrees(np.arccos(np.clip(d @ prev_dir, -1, 1))):.3f}" prev_dir = d print(f"{T:>10}{np.linalg.norm(th):>14.3f}{logistic_loss(th, Xi, y):>14.2e}{angle:>16}" f"{np.linalg.norm(th_r):>16.4f}") print(""" 观察 Observations： · λ = 0：步数每多 10 倍，‖θ‖ 只增加 3~4（对数增长），loss 缩小约 7 倍 —— 永远在下降，永远到不了底。 原因：所有样本都分对了，把 θ 放大 c 倍会让每个 margin 都放大 c 倍、loss 更小。最优点在无穷远处。 · 但 θ 的**方向**稳定下来了（方向变化越来越小）。Soudry et al.（2018）证明：GD 在可分数据上的方向 收敛到 max-margin（硬间隔 SVM）的方向 —— 这就是第 2 讲 SVM 和第 3 讲 implicit regularization 的连接点。 · λ = 0.01：‖θ‖ 很快稳定在一个有限值，问题有唯一的解。实践中（sklearn 的 LogisticRegression 默认 C = 1） 总是带正则化，部分原因就在这里。 · 预测出的"概率"在 λ = 0 时会越来越接近 0 / 1 —— 模型过度自信（overconfident），概率失去意义。""") # ============================================================================= # 演示 4：感知机 # ============================================================================= def margin(theta, X, y): """几何间隔 γ = min_i y'ᵢ θᵀxᵢ / ‖θ‖（y' = 2y − 1）。γ > 0 表示 θ 把数据完全分开。""" return float(np.min((2 * y - 1) * (X @ theta)) / np.linalg.norm(theta)) def demo_perceptron(): print_header("演示 4：Perceptron —— 可分时有限步收敛，不可分时永不停止") rng = np.random.default_rng(35) X, y = make_blobs(50, [[0.0, 0.0], [4.0, 4.0]], rng) Xi = add_intercept(X) th, mistakes = perceptron(Xi, y, max_epochs=1000) R = float(np.max(np.linalg.norm(Xi, axis=1))) th_ref, _ = fit_logistic_newton(Xi, y, l2=1e-3) # 一个"参考"分隔方向，用来算 Novikoff 上界 g_ref = margin(th_ref, Xi, y) print(f"(a) 可分数据（中心距离 4√2，n = 100）：每个 epoch 的犯错次数 = {mistakes}") print(f" 总犯错次数 = {sum(mistakes)}；R = max‖x‖ = {R:.3f}") if g_ref > 0: print(f" 用正则化逻辑回归的方向当作参考分隔面 u：γ_u = {g_ref:.4f}，Novikoff 上界 (R/γ_u)² = {(R / g_ref) ** 2:.0f}") print(f" 感知机自己最后的 θ 的间隔 = {margin(th, Xi, y):.4f}（> 0，分开了，但比参考方向的间隔小）") rng2 = np.random.default_rng(36) X2, y2 = make_blobs(100, [[0.0, 0.0], [2.0, 1.0]], rng2) X2i = add_intercept(X2) th_lr, _ = fit_logistic_newton(X2i, y2) print(f"\n(b) 不可分数据（和演示 2 同样的分布）：逻辑回归的训练准确率 = {accuracy(y2, (X2i @ th_lr >= 0).astype(int)):.3f}") print(" 感知机每个 epoch 结束时 θ 的训练准确率（以及这个 epoch 的犯错次数）：") accs = [] for e in range(12): th_e, m = perceptron(X2i, y2, max_epochs=e + 1) # 每次从头跑 e+1 个 epoch（顺序固定，结果确定） accs.append((accuracy(y2, (X2i @ th_e >= 0).astype(int)), m[-1])) print(" " + " ".join(f"e{e + 1}: {a:.3f}({m})" for e, (a, m) in enumerate(accs))) print(""" 观察 Observations： · 可分时感知机在有限步内停下，实际犯错次数远小于 Novikoff 上界 (R/γ)²（上界是最坏情况）。 上界只和"数据有多容易分"（间隔 γ 相对于数据半径 R）有关，和样本数、维度无关 —— 这是最早的 "margin 决定泛化 / 收敛"的结果之一，SVM（第 2 讲）把这个思想推到了极致。 · 感知机找到的是"某一个"分隔面，不是间隔最大的那个；它停下来的地方取决于样本顺序。 · 不可分时它每个 epoch 都在犯错、都在更新，θ 在不同的解之间跳，最后一个 epoch 的准确率时高时低。 逻辑回归给出唯一的、带概率的答案。（averaged perceptron —— 对所有 epoch 的 θ 取平均 —— 是一个常用的修补。） · 历史：Rosenblatt 1958 年的感知机让人们对神经网络非常乐观；Minsky & Papert 1969 年指出单层感知机 连 XOR 都学不了（不是线性可分的），神经网络研究随之沉寂了十几年。解决办法是多层 + 可导的激活函数 + backprop（第 4 讲）。""") def main(): demo_why_not_linear_regression() demo_gd_vs_newton() demo_separable() demo_perceptron() if __name__ == "__main__": main()