线性回归：第一个完整模型

在导读里，我们曾经给打印耗时挑过一条直线：把几个候选参数挨个试一遍，看看谁预测得更准。现在有了梯度，就不必再靠枚举寻找参数了。

不过，学会更新参数只完成了一半工作。数据怎样分？输入要不要缩放？参数越大越好吗？选好模型以后，怎样判断它对新任务有没有用？这一章把这些问题放进同一个小项目。

1. 先说清楚：预测时能看到什么

任务仍然是：收到一份文档，在打印开始前预测耗时。输入页数 $x$ 叫作特征，实际耗时 $y$ 叫作标签。模型输出 $\hat y$，帽子提醒我们“这是预测值”。

我们沿用导读的教学数据，单位分别是页和秒：

训练集用来调整参数；验证集用来比较方案；测试集留到方案确定之后。它们只有几条记录，适合把计算看清楚，不足以估计真实打印服务的可靠性。

“预测时能看到什么”比选算法更早。页数在开始前可知，可以当特征；任务完成时间虽然与耗时关系密切，却要等打印结束才知道。把后者放进输入，练出来的是一道已经泄露答案的题。

先设一个很朴素的对手：无论多少页，都猜训练耗时的平均数 72.5 秒。这个基线在验证集上的平均绝对误差为 21.167 秒。后面的模型如果连它都赢不了，就值得回头检查数据和实现。

2. 换一把尺子，不换预测问题

直接写 $\hat y=ax+b$ 没有错。不过，页数在 5～30 之间，乘上页数的梯度与固定项的梯度尺度不同，选学习率会稍麻烦一些。

先只用训练页数，计算中心和尺度：

$$\mu=\frac1n\sum_{i=1}^n x_i,\qquad s=\sqrt{\frac1n\sum_{i=1}^n(x_i-\mu)^2}$$

本例 $\mu=17.5$，$s\approx8.539$。把页数改写为 $z=(x-\mu)/s$，训练：

$$\hat y=wz+b$$

这里的 $s$ 用 $n$ 作分母，目的是把训练特征缩放到方便计算的尺度；第四章用 $n-1$ 估计总体方差，解决的是另一个问题。

验证、测试和以后新来的任务，都沿用训练时的 $\mu$ 和 $s$。不能每来一批数据就重新计算一把尺子，否则同样的 12 页会变成不同的 $z$，原有参数也就接不上了。

想把结果换回页数单位，展开即可：

$$w\frac{x-\mu}{s}+b =\frac ws x+\left(b-\frac{w\mu}{s}\right)$$

所以原始斜率是 $a=w/s$，截距是 $b-w\mu/s$。标准化没有让直线变成别的形状，只是换了坐标。

3. 让参数沿着误差往下走

把第 $i$ 条记录的预测误差记为 $e_i=wz_i+b-y_i$。训练目标仍然是均方误差：

$$L(w,b)=\frac1n\sum_{i=1}^n e_i^2$$

按照第二章的链式法则：

$$\begin{aligned} \frac{\partial L}{\partial w} &=\frac2n\sum_{i=1}^n e_i z_i\\ \frac{\partial L}{\partial b} &=\frac2n\sum_{i=1}^n e_i \end{aligned}$$

每一轮先用同一组旧参数算出全部误差，再一起更新 $w,b$。如果先更新 $w$，再用新 $w$ 算 $b$ 的梯度，就不是这里推导的这一步了。

图 1：同一份训练数据上的优化过程。曲线下降说明训练目标被改善，还没有回答新任务预测得怎样。

这组数据最后得到：

$$\hat y=3x+20$$

于是 24 页预测为 92 秒。斜率 3 表示：在这个模型中，多一页，预测增加 3 秒。截距 20 是直线在零页处的取值。我们没有观测零页任务，也没有做因果实验，不能单凭这条线断言“真实设备固定启动时间就是 20 秒”。

4. 为什么有时要约束参数

假如有几百个特征，却只有几十条记录，模型可能用很大的正负权重互相抵消，把训练数据中的偶然波动也追得很紧。给权重加一点约束，常常能让拟合更稳一些。

本章采用 L2 正则化：

