"""
=============================================================================
第 2 讲 · 代码 4：支持向量机与 SMO（SVMs: max margin, soft margin, kernels, implicit bias of GD）
=============================================================================

对应讲义：lectures/lecture02_generative_kernels_svm.md 第 5 节

-----------------------------------------------------------------------------
【这个文件要回答的问题】
-----------------------------------------------------------------------------
    演示 1：Hard-margin SVM。functional vs geometric margin；support vectors 恰好在间隔边界上；
            只用 support vectors 重新训练得到同一个解；和感知机、逻辑回归（第 1 讲）的间隔比较。
    演示 2：Soft margin 与 C。C 从 0.001 到 100：support vector 个数、间隔宽度、训练 / 测试准确率；
            primal（hinge loss + L2）与 dual 的目标值相等（duality gap ≈ 0）；和逻辑回归比较。
    演示 3：Kernel SVM。同心圆数据上线性 SVM 失败，二次多项式核和 RBF 核成功；RBF 带宽 σ 的影响。
    演示 4：第 1 讲 §3.5 的承诺 —— 可分数据上，逻辑回归的 GD 方向收敛到 max-margin 方向（Soudry et al. 2018）。
            截距被当作 x 的一个分量时，对应的是一个**截距也被惩罚**的 SVM；什么时候它和标准 SVM 不一样？

-----------------------------------------------------------------------------
【运行方式】
-----------------------------------------------------------------------------
    python3 04_svm_smo.py          # 约 10 秒（演示 4 跑两次 10⁶ 步 GD，占大部分时间）

依赖本目录的 kernel_methods.py，以及第 1 讲的 linear_models.py / datasets.py（kernel_methods 导入时自动加入搜索路径，
所以下面要先 import kernel_methods）。
"""

import numpy as np

from kernel_methods import (dcd_svm_nobias, linear_kernel, make_circles, poly_kernel, rbf_kernel, smo_train,
                            svm_decision, svm_dual_objective, svm_primal_objective, to_pm1)
from datasets import make_blobs
from linear_models import (add_intercept, fit_glm_newton, fit_logistic_newton, logistic_grad, perceptron,
                           print_header)


def geometric_margin(w, b, X, y):
    """
    几何间隔 γ = minᵢ yᵢ (wᵀx⁽ⁱ⁾ + b) / ‖w‖（y ∈ {−1, +1}）：最近的点到超平面的（带符号的）欧氏距离。

    【和 functional margin 的区别】functional margin γ̂ᵢ = yᵢ(wᵀx⁽ⁱ⁾ + b) 可以靠把 (w, b) 乘以 10 变大 10 倍，
        超平面却没变；除以 ‖w‖ 之后才是真正的几何距离，和缩放无关。γ < 0 表示有点被分错了。
    """
    return float(np.min(y * (X @ w + b)) / np.linalg.norm(w))


def angle_deg(a, b):
    """两个向量的夹角（度）。用 clip 防止舍入让 cos 略大于 1 导致 arccos 返回 nan。"""
    c = a @ b / (np.linalg.norm(a) * np.linalg.norm(b))
    return float(np.degrees(np.arccos(np.clip(c, -1.0, 1.0))))


