"""
=============================================================================
第 1 讲 · 代码 1：线性回归 —— 三条路到同一个答案，以及每条路的坑
=============================================================================

对应讲义：lectures/lecture01_linear_models.md 第 2.1–2.4 节

-----------------------------------------------------------------------------
【这个文件要回答的问题 What this file answers】
-----------------------------------------------------------------------------
你大概已经知道线性回归可以用 normal equation 一步解出来，也可以用 gradient descent 迭代。
这里不重复"怎么做"，而是用实验回答几个"为什么 / 什么时候会出问题"：

    演示 1：Normal equation、batch GD、SGD 在房价数据上给出同一个 θ，并且接近生成数据的真值。
    演示 2：学习率的上限是 2/λ_max —— 越过它一点点，GD 就从收敛变成发散。
    演示 3：Feature scaling 为什么重要：条件数 κ = λ_max/λ_min 决定 GD 要走多少步。
    演示 4：常数学习率的 SGD 停在最优点附近"抖动"；学习率衰减才能真正收敛。
    演示 5：概率解释 —— 最小二乘 = 高斯噪声下的最大似然；所以它也继承了高斯的弱点：怕离群点（outliers）。
            换成 Laplace 噪声（L1 损失）就稳健得多。
    演示 6：XᵀX 不可逆时（重复特征、d > n）normal equation 失效；lstsq 给出最小范数解。

-----------------------------------------------------------------------------
【运行方式】
-----------------------------------------------------------------------------
    python3 01_linear_regression.py          # 约 1 秒

需要 numpy。依赖本目录的 datasets.py 和 linear_models.py。
"""

import numpy as np

from datasets import HOUSING_NOISE_STD, HOUSING_TRUE_THETA, make_housing
from linear_models import (Standardizer, add_intercept, fit_gd_linear, fit_least_squares, fit_normal_equation,
                           fit_sgd_linear, linreg_cost, print_header)


def fmt(v, spec=".4g"):
    """把向量格式化成一行字符串，例如 [50.1, 0.12, 9.87]。"""
    return "[" + ", ".join(format(float(x), spec) for x in v) + "]"


def to_original_units(theta_std, sc):
    """
    把"标准化特征上的 θ"换算回"原始单位上的 θ"。

    【思路】标准化后 x̃_j = (x_j − m_j)/s_j。模型 θ̃₀ + Σ θ̃_j x̃_j = (θ̃₀ − Σ θ̃_j m_j/s_j) + Σ (θ̃_j/s_j) x_j，
        所以原始单位下：θ_j = θ̃_j / s_j，截距 θ₀ = θ̃₀ − Σ θ̃_j m_j / s_j。
        这让我们可以把学到的系数和生成数据用的真值（以 k$/ft² 等为单位）直接比较。
    """
    slopes = theta_std[1:] / sc.std
    intercept = theta_std[0] - np.sum(theta_std[1:] * sc.mean / sc.std)
    return np.concatenate([[intercept], slopes])