$$J(w,b)=L(w,b)+\lambda w^2,\qquad\lambda\ge0$$

平方项会惩罚较大的权重。这里不惩罚固定项 $b$，否则数据整体偏高也会受到额外惩罚。新梯度只需给 $w$ 的梯度加上 $2\lambda w$；$b$ 的梯度保持原样。

$w,b$ 从训练数据中学出来，叫模型参数；$\lambda$ 决定怎样训练，叫超参数。我们事先列出 0、0.1、1 三个候选值，用验证 MAE 比较它们。

图 2：本例三个模型都通过训练数据中心；增加正则化使斜率缩小，但未必改善验证误差。

这一次 $\lambda=0$ 最好。这完全合理：正则化是一种可供检验的约束，不是每次加上都会涨分的开关。

多特征时，惩罚项通常变为 $\lambda\sum_j w_j^2$。这也解释了为什么尺度很重要：同一个物理量从“页”换成“千页”，所需权重会放大一千倍；若直接惩罚权重，单位选择就会改变约束的力度。标准化能帮助我们用较一致的尺度比较权重，但不会自动发现有用特征。

5. 完整运行一次

下载完整脚本。它只使用 Python 标准库，在仓库根目录运行：

python3 src/examples/linear_regression.py

"""第 6 章：从零训练带 L2 正则项的一元线性回归；只用标准库。"""

from math import sqrt
from statistics import mean


# 与导读使用同一份教学数据，按同一固定划分进行实验。
TRAIN = [(5, 37), (10, 47), (15, 66), (20, 81), (25, 92), (30, 112)]
VALIDATION = [(8, 43), (18, 76), (28, 103)]
TEST = [(12, 58), (22, 85), (26, 101)]


def predict(model, pages):
    z = (pages - model["center"]) / model["scale"]
    return model["weight"] * z + model["bias"]


def mae(model, data):
    return mean(abs(predict(model, x) - y) for x, y in data)


def mse(model, data):
    return mean((predict(model, x) - y)**2 for x, y in data)


def fit(data, regularization=0.0, rate=0.1, steps=300):
    center = mean(x for x, _ in data)
    scale = sqrt(mean((x-center)**2 for x, _ in data))
    if scale == 0:
        raise ValueError("训练页数没有变化，无法据此确定斜率")
    z_values = [(x-center)/scale for x, _ in data]
    weight, bias = 0.0, 0.0
    model = {"center": center, "scale": scale, "weight": weight, "bias": bias}
    history = [mse(model, data)]
    for _ in range(steps):
        errors = [weight*z + bias - y for z, (_, y) in zip(z_values, data)]
        dw = 2*mean(e*z for e, z in zip(errors, z_values)) + 2*regularization*weight
        db = 2*mean(errors)
        weight, bias = weight-rate*dw, bias-rate*db
        model.update(weight=weight, bias=bias)
        history.append(mse(model, data) + regularization*weight**2)
    return model, history


def raw_parameters(model):
    slope = model["weight"] / model["scale"]
    intercept = model["bias"] - slope * model["center"]
    return slope, intercept


def main():
    baseline = mean(y for _, y in TRAIN)
    baseline_mae = mean(abs(baseline-y) for _, y in VALIDATION)
    print(f"均值基线的验证 MAE：{baseline_mae:.3f} 秒")
    candidates = []
    for strength in [0.0, 0.1, 1.0]:
        model, _ = fit(TRAIN, regularization=strength)
        score = mae(model, VALIDATION)
        candidates.append((score, strength, model))
        slope, intercept = raw_parameters(model)
        print(f"λ={strength:.1f}：斜率={slope:.4f}，截距={intercept:.4f}，验证 MAE={score:.3f}")
    _, strength, selected = min(candidates, key=lambda item: item[0])
    print(f"选择 λ={strength:.1f}；参数到此固定")
    print(f"训练 MSE：{mse(selected, TRAIN):.3f} 秒²")
    print(f"测试 MAE：{mae(selected, TEST):.3f} 秒")
    for x, y in TEST:
        predicted = predict(selected, x)
        print(f"{x} 页：预测 {predicted:.1f}，实际 {y}，误差 {predicted-y:+.1f} 秒")
    print(f"24 页的预测耗时：{predict(selected, 24):.1f} 秒")