# =============================================================================
# 演示 1：Hard-margin SVM
# =============================================================================
def demo_hard_margin():
    print_header("演示 1：Hard-margin SVM —— 间隔、support vectors、稀疏性")
    rng = np.random.default_rng(35)                     # 和第 1 讲代码 03 演示 4(a)（感知机）同一份可分数据
    X, y01 = make_blobs(50, [[0.0, 0.0], [4.0, 4.0]], rng)
    y = to_pm1(y01)
    K = linear_kernel(X, X)
    model = smo_train(K, y, C=np.inf, tol=1e-10)
    alpha = model["alpha"]
    w, b = X.T @ (alpha * y), model["b"]
    sv = alpha > 1e-8
    fm = y * (X @ w + b)                                # 每个点的 functional margin
    print(f"第 1 讲感知机演示用的可分数据（两团，中心距离 4√2，n = 100）。SMO（C = ∞）更新了 {model['n_iter']} 对 α。")
    print(f"w = [{w[0]:.4f}, {w[1]:.4f}]，b = {b:.4f}，几何间隔 γ = 1/‖w‖ = {1 / np.linalg.norm(w):.4f}")
    print(f"support vectors（α > 0）：{int(sv.sum())} 个，下标 {np.flatnonzero(sv).tolist()}，α = "
          f"{np.round(alpha[sv], 4).tolist()}")
    print(f"它们的 functional margin yᵢ(wᵀxᵢ + b) = {np.round(fm[sv], 6).tolist()}；"
          f"其余点的最小 functional margin = {fm[~sv].min():.4f}")
    print(f"KKT / 对偶可行性：Σ αᵢyᵢ = {np.sum(alpha * y):.1e}；primal ½‖w‖² = {0.5 * w @ w:.6f}，"
          f"dual W(α) = {svm_dual_objective(K, y, alpha):.6f}")
    print(f"把 (w, b) 乘以 10：functional margin 的最小值 {fm.min():.3f} → {(10 * fm).min():.3f}，"
          f"几何间隔 {geometric_margin(w, b, X, y):.4f} → {geometric_margin(10 * w, 10 * b, X, y):.4f}")

    # 只用 support vectors 重新训练
    Xs, ys = X[sv], y[sv]
    m2 = smo_train(linear_kernel(Xs, Xs), ys, C=np.inf, tol=1e-10)
    w2, b2 = Xs.T @ (m2["alpha"] * ys), m2["b"]
    print(f"只用这 {int(sv.sum())} 个 support vectors 重新训练：(w, b) 的最大差 = "
          f"{max(np.max(np.abs(w2 - w)), abs(b2 - b)):.1e}")

    # 和第 1 讲的感知机、逻辑回归比较间隔。两种量法：
    #   x 空间的几何间隔 min yᵢ(wᵀxᵢ + b)/‖w‖（SVM 最大化的就是它）；
    #   增广空间 [1, x] 里的间隔 min yᵢ θᵀx̃ᵢ/‖θ‖，θ = [b, w] —— 第 1 讲代码 03 演示 4 的 γ 就是这样算的
    #   （感知机把截距当作特征，Novikoff 界里的 R、γ 都在增广空间里）。
    Xi = add_intercept(X)
    th_p, _ = perceptron(Xi, y01, max_epochs=1000)
    th_l, _ = fit_logistic_newton(Xi, y01, l2=1e-3)
    th_nb = dcd_svm_nobias(Xi, y, C=np.inf, tol=1e-12)["w"]      # 增广空间里间隔最大的方向（演示 4 再细说）
    print(f"\n{'separator':<38}{'margin in x space':>19}{'margin in [1,x] space':>23}")
    rows = [("perceptron (Lecture 1)", th_p),
            ("logistic, lambda = 1e-3 (Lecture 1)", th_l),
            ("hard-margin SVM (SMO, b free)", np.concatenate([[b], w])),
            ("max margin in [1,x] space (DCD)", th_nb)]
    for name, th in rows:
        aug = float(np.min(y * (Xi @ th)) / np.linalg.norm(th))
        print(f"{name:<38}{geometric_margin(th[1:], th[0], X, y):>19.4f}{aug:>23.4f}")
    R = float(np.max(np.linalg.norm(Xi, axis=1)))
    g_aug = 1.0 / np.linalg.norm(th_nb)
    print(f"用增广空间里最大的间隔 γ* = {g_aug:.4f}，Novikoff 上界 (R/γ*)² = ({R:.3f}/{g_aug:.4f})² = {(R / g_aug) ** 2:.0f}"
          f"（第 1 讲用逻辑回归的方向得到 759；感知机实际犯错 10 次）")
    print(f"""
观察 Observations：
  · 100 个点里只有 {int(sv.sum())} 个 support vector，它们的 functional margin 恰好是 1（KKT 的互补松弛条件：
    αᵢ > 0 ⇒ 约束取等号）；其余点都在间隔外面（最小 {fm[~sv].min():.3f} > 1），α = 0，对解没有任何贡献。
  · 稀疏性的直接证明：扔掉其他 {100 - int(sv.sum())} 个点，只用 support vectors 重新训练，得到同一个超平面。
    这也是 SVM 预测快的原因：f(x) = Σᵢ αᵢyᵢ K(xᵢ, x) + b 只需要对 support vectors 求和。
  · functional margin 可以靠缩放任意放大，几何间隔不变 —— 所以 SVM 固定“最小 functional margin = 1”，
    再最小化 ‖w‖，等价于最大化几何间隔 1/‖w‖。
  · 四个分隔面都把数据完全分开，但感知机的间隔最小（它停在“某一个”分得开的地方），
    正则化的逻辑回归接近、但不等于 max margin（λ → 0 时会收敛过去，见演示 4）。
    “最大”要说清楚在哪个空间：b 不受惩罚的 SVM 在 x 空间里间隔最大；把截距当作特征时，
    增广空间 [1, x] 里间隔最大的是另一个问题（截距也进入 ‖·‖）。这份数据上两者恰好相同（演示 4 解释为什么）。
  · 第 1 讲的 Novikoff 上界 759 用的是逻辑回归方向的间隔；换成最大间隔 γ*，上界缩小到 {(R / g_aug) ** 2:.0f}，
    仍远大于实际的 10 次 —— 上界是最坏情况。
  · primal 和 dual 的目标值相等（强对偶）：这是检查求解器是否收敛到最优的好办法。""")