# =============================================================================
# 演示 1：三条路，同一个答案
# =============================================================================
def demo_three_solvers():
    print_header("演示 1：Normal equation vs batch GD vs SGD —— same θ, different cost")
    rng = np.random.default_rng(1)
    X, y = make_housing(200, rng)
    sc = Standardizer().fit(X)
    Xs = add_intercept(sc.transform(X))

    th_ne = fit_normal_equation(Xs, y)
    th_gd, it_gd, status = fit_gd_linear(Xs, y, lr=0.5, tol=1e-8)
    th_sgd = fit_sgd_linear(Xs, y, lr=0.02, epochs=50, rng=np.random.default_rng(11), decay=True)

    print(f"数据：n = 200 套房子，特征 = [面积 ft², 卧室数]，价格单位 k$，噪声标准差 = {HOUSING_NOISE_STD:g}")
    print(f"真值 true θ (原始单位)       = {fmt(HOUSING_TRUE_THETA)}   ← [截距, 每 ft² 的价格, 每间卧室的价格]\n")
    print(f"{'solver':<22}{'theta (standardized)':<30}{'theta (original units)':<28}{'train J':>10}")
    for name, th in [("normal equation", th_ne), (f"batch GD ({it_gd} it)", th_gd), ("SGD (50 epochs)", th_sgd)]:
        print(f"{name:<22}{fmt(th):<30}{fmt(to_original_units(th, sc)):<28}{linreg_cost(th, Xs, y):>10.2f}")
    print(f"\n  GD 状态：{status}；GD 与 normal equation 的最大差 = {np.max(np.abs(th_gd - th_ne)):.1e}")
    print(f"  SGD 与 normal equation 的最大差 = {np.max(np.abs(th_sgd - th_ne)):.3f}"
          f"（SGD 的终点离最优点还有一点距离，见演示 4）")
    print(f"  J(θ*) = {linreg_cost(th_ne, Xs, y):.1f}，而 σ²/2 = {HOUSING_NOISE_STD ** 2 / 2:.0f}："
          "最小二乘的残差几乎全是噪声，模型没什么可再学的了。")

    # 标准误（standard errors）：Cov(θ̂) = σ² (XᵀX)⁻¹，σ² 用 RSS/(n − d) 估计。在原始单位的设计矩阵上算。
    Xr = add_intercept(X)
    rss = float(np.sum((y - Xr @ fit_least_squares(Xr, y)) ** 2))
    sigma2 = rss / (len(y) - Xr.shape[1])
    se = np.sqrt(sigma2 * np.diag(np.linalg.inv(Xr.T @ Xr)))
    th_orig = to_original_units(th_ne, sc)
    corr = np.corrcoef(X[:, 0], X[:, 1])[0, 1]
    print(f"\n  95% 置信区间 θ̂ ± 1.96·SE（原始单位）；面积和卧室数的相关系数 = {corr:.2f}")
    for name, est, s, truth in zip(["intercept", "area", "bedrooms"], th_orig, se, HOUSING_TRUE_THETA):
        inside = "包含真值" if abs(est - truth) < 1.96 * s else "不包含真值！"
        print(f"    {name:<10} {est:>9.4g} ± {1.96 * s:<9.3g} → [{est - 1.96 * s:>8.4g}, {est + 1.96 * s:>8.4g}]"
              f"   真值 {truth:<5g} {inside}")

    # 这份样本是不是"不走运"？重新生成 2000 份数据，每份都拟合一次、算一次置信区间。
    ests, hits = [], np.zeros(3)
    for rep in range(2000):
        Xb, yb = make_housing(200, np.random.default_rng(10000 + rep))
        Ab = add_intercept(Xb)
        tb = fit_least_squares(Ab, yb)
        rb = yb - Ab @ tb
        seb = np.sqrt(float(rb @ rb) / (len(yb) - 3) * np.diag(np.linalg.inv(Ab.T @ Ab)))
        ests.append(tb)
        hits += np.abs(tb - HOUSING_TRUE_THETA) < 1.96 * seb
    ests = np.array(ests)
    print("\n  重复 2000 次实验（每次重新生成 200 套房子）：")
    print(f"    估计值的平均 = {fmt(ests.mean(axis=0))}   ← 和真值 {fmt(HOUSING_TRUE_THETA)} 几乎一样：无偏（unbiased）")
    print(f"    估计值的标准差 = {fmt(ests.std(axis=0))}   ← 卧室系数的标准差约 4，相对真值 10 很大")
    print(f"    95% 置信区间包含真值的比例 = {fmt(hits / 2000, '.3f')}")
    print("""
观察 Observations：
  · 三种方法到达同一个最优点（J 是凸二次函数，只有一个最低点）。区别只在代价：
    normal equation 一次 O(nd² + d³)；GD 每步 O(nd)；SGD 每步 O(d)，但需要很多步、而且只能"接近"最优点。
  · "卧室"的系数估出来甚至是负的（−3.3，真值 10），面积的系数偏高（0.139，真值 0.12），两个置信区间都没有包含真值。
    这是代码错了吗？重复实验说明没有：估计是无偏的，置信区间 95% 的时候包含真值 —— 这份样本恰好是不走运的
    那一小部分（两个系数都偏了约 3 个标准误，而且方向相反）。
  · 为什么偏得这么"有默契"？卧室数和面积的相关系数是 0.88：数据很难分清"贵是因为面积大还是卧室多"，
    一个系数高估、另一个就低估来补偿。这个方向上的估计方差（estimation variance）很大：
    卧室系数的标准差约 4，而真值才 10。这不是算法的错 —— 三种算法给出同一个数 —— 而是数据本身的信息量不够。
    第 3 讲专门讨论这件事。
  · 所以解读线性回归的单个系数（"多一间卧室值多少钱"）要非常小心。预测（prediction）受的影响小得多：
    只要新房子的"面积–卧室"组合和训练数据相似，两个系数的误差会互相抵消。
  · 标准化之后的 θ 不能直接解读成"每 ft² 多少钱"，要换算回原始单位。""")