if __name__ == "__main__":
    main()

输出为：

均值基线的验证 MAE：21.167 秒
λ=0.0：斜率=3.0000，截距=20.0000，验证 MAE=1.333
λ=0.1：斜率=2.7273，截距=24.7727，验证 MAE=2.530
λ=1.0：斜率=1.5000，截距=46.2500，验证 MAE=10.917
选择 λ=0.0；参数到此固定
训练 MSE：4.667 秒²
测试 MAE：2.000 秒
12 页：预测 56.0，实际 58，误差 -2.0 秒
22 页：预测 86.0，实际 85，误差 +1.0 秒
26 页：预测 98.0，实际 101，误差 -3.0 秒
24 页的预测耗时：92.0 秒

训练使用 MSE，报告使用 MAE，并不矛盾。MSE 光滑，较大误差受到更重惩罚；MAE 与耗时同为秒，更容易向使用者解释。训练 MSE 的单位是秒的平方，不能直接拿 4.667 与 MAE 的 2.000 比大小。

保存模型时，除了权重和固定项，也要保存中心、尺度、特征顺序及单位。只留两个参数，却忘记输入怎样转换，恢复出来的就不是同一个预测过程了。

6. 平均分之外，再看看错在哪里

三条测试误差分别为 −2、+1、−3 秒，负号表示低估，正号表示高估。平均绝对值是 2 秒，但平均数没有告诉我们错误是否集中在某一类任务上。

如果长文档总被低估，也许关系开始弯曲；如果彩印任务普遍更慢，也许缺了打印模式这个特征。先观察错误，再决定增加什么信息，比盲目增加模型复杂度更有方向。

训练页数最大只有 30。拿直线去预测 3000 页属于外推：公式当然能算，但没有数据保证同样的关系仍然成立。设备缓存、暂停散热和补纸，都可能改变规律。

还有一个容易忽略的极端：如果训练记录全是 10 页，我们无法从中分辨斜率与截距。很多不同直线都能给 10 页相同的预测。程序在尺度为零时主动报错，是因为数据没提供足够信息；继续除以零并不能补上这些信息。

多特征模型可以写成 $\hat y=\boldsymbol w^\mathsf T\boldsymbol x+b$，整批预测为：

$$\hat{\boldsymbol y}=X\boldsymbol w+b\boldsymbol1$$

$X$ 的形状是 $n\times d$，每行一个样本、每列一个特征；$\boldsymbol w$ 有 $d$ 个分量。第一章的矩阵运算，到这里就成为实际模型的计算方式了。

7. 自测与小实验

1. 训练、验证、测试三组数据分别负责什么？
2. 模型对 22 页预测 86 秒，实际为 85 秒。误差、平方误差、绝对误差分别是多少？
3. 若训练中心为 10、尺度为 2，模型为 $\hat y=6z+50$，换回页数后是什么？
4. 为什么不能用测试 MAE 选择 $\lambda$，然后还把同一个 MAE 当作独立评估？
5. 本例正则化越强，验证误差越大，这说明正则化没有用吗？
6. 将学习率改成 0.001，仍训练 300 步。结果变差一定是模型表达能力不够吗？
7. 给模型增加“实际完成时间”特征，训练误差大幅下降。为什么仍不能用于开始前预测？

选做：自己用训练数据计算一元最小二乘的闭式解，和梯度训练比较。令 $\bar x,\bar y$ 为训练均值：

$$a=\frac{\sum_i(x_i-\bar x)(y_i-\bar y)} {\sum_i(x_i-\bar x)^2}, \qquad b=\bar y-a\bar x$$

本例应得到 $a=3,b=20$。闭式解给了我们独立核对训练程序的办法；多特征时，实际软件通常使用数值线性代数方法求解，未必要迭代走梯度。

接下来，如果问题变成“这个任务会超时吗”，输出就不再是一段耗时，而是一个类别。我们将在分类中给直线的输出加上概率解释。

想了解常用实现，可查阅 scikit-learn 的线性模型文档。比较正则化参数时要同时查看损失的定义：有的实现对误差求和，有的取平均，参数数值不能直接照搬。