# =============================================================================
# 演示 2：Soft margin 与 C
# =============================================================================
def demo_soft_margin():
    print_header("演示 2：Soft margin —— C 的作用；hinge loss + L2；和逻辑回归比较")
    rng = np.random.default_rng(62)
    X, y01 = make_blobs(100, [[0.0, 0.0], [2.0, 1.0]], rng)
    Xte, yte01 = make_blobs(2500, [[0.0, 0.0], [2.0, 1.0]], rng)
    mean = X.mean(axis=0)
    X, Xte = X - mean, Xte - mean                        # 中心化（原因见下面的说明）
    y, yte = to_pm1(y01), to_pm1(yte01)
    n = len(y)
    K, Kte = linear_kernel(X, X), linear_kernel(Xte, X)
    print("两团重叠的高斯（和第 1 讲 GD vs Newton 同一种分布），训练 200 个、测试 5000 个；特征先中心化。")
    print("soft margin：min ½‖w‖² + C Σᵢ max(0, 1 − yᵢ(wᵀxᵢ + b)) ⇔ min (1/n)Σ hinge + (λ/2)‖w‖²，λ = 1/(nC)。")
    print("最后一列是同样 λ 的 L2 逻辑回归（第 1 讲的实现连截距一起惩罚；特征中心化、两类各半后，最优截距本来就接近 0，")
    print("所以这个差别影响很小）。\n")
    print(f"{'C':>8}{'#pairs':>9}{'#SV':>6}{'#(a=C)':>8}{'2/|w|':>8}{'train acc':>11}{'test acc':>10}"
          f"{'rel. gap':>10}{'LR test acc':>13}")
    out = {}
    for C in [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]:
        m = smo_train(K, y, C, tol=1e-6, max_iter=1000000)
        a = m["alpha"]
        w = X.T @ (a * y)
        p_obj = svm_primal_objective(K, y, a, m["b"], C)
        gap = (p_obj - svm_dual_objective(K, y, a)) / p_obj
        tr = np.mean(np.sign(svm_decision(K, y, m)) == y)
        te = np.mean(np.sign(svm_decision(Kte, y, m)) == yte)
        th, _ = fit_logistic_newton(add_intercept(X), y01, l2=1.0 / (n * C))
        te_lr = np.mean((add_intercept(Xte) @ th >= 0) == (yte > 0))
        nsv, nb = int(np.sum(a > 1e-8)), int(np.sum(a >= C * (1 - 1e-9)))
        out[C] = (nsv, nb, 2 / np.linalg.norm(w), tr, te, te_lr, m["n_iter"])
        print(f"{C:>8g}{m['n_iter']:>9}{nsv:>6}{nb:>8}{2 / np.linalg.norm(w):>8.3f}{tr:>11.3f}{te:>10.4f}"
              f"{gap:>10.1e}{te_lr:>13.4f}")
    Cs = sorted(out)
    max_diff = max(abs(out[c][4] - out[c][5]) for c in Cs)
    # 对照：不中心化（原始坐标）时，C = 0.001 对应的强 L2 会把截距也压向 0
    th_raw, _ = fit_logistic_newton(add_intercept(X + mean), y01, l2=1.0 / (n * 0.001))
    te_raw = np.mean((add_intercept(Xte + mean) @ th_raw >= 0) == (yte > 0))
    print(f"对照：**不**中心化时，C = 0.001 对应的逻辑回归（λ = {1 / (n * 0.001):g}）测试准确率 = {te_raw:.4f}")
    print(f"""
观察 Observations：
  · C 小（强正则化）：间隔很宽（C = 0.001 时 2/‖w‖ = {out[0.001][2]:.2f}），几乎所有点都在间隔内、都是 support vector
    （{out[0.001][0]} / 200，而且全部 α = C）。C 大（弱正则化）：间隔收窄到 {out[100.0][2]:.2f} 左右后不再变，
    support vector 剩 {out[100.0][0]} 个 —— 数据不可分，hinge loss 的最小值本身是有限的，C 再大也只是在逼近它。
  · 所有 C 的测试准确率都在 {min(out[c][4] for c in Cs):.4f}～{max(out[c][4] for c in Cs):.4f} 之间，最强的正则化（C = 0.001）反而略好：
    这份二维数据很简单、Bayes 边界就是一条直线，正则化强弱影响不大。C 真正重要的场合是高维、噪声大的数据
    （第 3 讲用交叉验证选它）。
  · rel. gap = (primal − dual)/primal ≈ 0：SMO 的解满足强对偶。#pairs 是 SMO 更新了多少对 α：C 越大越慢
    （C = 100 要 {out[100.0][6]} 对，C = 1 只要 {out[1.0][6]} 对）—— 很多 α 要一小步一小步地爬到很大的值。
    实践中常把特征标准化、用不太大的 C。
  · 同样 λ 的逻辑回归，测试准确率和 SVM 最多差 {max_diff:.4f}。两者只是损失函数不同：hinge 在 margin ≥ 1 后恰好为 0
    （稀疏的 α），logistic 永远不为 0（每个点都有一点贡献，但给出概率）。
    不中心化时，C = 0.001 对应的 λ = {1 / (n * 0.001):g} 会把第 1 讲实现里的截距也压向 0，逻辑回归掉到 {te_raw:.4f}
    —— SVM 的 b 不受惩罚。比较两个模型时，要先确认它们解的是“同一个问题”。""")