# =============================================================================
# 演示 2：学习率的上限 2/λ_max
# =============================================================================
def demo_learning_rate():
    print_header("演示 2：GD 的学习率 —— 收敛的条件是 α < 2/λ_max")
    rng = np.random.default_rng(2)
    X, y = make_housing(200, rng)
    Xs = add_intercept(Standardizer().fit(X).transform(X))
    eig = np.linalg.eigvalsh(Xs.T @ Xs / len(y))          # 从小到大
    lmin, lmax = eig[0], eig[-1]
    J_star = linreg_cost(fit_normal_equation(Xs, y), Xs, y)
    print(f"XᵀX/n 的特征值 eigenvalues = {fmt(eig)}；λ_max = {lmax:.4f}，2/λ_max = {2 / lmax:.4f}\n")
    print(f"{'α (× 1/λ_max)':>14}{'α':>10}{'status':>12}{'iters':>8}{'J − J*':>14}")
    for c in [0.1, 0.5, 1.0, 2 / (1 + lmin / lmax), 1.9, 1.99, 2.01, 2.2]:
        lr = c / lmax
        th, it, status = fit_gd_linear(Xs, y, lr=lr, max_iters=5000, tol=1e-8)
        gap = linreg_cost(th, Xs, y) - J_star
        gap_s = "-" if status == "diverged" else ("≈ 0" if abs(gap) < 1e-9 else f"{gap:.3g}")
        print(f"{c:>14.3f}{lr:>10.4f}{status:>12}{it:>8}{gap_s:>14}")
    print(f"""
观察 Observations：
  · J 是二次函数，沿第 k 个特征方向，每一步误差乘以 (1 − α λ_k)。只要某个方向 |1 − α λ_k| > 1，
    那个方向就指数爆炸 —— 所以上限由最大的特征值 λ_max 决定：α < 2/λ_max。
  · 最快的固定步长是 α = 2/(λ_max + λ_min)（表中第 4 行，{2 / (1 + lmin / lmax):.3f}/λ_max），它让最快和最慢的方向收缩得一样快。
  · α 接近 2/λ_max 时（1.99）虽然收敛，但在最陡的方向上来回震荡（oscillate），反而很慢；
    越过上限一点点（2.01）时发散得也很慢，要上千步 J 才涨到初始值的 10⁸ 倍 —— 所以"loss 慢慢变大"也可能是学习率太大。
  · 实践中不知道 λ_max，所以用 learning-rate sweep：从大往小试，选"刚好不发散"再小一点的。""")


# =============================================================================
# 演示 3：Feature scaling 与条件数
# =============================================================================
def demo_feature_scaling():
    print_header("演示 3：Feature scaling —— 条件数 κ 决定 GD 要走多少步")
    rng = np.random.default_rng(3)
    X, y = make_housing(200, rng)
    versions = {
        "raw (ft², #rooms)": add_intercept(X),
        "area in 1000 ft²": add_intercept(X / np.array([1000.0, 1.0])),
        "standardized": add_intercept(Standardizer().fit(X).transform(X)),
    }
    print(f"{'features':<20}{'λ_min':>10}{'λ_max':>12}{'κ = λmax/λmin':>15}{'GD status':>12}{'iters':>8}{'J − J*':>12}")
    for name, Xv in versions.items():
        eig = np.linalg.eigvalsh(Xv.T @ Xv / len(y))
        J_star = linreg_cost(fit_least_squares(Xv, y), Xv, y)
        th, it, status = fit_gd_linear(Xv, y, lr=1.0 / eig[-1], max_iters=20000, tol=1e-6)
        print(f"{name:<20}{eig[0]:>10.3g}{eig[-1]:>12.3g}{eig[-1] / eig[0]:>15.3g}{status:>12}{it:>8}"
              f"{linreg_cost(th, Xv, y) - J_star:>12.3g}")
    print("""
观察 Observations：
  · 原始特征里"面积"的数量级是 10³，"截距列"是 1：λ_max 被面积主导（≈ 7×10⁶），λ_min 很小，κ ≈ 6×10⁷。
    学习率必须小于 2/λ_max 才不发散，可这么小的步长在"截距"方向上几乎走不动 —— 20000 步后 J 还差得远。
  · 只把面积换成"千平方英尺"，κ 降到约 200，三千多步收敛；完全标准化后 κ ≈ 12，不到 200 步收敛。
  · GD 需要的步数大约 ∝ κ·log(1/ε)。Feature scaling 不改变最优解（只是换了坐标），但能让优化快几个数量级。
    Normal equation / Newton 法不受 κ 影响（它们用了二阶信息），这是第 3 节 Newton 法的伏笔。
  · 深度学习里的 BatchNorm / LayerNorm、Adam 的逐坐标步长，都可以理解为"在线地改善条件数"。""")