# =============================================================================
# 演示 3：Kernel SVM
# =============================================================================
def demo_kernel_svm():
    print_header("演示 3：Kernel SVM —— 同心圆：线性模型做不到，换个核就行")
    rng = np.random.default_rng(63)
    X, y01 = make_circles(100, rng)
    Xte, yte01 = make_circles(1000, rng)
    y, yte = to_pm1(y01), to_pm1(yte01)
    print("内圈（半径 1）是负类，外圈（半径 2）是正类，半径噪声 0.25；训练 200 个，测试 2000 个。C = 1。\n")
    kernels = [("linear", lambda A, B: linear_kernel(A, B)),
               ("poly, degree 2", lambda A, B: poly_kernel(A, B, 2, 1.0))]
    for s in [0.03, 0.1, 0.3, 1.0, 3.0, 10.0]:
        kernels.append((f"RBF, sigma = {s:g}", (lambda s_: lambda A, B: rbf_kernel(A, B, s_))(s)))
    print(f"{'kernel':<18}{'#pairs':>8}{'#SV':>6}{'train acc':>11}{'test acc':>10}")
    res = {}
    for name, kf in kernels:
        K = kf(X, X)
        m = smo_train(K, y, 1.0, tol=1e-6)
        tr = np.mean(np.sign(svm_decision(K, y, m)) == y)
        te = np.mean(np.sign(svm_decision(kf(Xte, X), y, m)) == yte)
        res[name] = (int(np.sum(m["alpha"] > 1e-8)), tr, te)
        print(f"{name:<18}{m['n_iter']:>8}{res[name][0]:>6}{tr:>11.3f}{te:>10.4f}")
    print(f"""
观察 Observations：
  · 线性 SVM 只有 {res['linear'][2]:.3f} 的测试准确率（接近瞎猜）：没有任何直线能切开两个同心圆。
  · 二次多项式核的特征里有 x₁²、x₂²，在这些特征上“半径的平方 = 常数”是一个超平面，
    于是测试准确率 {res['poly, degree 2'][2]:.3f}，而且只用了 {res['poly, degree 2'][0]} 个 support vectors。
    模型 + 优化器一个字没改，只换了 K —— 这就是 kernel trick 的威力：算法只通过内积接触数据。
  · RBF 的 σ 又是复杂度旋钮：σ = 0.03 时训练准确率 {res['RBF, sigma = 0.03'][1]:.3f}、测试只有 {res['RBF, sigma = 0.03'][2]:.3f}，
    每个点都是 support vector（{res['RBF, sigma = 0.03'][0]} / 200）—— 每个点只“照亮”自己周围一小块，等于在背训练集；
    σ = 0.3～3 时测试准确率约 0.97；σ = 10 时所有点之间的核值都接近 1（“一样相似”），在 C = 1 下模型几乎弯不动：
    欠拟合，训练准确率只有 {res['RBF, sigma = 10'][1]:.3f}，所有 {res['RBF, sigma = 10'][0]} 个点都落在间隔里成了 support vector。
  · support vector 的个数是一个粗略的复杂度指标：太多（接近 n）通常意味着过拟合或者 C / σ 没选好。
    （留一法误差的经典上界：LOO error ≤ #SV / n。）""")


# =============================================================================
# 演示 4：GD 的隐式偏置 → max margin
# =============================================================================
CHECKPOINTS = [10, 100, 1000, 10000, 100000, 1000000]


def gd_direction_trace(Xi, y01, refs):
    """
    在 (Xi, y01) 上跑逻辑回归的 GD（和第 1 讲 fit_logistic_gd 一样：θ = 0 出发，步长 1/L，λ = 0，不提前停止），
    在 CHECKPOINTS 处记录 ‖θ‖ 和 θ 与每个参考方向的夹角。

    【为什么不直接调用 fit_logistic_gd】它每次从头跑；这里要在一条轨迹上的多个时刻“拍照”，自己写循环更省时间。
        更新公式逐位相同（都用 linear_models.logistic_grad），所以第 1 讲表格里的 ‖θ‖ 可以直接对上。
    【返回】列表，每个元素是 (t, ‖θ_t‖, [夹角1, 夹角2, …])
    """
    n = len(y01)
    lr = 1.0 / (np.linalg.eigvalsh(Xi.T @ Xi / n)[-1] / 4.0)
    theta = np.zeros(Xi.shape[1])
    trace = []
    for t in range(1, CHECKPOINTS[-1] + 1):
        theta = theta - lr * logistic_grad(theta, Xi, y01)
        if t in CHECKPOINTS:
            trace.append((t, float(np.linalg.norm(theta)), [angle_deg(theta, r) for r in refs]))
    return trace