# =============================================================================
# 演示 4：SGD 的噪声地板
# =============================================================================
def demo_sgd():
    print_header("演示 4：SGD —— constant step size 停在噪声地板上，decaying step size 才收敛")
    rng = np.random.default_rng(4)
    X, y = make_housing(200, rng)
    Xs = add_intercept(Standardizer().fit(X).transform(X))
    J_star = linreg_cost(fit_normal_equation(Xs, y), Xs, y)
    checkpoints = [1, 5, 20, 100]
    print(f"J* = {J_star:.2f}。表中是 J(θ) − J*（越小越好），每行用同一个打乱顺序的种子：\n")
    print(f"{'schedule':<24}" + "".join(f"{'epoch ' + str(e):>12}" for e in checkpoints)
          + f"{'mean 51-100':>14}{'max 51-100':>13}")
    for name, lr, decay in [("constant α = 0.05", 0.05, False), ("constant α = 0.005", 0.005, False),
                            ("decay α = 0.05/(1+t/n)", 0.05, True)]:
        hist = []
        fit_sgd_linear(Xs, y, lr=lr, epochs=checkpoints[-1], rng=np.random.default_rng(44), decay=decay,
                       history=hist)
        tail = np.array(hist[50:]) - J_star                # 后 50 个 epoch：已经到达"地板"的部分
        print(f"{name:<24}" + "".join(f"{hist[e - 1] - J_star:>12.3g}" for e in checkpoints)
              + f"{tail.mean():>14.3g}{tail.max():>13.3g}")
    print("""
观察 Observations：
  · 常数学习率 0.05：一个 epoch 就到了最优点附近，然后停在一个"噪声地板"上（后 50 个 epoch 平均差 49）
    不再下降 —— 每一步只用一个样本的梯度，它在最优点处也不为 0（单个样本的残差不为 0），θ 一直被推来推去。
  · 常数学习率 0.005：前期慢得多（第 1 个 epoch 后还差 10⁴），但地板低得多（平均 0.8）。这就是 SGD 的基本权衡。
    （经典分析说地板 ∝ α；这里 α 小 10 倍、地板低了约 60 倍，因为我们每个 epoch 随机重排、不放回地取样本，
    并且在 epoch 结尾测量 —— 这种 random reshuffling 的地板下降得比 α 更快。细节不重要，方向才重要。）
  · 学习率衰减兼顾两者：前期大步走，后期小步收敛。深度学习的 learning-rate schedule 也是这个逻辑。
  · 为什么还要用 SGD？当 n 是几十亿时，batch GD 走一步要看完所有数据；SGD（或 mini-batch）看一点就走一步。
    对于"只需要大致最优"的机器学习来说，这笔账非常划算。""")


# =============================================================================
# 演示 5：概率解释 —— 最小二乘 = 高斯 MLE，以及离群点
# =============================================================================
def gaussian_loglik(theta, X, y, sigma):
    """ℓ(θ, σ) = −(n/2) log(2πσ²) − Σ(y − θᵀx)² / (2σ²)：噪声 ε ~ N(0, σ²) 时的对数似然。"""
    r = y - X @ theta
    n = len(y)
    return -0.5 * n * np.log(2 * np.pi * sigma ** 2) - (r @ r) / (2 * sigma ** 2)