def svm_pair(X, y):
    """
    同一份数据上的两个 hard-margin SVM，都写成 [b, w₁, w₂, …] 的形式（和 add_intercept 的约定一致）：
        标准 SVM：     min ½‖w‖²         s.t. yᵢ(wᵀxᵢ + b) ≥ 1   —— b 不受惩罚（SMO 求解）
        截距被惩罚的：  min ½(b² + ‖w‖²)  s.t. 同上              —— 等价于把 1 当作 x 的一个分量（DCD 求解）
    【返回】(std, nobias)，两个 (d+1,) 数组
    """
    m = smo_train(linear_kernel(X, X), y, C=np.inf, tol=1e-12)
    std = np.concatenate([[m["b"]], X.T @ (m["alpha"] * y)])
    nobias = dcd_svm_nobias(add_intercept(X), y, C=np.inf, tol=1e-12, max_epochs=100000)["w"]
    return std, nobias


def demo_implicit_bias():
    print_header("演示 4：可分数据上 GD 的方向 → max-margin 方向（第 1 讲 §3.5 的承诺）")
    # (a) 第 1 讲代码 03 演示 3 的数据
    rng = np.random.default_rng(34)
    X, y01 = make_blobs(50, [[0.0, 0.0], [5.0, 5.0]], rng)
    y = to_pm1(y01)
    Xi = add_intercept(X)
    std, nob = svm_pair(X, y)
    print("(a) 第 1 讲代码 03 演示 3 的可分数据（截距是 x 的第 0 个分量，GD 对 [θ₀, θ₁, θ₂] 一视同仁）。")
    print(f"    标准 SVM（b 不惩罚）[b, w] = {np.round(std, 4).tolist()}")
    print(f"    截距也被惩罚的 SVM [b, w] = {np.round(nob, 4).tolist()}；两者夹角 = {angle_deg(std, nob):.4f}°")
    n_sv = int(np.sum(np.isclose(y * (Xi @ nob), 1.0, atol=1e-6)))
    print(f"    （间隔恰好为 1 的点有 {n_sv} 个：3 个等式正好确定 3 个参数 [b, w₁, w₂]，所以两种 SVM 在这里给出同一个解。）\n")
    trace_a = gd_direction_trace(Xi, y01, [nob])
    print(f"{'GD steps':>10}{'|theta|':>10}{'angle to SVM (deg)':>20}")
    for t, nrm, angs in trace_a:
        print(f"{t:>10}{nrm:>10.3f}{angs[0]:>20.4f}")
    print(f"\n    对照：L2 正则化逻辑回归（Newton + line search）当 λ → 0 时的方向：")
    print(f"{'lambda':>10}{'|theta|':>10}{'angle to SVM (deg)':>20}")
    path_a = []
    for lam in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]:
        th, _ = fit_glm_newton(Xi, y01, "bernoulli", l2=lam, max_iters=200)
        path_a.append((lam, angle_deg(th, nob)))
        print(f"{lam:>10g}{np.linalg.norm(th):>10.3f}{path_a[-1][1]:>20.4f}")

    # (b) 一个手工构造的小数据集：标准 SVM 只有 2 个 support vectors，截距没被“钉死”
    Xb = np.array([[1.0, 2.0], [0.0, 4.0], [0.0, 0.0], [3.0, 2.0], [4.0, 4.0], [4.0, 0.0]])
    yb01 = np.array([0, 0, 0, 1, 1, 1])
    yb = to_pm1(yb01)
    Xbi = add_intercept(Xb)
    std_b, nob_b = svm_pair(Xb, yb)
    print("\n(b) 手工构造的 6 个点：负类 (1,2) (0,4) (0,0)，正类 (3,2) (4,4) (4,0)。")
    print(f"    标准 SVM [b, w] = {np.round(std_b, 4).tolist()}：竖线 x₁ = 2，几何间隔 "
          f"{geometric_margin(std_b[1:], std_b[0], Xb, yb):.4f}，只有 2 个 support vectors")
    print(f"    截距也被惩罚的 SVM [b, w] = {np.round(nob_b, 4).tolist()}：几何间隔 "
          f"{geometric_margin(nob_b[1:], nob_b[0], Xb, yb):.4f}；两者夹角 = {angle_deg(std_b, nob_b):.2f}°\n")
    trace_b = gd_direction_trace(Xbi, yb01, [nob_b, std_b])
    print(f"{'GD steps':>10}{'|theta|':>10}{'angle to penalized-b SVM':>26}{'angle to standard SVM':>23}")
    for t, nrm, angs in trace_b:
        print(f"{t:>10}{nrm:>10.3f}{angs[0]:>26.4f}{angs[1]:>23.4f}")
    th_b, _ = fit_glm_newton(Xbi, yb01, "bernoulli", l2=1e-8, max_iters=200)
    print(f"    对照：λ = 1e-8 的 L2 逻辑回归：与截距被惩罚的 SVM 夹角 {angle_deg(th_b, nob_b):.4f}°，"
          f"与标准 SVM 夹角 {angle_deg(th_b, std_b):.4f}°")
    a_end, b_end = trace_b[-1][2]
    print(f"""
观察 Observations：
  · (a) GD 的方向确实收敛到 max-margin 方向：夹角从 {trace_a[0][2][0]:.1f}°（10 步）降到 {trace_a[-1][2][0]:.3f}°（10⁶ 步），
    而 ‖θ‖ 像 log t 一样慢慢长（第 1 讲 §3.5 的表格：10⁵ 步时 11.897，这里逐位相同）。
    λ → 0 的正则化路径也收敛到同一个方向（λ = 1e−6 时 {path_a[-1][1]:.3f}°，Rosset, Zhu & Hastie, NIPS 2003）：
    “不加正则化跑很久的 GD”和“加一点点 L2”殊途同归 —— 这就是 implicit regularization（第 3 讲）。
  · 注意第 1 讲把截距当作 x 的一个分量，GD 对 θ₀ 和其他分量一视同仁，所以它对应的 SVM 是
    min ½(b² + ‖w‖²)，**截距也被惩罚**。在 (a) 里这无所谓：3 个 support vectors 恰好钉死 3 个参数，两种 SVM 相同。
  · (b) 这个差别就显现了：标准 SVM 是竖线 x₁ = 2（间隔 1，只有 2 个 support vectors，b = −2）；
    截距被惩罚的 SVM 宁可把间隔降到 {geometric_margin(nob_b[1:], nob_b[0], Xb, yb):.3f}，换一个小一点的 |b|，两者差 {angle_deg(std_b, nob_b):.1f}°。
    GD 走向的是后者：10⁶ 步后和它差 {a_end:.2f}°，和标准 SVM 差 {b_end:.1f}°。
    所以 GD 最大化的是**增广空间 [1, x] 里**的间隔，不是 x 空间里的几何间隔 —— 它依赖于坐标原点放在哪里。
  · 收敛非常慢：(b) 里从 10⁴ 步到 10⁶ 步，夹角只从 {trace_b[3][2][0]:.2f}° 降到 {a_end:.2f}°（缩小 {trace_b[3][2][0] / a_end:.2f} 倍）。
    Soudry et al. (2018) 证明方向误差一般按 O(1/log t) 下降，1/log t 在这段的预测是缩小 log(10⁶)/log(10⁴) = 1.5 倍，
    和实际吻合。(a) 在 10³～10⁶ 步之间降得比 1/log t 快得多（每 10 倍步数缩小 2～5 倍）：那里的几何比较“友好”，
    O(1/log t) 是一般情况下的速度，不是每个数据集都这么慢。别指望靠“多跑几步”得到 SVM。""")


def main():
    demo_hard_margin()
    demo_soft_margin()
    demo_kernel_svm()
    demo_implicit_bias()


if __name__ == "__main__":
    main()