def fit_lad_irls(X, y, iters=100, eps=1e-6):
    """
    Least absolute deviations：minimize Σ |y − θᵀx|，即 Laplace 噪声下的最大似然。

    【思路】IRLS（iteratively reweighted least squares）：|r| = r²/|r|，所以在当前残差 r 处，
        L1 损失可以近似成"权重 wᵢ = 1/|rᵢ| 的加权平方损失"。反复：用当前权重解加权最小二乘 → 更新残差 → 更新权重。
    【边界】残差恰好为 0 时 1/|r| = ∞，所以用 max(|r|, eps) 截断。L1 损失在 0 处不可导，没有 normal equation，
        这就是为什么需要迭代算法（实践中用线性规划或者 sklearn 的 QuantileRegressor）。
    """
    theta = fit_least_squares(X, y)
    for _ in range(iters):
        w = 1.0 / np.maximum(np.abs(y - X @ theta), eps)
        sw = np.sqrt(w)
        theta = fit_least_squares(X * sw[:, None], y * sw)
    return theta


def demo_probabilistic():
    print_header("演示 5：概率解释 —— least squares = Gaussian MLE；所以它怕 outliers")
    rng = np.random.default_rng(5)
    X, y = make_housing(200, rng)
    Xs = add_intercept(Standardizer().fit(X).transform(X))
    th = fit_least_squares(Xs, y)
    rss = float(np.sum((y - Xs @ th) ** 2))
    print("(a) 在 θ_LS 附近随机扰动 θ，看高斯对数似然 ℓ(θ, σ) 下降多少（取 3 个不同的 σ）：\n")
    print(f"{'扰动大小 ‖δ‖':>12}" + "".join(f"{'σ = ' + str(s):>16}" for s in [20, 40, 80]))
    drng = np.random.default_rng(55)
    for scale in [0.5, 2.0, 8.0]:
        delta = drng.normal(size=3)
        delta *= scale / np.linalg.norm(delta)
        row = [gaussian_loglik(th + delta, Xs, y, s) - gaussian_loglik(th, Xs, y, s) for s in [20, 40, 80]]
        print(f"{scale:>12.1f}" + "".join(f"{v:>16.3f}" for v in row))
    print(f"\n  σ 的最大似然估计 σ̂ = √(RSS/n) = {np.sqrt(rss / len(y)):.2f}（生成数据用的真值是 {HOUSING_NOISE_STD:g}）")
    print("  每一格都是负数：无论 σ 取多少，ℓ 都在 θ_LS 处最大。σ 只改变曲线的'陡峭程度'（变化量 ∝ 1/σ²），不移动最高点。")

    print("\n(b) 加入 5 个离群点（录入错误、传感器故障……）：y = 2 + 3x + N(0, 0.5²)，n = 50")
    orng = np.random.default_rng(56)
    x = orng.uniform(0, 1, size=50)
    yy = 2 + 3 * x + orng.normal(0, 0.5, size=50)
    x_out = orng.uniform(0.85, 1.0, size=5)
    y_out = 2 + 3 * x_out - 10.0                          # 比正常值低 10，相当于 20 个标准差
    print(f"\n{'data':<26}{'least squares (Gaussian)':>28}{'LAD (Laplace)':>22}")
    for name, xv, yv in [("clean", x, yy), ("+ 5 outliers", np.concatenate([x, x_out]), np.concatenate([yy, y_out]))]:
        Xv = add_intercept(xv)
        print(f"{name:<26}{fmt(fit_least_squares(Xv, yv), '.3f'):>28}{fmt(fit_lad_irls(Xv, yv), '.3f'):>22}")
    print("""                          （表中是 [截距, 斜率]，真值 [2, 3]）

观察 Observations：
  · (a) 说明了讲义 2.4 节的推导：ε ~ N(0, σ²) 独立同分布时，log-likelihood = 常数 − RSS/(2σ²)，
    最大化似然 ⇔ 最小化平方误差。σ 不影响 θ 的估计。
  · 这个视角的价值在于：损失函数不再是"拍脑袋选的"，而是"你对噪声的假设"。换一个假设就换一个损失：
    Laplace 噪声 ⇒ L1 损失（LAD）；Bernoulli ⇒ log loss（第 3 节）；Poisson ⇒ Poisson deviance（第 4 节）。
  · (b) 高斯分布的尾巴很薄：残差 10（20σ）在高斯假设下"几乎不可能"，模型宁可把整条直线拉过去也要减小它
    —— 平方损失下它的代价是 100，正常点只有 0.25。Laplace 的尾巴厚，|r| 只线性增长，几个离群点拉不动直线。
  · 所以"该用什么损失"的问题，等价于"你的数据噪声长什么样"。""")


# =============================================================================
# 演示 6：XᵀX 不可逆
# =============================================================================
def demo_singular():
    print_header("演示 6：XᵀX 不可逆 —— 重复特征与 d > n")
    rng = np.random.default_rng(6)
    X, y = make_housing(200, rng)
    X_dup = np.column_stack([X, X[:, 0] * 0.092903])     # 再加一列"面积（平方米）"：和第一列完全共线
    Xd = add_intercept(Standardizer().fit(X_dup).transform(X_dup))
    Xo = add_intercept(Standardizer().fit(X).transform(X))
    print(f"(a) 特征 = [面积 ft², 卧室, 面积 m²]。XᵀX 的条件数 cond = {np.linalg.cond(Xd.T @ Xd):.3g}")
    try:
        th_ne = fit_normal_equation(Xd, y)
        print(f"    normal equation 没有报错，给出 θ = {fmt(th_ne)}")
        print(f"      两个'面积'系数之和 = {th_ne[1] + th_ne[3]:.2f}，但怎么分是舍入误差随手决定的（换个 BLAS 可能完全不同）")
    except np.linalg.LinAlgError as err:
        print(f"    normal equation 报错：LinAlgError: {err}")
    th_ls = fit_least_squares(Xd, y)
    th_o = fit_least_squares(Xo, y)
    print(f"    lstsq（最小范数解）  θ = {fmt(th_ls)}")
    print(f"    去掉重复列后的解      θ = {fmt(th_o)}")
    print(f"    两者的预测的最大差 = {np.max(np.abs(Xd @ th_ls - Xo @ th_o)):.1e}；"
          f"lstsq 把面积的系数 {th_o[1]:.2f} 平分给了两个重复列（{th_ls[1]:.2f} + {th_ls[3]:.2f}）")

    print("\n(b) d > n：10 个样本、30 个随机特征（更多未知数，更少方程）")
    rng2 = np.random.default_rng(66)
    A = rng2.normal(size=(10, 30))
    b = rng2.normal(size=10)
    th = fit_least_squares(A, b)
    v = np.linalg.svd(A)[2][-1]                           # A 的零空间（null space）里的一个方向：A @ v ≈ 0
    th2 = th + 5.0 * v
    x_new = rng2.normal(size=30)
    print(f"    lstsq 解：训练 MSE = {np.mean((A @ th - b) ** 2):.1e}，‖θ‖ = {np.linalg.norm(th):.3f}")
    print(f"    θ + 5·v（v 在零空间里）：训练 MSE = {np.mean((A @ th2 - b) ** 2):.1e}，‖θ‖ = {np.linalg.norm(th2):.3f}")
    print(f"    两个解在一个新点上的预测：{x_new @ th:.3f} vs {x_new @ th2:.3f}")
    lam_max = np.linalg.eigvalsh(A.T @ A / len(b))[-1]
    th_gd, it, status = fit_gd_linear(A, b, lr=1.0 / lam_max, max_iters=200000, tol=1e-12)
    print(f"    GD 从 θ = 0 出发（{status}，{it} 步）：和 lstsq 解的最大差 = {np.max(np.abs(th_gd - th)):.1e}")
    print("""
观察 Observations：
  · 重复 / 共线特征让 XᵀX 奇异：有一个方向（"ft² 系数加 c，m² 系数减相应的量"）上 J 完全平坦，最优解有无穷多个。
    Normal equation 要么报错（LinAlgError），要么"不报错"但给出由舍入误差决定的任意分法（有时是巨大的、符号相反的数）。
    预测本身没问题，但系数没法解读。所以实践中用 lstsq / QR / SVD，或者加 ridge（第 3 讲）。
  · d > n 时，训练误差可以是 0，而且有无穷多个 θ 都做到 0 —— 它们在训练集上一模一样，在新数据上却不一样。
    "选哪一个"变成了一个真正的问题。lstsq 选范数最小的那个；GD 从 θ = 0 出发也收敛到同一个（上面验证了）——
    因为每一步的梯度 Xᵀ(Xθ − y)/n 都在 Xᵀ 的列空间（行空间）里，θ 永远不会长出零空间的分量。
    这叫 implicit regularization（隐式正则化），是第 3 讲 double descent 的核心。""")


def main():
    demo_three_solvers()
    demo_learning_rate()
    demo_feature_scaling()
    demo_sgd()
    demo_probabilistic()
    demo_singular()


if __name__ == "__main__":
    main()