用 Rust 从头实现机器学习算法·回归②：多元回归、特征缩放与正则化

第 3 篇结尾留过一个承诺：「等进入多元回归，把标量 $x_i$ 换成特征向量、把 Vec 换成 Matrix，梯度公式几乎原封不动」。本篇先兑现它——同一个损失函数、同一次链式法则求导，只是把写法从”逐个样本的循环”升级成”整个数据集的一次矩阵乘法。然后立刻用它做三个对照实验：学习率该取多大、特征量纲悬殊时怎么办、以及给参数戴上”紧箍咒”的 L2 正则化。

实验载体换成一个更像真实任务的玩具问题：房价预测。每条样本有两个量纲悬殊的特征——面积（50–150 ㎡）和房龄（0–30 年），真值模型是

$\text{price} = 0.8 \times \text{面积} - 1.2 \times \text{房龄} + 50 + \varepsilon, \qquad \varepsilon \sim U(-10, 10)$

用与第 3 篇相同的 xorshift 生成器（种子 42）造 300 个样本。面积数值是房龄的十倍量级，这个”别扭”是我们故意安排的——它正是本篇实验二、实验三要解决的矛盾。

核心思想：从一元到多元，公式几乎不变

一元线性回归的预测是 $\hat{y}_i = w x_i + b$：一个特征乘斜率加截距。换成多元，第 $i$ 个样本不再是一个数，而是一行特征 $\boldsymbol{x}_i = (x_{i1}, x_{i2}, \dots, x_{id})$——本篇里 $d = 3$：面积、房龄，外加一个恒为 1 的”虚拟特征”专门给截距用。预测写成

$\hat{y}_i = w_1 x_{i1} + w_2 x_{i2} + \cdots + w_d x_{id} = \boldsymbol{x}_i^{\top} \boldsymbol{w}$

参数也从两个标量 $(w, b)$ 变成一个向量 $\boldsymbol{w} = (w_1, \dots, w_d)^{\top}$。截距去哪了？被吸进了 $\boldsymbol{w}$ 的最后一个分量——只要样本向量最后一位固定填 1，$w_d \cdot 1$ 就是那个常数项。这一招叫截距列，它让”截距”和”斜率”在代码里无需区分。

把 $n$ 个样本的 $\boldsymbol{x}_i^{\top}$ 一行行叠起来，就得到第 2 篇末尾预告过的设计矩阵 $X$（$n \times d$）；$n$ 个真值排成列向量 $\boldsymbol{y}$（$n \times 1$）。至此，“什么叫拟合得好”的度量与第 3 篇逐字相同：

$L(\boldsymbol{w}) = \frac{1}{2n} \sum_{i=1}^{n} \big( \hat{y}_i - y_i \big)^2 = \frac{1}{2n} \lVert X \boldsymbol{w} - \boldsymbol{y} \rVert^2$

看看第 3 篇的训练循环：dw += diff * xs[i]。它在做的事——“每个样本的误差按特征加权后累加”——换成矩阵语言，就是误差向量与 $X$ 的各列分别做点积。承诺里的”几乎原封不动”，下一段推导就能看到有多原封不动。

数学推导

梯度：逐分量求导，再堆成向量

$L$ 对 $w_j$ 的偏导数，和第 3 篇对 $w$ 求导一模一样：外层平方、内层 $\boldsymbol{x}_i^{\top} \boldsymbol{w}$，链式法则剥两层，

$\frac{\partial L}{\partial w_j} = \frac{1}{n} \sum_{i=1}^{n} \big( \hat{y}_i - y_i \big) x_{ij} = \frac{1}{n} \sum_{i=1}^{n} \big( X \boldsymbol{w} - \boldsymbol{y} \big)_i \, X_{ij}$

读法：第 $j$ 个特征的整列 $X_{:j}$，与误差向量逐元素相乘再求和——正是向量点积。把 $j = 1, \dots, d$ 个偏导数从上到下堆成一列，就是那个 $d \times 1$ 的梯度向量：

$\nabla_{\boldsymbol{w}} L = \frac{1}{n} X^{\top} \big( X \boldsymbol{w} - \boldsymbol{y} \big)$

$X^{\top}$ 的作用一句话说清：它把”逐行看样本”的视角转成”逐列看特征”——乘出来的第 $j$ 行，恰好是第 $j$ 个特征列与误差向量的点积。第 3 篇那个内层循环，在这里变成一次矩阵乘法。

更新式：两行并一行

第 3 篇有两个更新式（$w \leftarrow w - \eta\, \partial L / \partial w$，$b \leftarrow b - \eta\, \partial L / \partial b$）。现在它们合并成一个向量更新：

$\boldsymbol{w} \leftarrow \boldsymbol{w} - \eta \, \nabla_{\boldsymbol{w}} L$

形状检查：Xw 是 $n \times 1$，减 $\boldsymbol{y}$ 合法；$X^{\top}$（$d \times n$）乘它得 $d \times 1$，与 $\boldsymbol{w}$ 同形。第 2 篇反复强调的形状规则，在这里是最后一道防线。

损失曲面的形状：为后面两个实验埋的伏笔

$L(\boldsymbol{w})$ 是 $\boldsymbol{w}$ 的二次函数，其二阶导数（Hesse 矩阵）是

$H = \frac{1}{n} X^{\top} X$

等高线是以最低点为中心的椭圆，椭圆的”扁率”由 $H$ 的特征值分布决定。特征值悬殊 → 椭圆狭长，像一条山谷；接近 → 椭圆接近正圆。我们的原始数据里，面积列的平方均值约 $10^4$，而房龄列与截距列张成的方向上特征值不到 $1$——相差四个数量级。这意味着：不同参数方向需要的学习率差了四个数量级，任何单一 $\eta$ 都顾此失彼。实验一会展示”顾此失彼”的惨状，实验二用标准化把山谷变圆。

岭回归：目标函数里加一项

L2 正则化（岭回归）在最小二乘目标上加一个”系数平方和”的惩罚项，强迫系数变小：

$L(\boldsymbol{w}) = \frac{1}{2n} \lVert X \boldsymbol{w} - \boldsymbol{y} \rVert^2 + \frac{\lambda}{2n} \lVert \boldsymbol{w} \rVert^2, \qquad \nabla_{\boldsymbol{w}} L = \frac{1}{n} \Big( X^{\top} \big( X \boldsymbol{w} - \boldsymbol{y} \big) + \lambda \boldsymbol{w} \Big)$

惩罚项对 $\boldsymbol{w}$ 求导就是 $\lambda \boldsymbol{w} / n$——更新时参数多受一股”往零点拽”的力，$\lambda$ 越大拽得越狠。极限行为：$\lambda \to 0$ 退化为普通最小二乘；$\lambda \to \infty$ 时系数全部被压向 0（清单原文）。

两个约定先声明。其一，本篇的 $\lambda$ 采用”除以 $n$“的写法（与损失项同尺度），不同教材约定不同，不影响机制；其二，闭式解 $\hat{\boldsymbol{\beta}} = (X^{\top} X + \lambda I)^{-1} X^{\top} \boldsymbol{y}$ 确实存在——而且 $X^{\top} X + \lambda I$ 一定可逆，天然解决多重共线性（同样是清单原文）——但我们不用它。用梯度下降统一求解，是为了让”正则化 = 梯度里多加一项”这个机制暴露得更彻底；贝叶斯视角（等价于给系数加高斯先验）留给回归⑤。

Rust 实现：数据、矩阵与训练器

完整代码是一个 696 行的 src/main.rs（配套工程在 series/rust-ml/rust-ml-04-gd-tips/），本节起分段展示，各段连起来与文件逐字一致。文件开头是系列共享的迷你 SVG 绘图器——与 series/rust-ml/_plotter-demo/ 里那份相同，本篇额外加了三个小工具：legend（图例）、hline（水平参考线）、note（图上注释）：

// 回归②：多元回归、特征缩放与正则化
// 三个对照实验：① 学习率 ② 特征缩放 ③ L2 正则化（岭回归）。
// 仅依赖标准库；`cargo run` 依次跑完三个实验，并把三张 SVG 图写入
// ../../../frontend/public/images/series/rust-ml-04-gd-tips/

use std::fmt::Write as _;

// ==================== 迷你 SVG 绘图器（系列共享实现） ====================

/// 一张图：坐标映射 + 已累积的 SVG 元素。
struct Canvas {
    w: f64,
    h: f64,
    xmin: f64,
    xmax: f64,
    ymin: f64,
    ymax: f64,
    pad_l: f64,
    pad_r: f64,
    pad_t: f64,
    pad_b: f64,
    body: String,
}

impl Canvas {
    /// 按数据范围自动留 5% 边距建图；xmin == xmax 时手动摊开避免除零。
    fn new(w: f64, h: f64, xmin: f64, xmax: f64, ymin: f64, ymax: f64) -> Canvas {
        let (xmin, xmax) = spread(xmin, xmax);
        let (ymin, ymax) = spread(ymin, ymax);
        Canvas {
            w,
            h,
            xmin,
            xmax,
            ymin,
            ymax,
            pad_l: 52.0,
            pad_r: 16.0,
            pad_t: 16.0,
            pad_b: 40.0,
            body: String::new(),
        }
    }

    fn px(&self, x: f64) -> f64 {
        self.pad_l + (x - self.xmin) / (self.xmax - self.xmin) * (self.w - self.pad_l - self.pad_r)
    }

    fn py(&self, y: f64) -> f64 {
        self.h - self.pad_b - (y - self.ymin) / (self.ymax - self.ymin) * (self.h - self.pad_t - self.pad_b)
    }

    /// 折线，点需按 x 升序传入。
    fn polyline(&mut self, pts: &[(f64, f64)], color: &str, width: f64) {
        let mut d = String::new();
        for (i, (x, y)) in pts.iter().enumerate() {
            let _ = write!(d, "{}{:.1},{:.1}", if i == 0 { "M" } else { "L" }, self.px(*x), self.py(*y));
        }
        let _ = write!(
            self.body,
            r##"<path d="{d}" fill="none" stroke="{color}" stroke-width="{width}" stroke-linejoin="round"/>"##
        );
    }

    /// 散点。（本篇只画折线，此工具留给复现的读者）
    #[allow(dead_code)]
    fn dots(&mut self, pts: &[(f64, f64)], color: &str, r: f64) {
        for (x, y) in pts {
            let _ = write!(
                self.body,
                r##"<circle cx="{:.1}" cy="{:.1}" r="{r}" fill="{color}" fill-opacity="0.65"/>"##,
                self.px(*x),
                self.py(*y)
            );
        }
    }

    /// 坐标轴 + 虚线网格 + 刻度标签。
    fn axes(&mut self, xlabel: &str, ylabel: &str) {
        let x0 = self.px(self.xmin);
        let x1 = self.px(self.xmax);
        let y0 = self.py(self.ymin);
        let y1 = self.py(self.ymax);
        for t in nice_ticks(self.xmin, self.xmax, 6) {
            let x = self.px(t);
            let _ = write!(
                self.body,
                r##"<line x1="{x:.1}" y1="{y0:.1}" x2="{x:.1}" y2="{y1:.1}" stroke="#16161d" stroke-opacity="0.08"/>"##
            );
            let _ = write!(
                self.body,
                r##"<text x="{x:.1}" y="{:.1}" text-anchor="middle" fill="#16161d" fill-opacity="0.45">{t:.6}</text>"##,
                y0 + 16.0,
                t = trim(t)
            );
        }
        for t in nice_ticks(self.ymin, self.ymax, 6) {
            let y = self.py(t);
            let _ = write!(
                self.body,
                r##"<line x1="{x0:.1}" y1="{y:.1}" x2="{x1:.1}" y2="{y:.1}" stroke="#16161d" stroke-opacity="0.08"/>"##
            );
            let _ = write!(
                self.body,
                r##"<text x="{:.1}" y="{:.1}" text-anchor="end" fill="#16161d" fill-opacity="0.45">{t:.6}</text>"##,
                x0 - 6.0,
                y + 4.0,
                t = trim(t)
            );
        }
        let _ = write!(
            self.body,
            r##"<path d="M{x0:.1},{y0:.1}H{x1:.1}M{x0:.1},{y0:.1}V{y1:.1}" fill="none" stroke="#16161d" stroke-opacity="0.6"/>"##
        );
        let _ = write!(
            self.body,
            r##"<text x="{:.1}" y="{:.1}" text-anchor="middle" fill="#16161d" fill-opacity="0.7">{xlabel}</text>"##,
            (x0 + x1) / 2.0,
            self.h - 8.0
        );
        let _ = write!(
            self.body,
            r##"<text x="14" y="{:.1}" text-anchor="middle" transform="rotate(-90 14 {:.1})" fill="#16161d" fill-opacity="0.7">{ylabel}</text>"##,
            (y0 + y1) / 2.0,
            (y0 + y1) / 2.0
        );
    }

    fn save(&self, path: &str) {
        let svg = format!(
            r##"<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {w} {h}" font-family="'JetBrains Mono',monospace" font-size="11"><rect width="{w}" height="{h}" fill="#fbfaf7"/>{body}</svg>"##,
            w = self.w,
            h = self.h,
            body = self.body
        );
        std::fs::write(path, svg).expect("写 SVG 失败");
    }
}

/// 把 fmt 出来的 f64 末尾零去掉，让刻度标签短一点。
fn trim(v: f64) -> String {
    let s = format!("{v:.6}");
    s.trim_end_matches('0').trim_end_matches('.').to_string()
}

/// 数据范围退化成一点时向两侧摊开。
fn spread(min: f64, max: f64) -> (f64, f64) {
    if min == max {
        (min - 1.0, max + 1.0)
    } else {
        let pad = (max - min) * 0.05;
        (min - pad, max + pad)
    }
}

/// “好看”的刻度：步长取 1/2/2.5/5 × 10^k。
fn nice_ticks(min: f64, max: f64, n: usize) -> Vec<f64> {
    let raw = (max - min) / n.max(1) as f64;
    if raw <= 0.0 || !raw.is_finite() {
        return vec![min];
    }
    let exp = raw.abs().log10().floor() as i32;
    let base = 10f64.powi(exp);
    let frac = raw / base;
    let step = if frac <= 1.0 {
        base
    } else if frac <= 2.0 {
        2.0 * base
    } else if frac <= 2.5 {
        2.5 * base
    } else if frac <= 5.0 {
        5.0 * base
    } else {
        10.0 * base
    };
    let mut ticks = Vec::new();
    let mut t = (min / step).ceil() * step;
    while t <= max + 1e-9 {
        ticks.push(if t.abs() < step * 1e-9 { 0.0 } else { t });
        t += step;
    }
    ticks
}

/// 图例：在数据坐标 (x, y) 起，逐行画「色线 + 文字」。
fn legend(c: &mut Canvas, x: f64, y: f64, entries: &[(&str, &str)]) {
    let x0 = c.px(x);
    let mut yy = c.py(y);
    for (color, text) in entries {
        let _ = write!(
            c.body,
            r##"<line x1="{x0:.1}" y1="{yy:.1}" x2="{:.1}" y2="{yy:.1}" stroke="{color}" stroke-width="2.5"/>"##,
            x0 + 28.0
        );
        let _ = write!(
            c.body,
            r##"<text x="{:.1}" y="{yy:.1}" dominant-baseline="middle" fill="#16161d" fill-opacity="0.8">{text}</text>"##,
            x0 + 34.0
        );
        yy += 17.0;
    }
}

/// 在 y 处画一条贯穿绘图区的水平虚线（参考线），右侧加标签。
fn hline(c: &mut Canvas, y: f64, label: &str) {
    let x0 = c.px(c.xmin);
    let x1 = c.px(c.xmax);
    let yy = c.py(y);
    let _ = write!(
        c.body,
        r##"<line x1="{x0:.1}" y1="{yy:.1}" x2="{x1:.1}" y2="{yy:.1}" stroke="#16161d" stroke-opacity="0.35" stroke-dasharray="5,4"/>"##
    );
    let _ = write!(
        c.body,
        r##"<text x="{:.1}" y="{:.1}" text-anchor="end" fill="#16161d" fill-opacity="0.45">{label}</text>"##,
        x1 - 4.0,
        yy - 5.0
    );
}

接下来是第 3 篇用过的 xorshift（原样保留，种子仍取 42）和第 1、2 篇 Matrix 的最小实现——matmul/transpose 逐字来自第 2 篇，另加 add/sub/scale/sq_norm 四个小运算，支撑梯度公式里出现的”乘标量、加减向量”：

// ==================== 伪随机数（与第 3 篇相同） ====================

/// xorshift64 伪随机数生成器：无需外部依赖，结果可复现。
struct XorShift(u64);

impl XorShift {
    fn new(seed: u64) -> Self {
        assert!(seed != 0, "种子不能为 0");
        XorShift(seed)
    }

    fn next_u64(&mut self) -> u64 {
        let mut x = self.0;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.0 = x;
        x
    }

    /// 返回 [0, 1) 区间的 f64。
    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }
}

// ==================== Matrix 最小实现（第 1、2 篇） ====================

/// 行优先存储的二维浮点矩阵：第 i 行第 j 列存放在 data[i * cols + j]。
#[derive(Debug, Clone)]
struct Matrix {
    data: Vec<f64>,
    rows: usize,
    cols: usize,
}

impl Matrix {
    /// 全零矩阵。
    fn zeros(rows: usize, cols: usize) -> Matrix {
        Matrix { data: vec![0.0; rows * cols], rows, cols }
    }

    /// 用一维向量按行优先构造矩阵。
    fn from_vec(data: Vec<f64>, rows: usize, cols: usize) -> Matrix {
        assert_eq!(data.len(), rows * cols, "数据长度与矩阵形状不一致");
        Matrix { data, rows, cols }
    }

    /// 读取第 i 行第 j 列的元素。
    fn get(&self, i: usize, j: usize) -> f64 {
        self.data[i * self.cols + j]
    }

    /// 矩阵乘法：(n×m) · (m×p) → (n×p)，逐字翻译 C_ij = Σ_k A_ik B_kj。
    fn matmul(&self, other: &Matrix) -> Matrix {
        assert_eq!(
            self.cols, other.rows,
            "内维不一致：{}×{} 无法乘 {}×{}",
            self.rows, self.cols, other.rows, other.cols
        );
        let mut out = Matrix::zeros(self.rows, other.cols);
        for i in 0..self.rows {
            for j in 0..other.cols {
                let mut sum = 0.0;
                for k in 0..self.cols {
                    sum += self.data[i * self.cols + k] * other.data[k * other.cols + j];
                }
                out.data[i * out.cols + j] = sum;
            }
        }
        out
    }

    /// 转置：n×m → m×n。
    fn transpose(&self) -> Matrix {
        let mut out = Matrix::zeros(self.cols, self.rows);
        for i in 0..self.rows {
            for j in 0..self.cols {
                out.data[j * self.rows + i] = self.data[i * self.cols + j];
            }
        }
        out
    }

    /// 同形矩阵逐元素相加。
    fn add(&self, other: &Matrix) -> Matrix {
        assert_eq!(self.rows, other.rows);
        assert_eq!(self.cols, other.cols);
        let data: Vec<f64> = self
            .data
            .iter()
            .zip(other.data.iter())
            .map(|(a, b)| a + b)
            .collect();
        Matrix { data, rows: self.rows, cols: self.cols }
    }

    /// 同形矩阵逐元素相减。
    fn sub(&self, other: &Matrix) -> Matrix {
        assert_eq!(self.rows, other.rows);
        assert_eq!(self.cols, other.cols);
        let data: Vec<f64> = self
            .data
            .iter()
            .zip(other.data.iter())
            .map(|(a, b)| a - b)
            .collect();
        Matrix { data, rows: self.rows, cols: self.cols }
    }

    /// 全体元素乘一个标量。
    fn scale(&self, k: f64) -> Matrix {
        let data: Vec<f64> = self.data.iter().map(|a| a * k).collect();
        Matrix { data, rows: self.rows, cols: self.cols }
    }

    /// 全体元素的平方和（列向量 n×1 时即 ‖v‖²）。
    fn sq_norm(&self) -> f64 {
        self.data.iter().map(|a| a * a).sum()
    }
}

造数据、算均值方差、拼设计矩阵。standardize 就是 $z = (x - \mu)/\sigma$，逐元素 map 一遍即可——mean/std_dev/standardize 三行的写法正是本节后面”Rust 语法角”的主角：

// ==================== 数据：房价预测 ====================

/// 生成 n 条房价样本。真值模型：price = 0.8·面积 − 1.2·房龄 + 50 + ε，
/// 面积 ∈ [50, 150)（㎡），房龄 ∈ [0, 30)（年），ε ~ U(-10, 10)。
/// 面积与房龄量纲悬殊（数值范围差一个数量级），正是本篇的实验素材。
fn gen_data(n: usize, seed: u64) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
    let mut rng = XorShift::new(seed);
    let mut area = Vec::with_capacity(n);
    let mut age = Vec::with_capacity(n);
    let mut price = Vec::with_capacity(n);
    for _ in 0..n {
        let a = 50.0 + rng.next_f64() * 100.0; // 面积
        let g = rng.next_f64() * 30.0;         // 房龄
        let eps = (rng.next_f64() - 0.5) * 20.0; // 噪声 ε ∈ [-10, 10)
        area.push(a);
        age.push(g);
        price.push(0.8 * a - 1.2 * g + 50.0 + eps);
    }
    (area, age, price)
}

// ==================== 特征缩放 ====================

/// 均值：所有元素求和除以个数。
fn mean(xs: &[f64]) -> f64 {
    xs.iter().sum::<f64>() / xs.len() as f64
}

/// 标准差（总体标准差，除以 n）。
fn std_dev(xs: &[f64], mu: f64) -> f64 {
    let var = xs.iter().map(|x| { let d = x - mu; d * d }).sum::<f64>() / xs.len() as f64;
    var.sqrt()
}

/// 标准化：z = (x − μ) / σ。
fn standardize(xs: &[f64], mu: f64, sigma: f64) -> Vec<f64> {
    xs.iter().map(|x| (x - mu) / sigma).collect()
}

/// 拼 n×3 设计矩阵：每一行是 [特征1, 特征2, 1]，最后一列是截距项的“系数 1”。
fn design_matrix(f1: &[f64], f2: &[f64]) -> Matrix {
    let n = f1.len();
    let mut data = Vec::with_capacity(n * 3);
    for i in 0..n {
        data.push(f1[i]);
        data.push(f2[i]);
        data.push(1.0);
    }
    Matrix::from_vec(data, n, 3)
}

最后是训练器 train：把上一段的梯度公式逐符号翻译成矩阵运算，λ 和衰减都做成参数——三个实验共用这同一段代码，改动只在传入的参数。mse 是评估用的均方误差，restore 把标准化空间学到的系数换算回原始量纲，后面验证”学对了没有”要靠它：

// ==================== 训练器：全批量梯度下降 ====================

/// 在 X（n×d）、y（n×1）上训练线性回归，目标函数
///   L(w) = ‖Xw − y‖² / (2n) + λ·‖w‖² / (2n)    （λ = 0 即普通最小二乘）
/// 梯度 ∇w = ( Xᵀ(Xw − y) + λw ) / n，每轮更新 w ← w − ηₜ·∇w。
/// decay = Some((eta0, tau)) 时学习率随轮数衰减：ηₜ = eta0 / (1 + t/tau)。
/// 返回 (最终 w, 每轮损失)。损失一旦出现 NaN/inf（数值发散）就提前停。
fn train(
    x: &Matrix,
    y: &Matrix,
    eta0: f64,
    epochs: usize,
    lambda: f64,
    decay: Option<(f64, f64)>,
) -> (Matrix, Vec<f64>) {
    let n = x.rows;
    let mut w = Matrix::zeros(x.cols, 1);
    let mut losses = Vec::with_capacity(epochs);
    let xt = x.transpose(); // 训练前转置一次，循环内复用
    for t in 0..epochs {
        let eta = match decay {
            Some((_, tau)) => eta0 / (1.0 + t as f64 / tau),
            None => eta0,
        };
        let diff = x.matmul(&w).sub(y); // Xw − y（n×1）
        let grad = xt
            .matmul(&diff)
            .add(&w.scale(lambda)) // Xᵀ(Xw−y) + λw（λ = 0 时加零向量）
            .scale(1.0 / n as f64);
        let mut loss = diff.sq_norm() / (2.0 * n as f64);
        if lambda > 0.0 {
            loss += lambda * w.sq_norm() / (2.0 * n as f64);
        }
        w = w.sub(&grad.scale(eta));
        if !loss.is_finite() {
            break; // 已发散，后面的轮次没有意义
        }
        losses.push(loss);
    }
    (w, losses)
}

/// 均方误差（不含正则项）：MSE = ‖Xw − y‖² / n。
fn mse(x: &Matrix, y: &Matrix, w: &Matrix) -> f64 {
    let diff = x.matmul(&w).sub(y);
    diff.sq_norm() / x.rows as f64
}

/// 把标准化空间学到的系数还原成原始量纲：
/// 斜率 w_orig = w_std/σ；截距补回标准化时减掉的 μ。
fn restore(w: &Matrix, mu1: f64, sd1: f64, mu2: f64, sd2: f64) -> (f64, f64, f64) {
    let w1 = w.get(0, 0);
    let w2 = w.get(1, 0);
    let w3 = w.get(2, 0);
    (w1 / sd1, w2 / sd2, w3 - w1 * mu1 / sd1 - w2 * mu2 / sd2)
}

实现要点：

• xt 在循环外转置一次、循环内反复复用——$X^{\top}$ 不随 $t$ 变，这是”矩阵化”相比第 3 篇循环累加的另一重好处：重复计算被显式暴露出来。
• 梯度链 matmul → add(λw) → scale(1/n) 与公式 $\big( X^{\top}(X \boldsymbol{w} - \boldsymbol{y}) + \lambda \boldsymbol{w} \big) / n$ 逐符号对应。$\lambda = 0$ 时 .add(&w.scale(0.0)) 加的是零向量，不必为此分支。
• !loss.is_finite() 防御发散：学习率太大时 $w$ 会指数级冲出 f64 的表示范围变成 inf/NaN，提前 break 让曲线停在”爆炸前夜”。
• 对照第 3 篇的 dw += diff * xs[i]：那一行累加在多元版里整个消失了——它正是 $X^{\top}(X\boldsymbol{w}-\boldsymbol{y})$ 的逐分量手工版。

分析与画图的小工具先备齐——log_curve 把损失序列变成对数坐标折线点（等距下落代表指数级改善），clip_y 把冲出画布的曲线截在上缘（视觉上表示”已发散”），decimate 给几万轮的长曲线抽稀，epochs_to_reach 定量回答”谁更快”：

/// 把每轮损失转成 (epoch, log10 loss) 折线点；log 纵轴下等距下落 = 指数级改善。
fn log_curve(losses: &[f64]) -> Vec<(f64, f64)> {
    losses
        .iter()
        .enumerate()
        .filter(|(_, l)| l.is_finite() && **l > 0.0)
        .map(|(t, l)| (t as f64 + 1.0, l.log10()))
        .collect()
}

/// 把折线点的 y 值截断到上限：冲出去的曲线贴在上缘，表示“已发散”。
fn clip_y(pts: &[(f64, f64)], ymax: f64) -> Vec<(f64, f64)> {
    pts.iter().map(|&(x, y)| (x, y.min(ymax))).collect()
}

/// 在数据坐标 (x, y) 处画一行灰色小字注释。
fn note(c: &mut Canvas, x: f64, y: f64, text: &str) {
    let _ = write!(
        c.body,
        r##"<text x="{:.1}" y="{:.1}" fill="#16161d" fill-opacity="0.6">{text}</text>"##,
        c.px(x),
        c.py(y)
    );
}

/// 损失序列首次降到 threshold 以下所用的轮数（用于“谁更快”的定量对比）。
fn epochs_to_reach(losses: &[f64], threshold: f64) -> Option<usize> {
    losses.iter().position(|&l| l < threshold).map(|i| i + 1)
}

/// 抽稀折线点：长曲线（几万轮）画出来太密，每 step 个点留 1 个，首尾保留。
fn decimate(pts: &[(f64, f64)], step: usize) -> Vec<(f64, f64)> {
    if pts.len() <= 2 || step <= 1 {
        return pts.to_vec();
    }
    let mut out: Vec<(f64, f64)> = pts.iter().step_by(step).copied().collect();
    if out.last() != pts.last() {
        out.push(*pts.last().unwrap());
    }
    out
}

实验一：学习率——太大震荡，太小磨蹭

先看不做任何缩放的原始数据。$w$ 从 0 出发，三种固定学习率各跑 3000 轮，再加一条”衰减曲线”作对照。loss 的”地板”是多少？噪声 $\varepsilon$ 的方差是 $20^2/12 \approx 33.3$，而 $L$ 带一个 $1/2$，所以地板 ≈ $16.7$——和第 3 篇”loss 停在噪声地板”的结论一致。

// ==================== main：三个实验 ====================

fn main() {
    let out_dir = "../../../frontend/public/images/series/rust-ml-04-gd-tips";
    std::fs::create_dir_all(out_dir).expect("创建输出目录失败");

    // ---- 数据：300 条房价样本（种子 42，每次运行结果完全一致） ----
    let n = 300usize;
    let (area, age, price) = gen_data(n, 42);
    let y = Matrix::from_vec(price.clone(), n, 1);
    let mu_a = mean(&area);
    let sd_a = std_dev(&area, mu_a);
    let mu_g = mean(&age);
    let sd_g = std_dev(&age, mu_g);
    println!("=== 数据 ===");
    println!("n = {}，真值模型 price = 0.8·面积 − 1.2·房龄 + 50 + ε（ε ~ U(-10,10)）", n);
    println!("面积: μ = {:.2}, σ = {:.2}；房龄: μ = {:.2}, σ = {:.2}", mu_a, sd_a, mu_g, sd_g);

    // ---- 实验一：学习率（原始量纲，不缩放） ----
    println!("\n=== 实验一：学习率（未标准化数据，各跑 3000 轮） ===");
    let x_raw = design_matrix(&area, &age);
    let mut curves: Vec<(Vec<(f64, f64)>, &str)> = Vec::new();
    for (eta, color, label) in [
        (0.0001, "#0f766e", "η = 0.0001"),
        (0.001, "#d97706", "η = 0.001"),
        (0.01, "#dc2626", "η = 0.01"),
    ] {
        let (_, losses) = train(&x_raw, &y, eta, 3000, 0.0, None);
        let final_loss = *losses.last().unwrap();
        if losses.len() < 3000 {
            println!("{}：第 {} 轮数值发散（此前 loss 已到 {:.1e}），提前中断", label, losses.len() + 1, final_loss);
        } else {
            println!("{}：没爆，但 3000 轮后 loss 仍有 {:.4}（地板 ≈ 16.7）", label, final_loss);
        }
        curves.push((log_curve(&losses), color));
    }
    // 学习率衰减：固定步长的 η₀ = 0.001 本来发散，让它按 ηₜ = η₀/(1+t/3) 快速衰减
    let epochs1 = 3000usize;
    let (_, decay_losses) = train(&x_raw, &y, 0.001, epochs1, 0.0, Some((0.001, 3.0)));
    println!(
        "{}：η₀ = 0.001 按 1/(1+t/3) 衰减，{}",
        "带衰减",
        if decay_losses.len() < epochs1 {
            format!("第 {} 轮仍发散", decay_losses.len() + 1)
        } else {
            format!("最终 loss = {:.4}", decay_losses.last().unwrap())
        }
    );
    curves.push((log_curve(&decay_losses), "#7c3aed"));

    let ymax_log = 7.9; // 纵轴截断：冲过 log10(loss)=7.9 的曲线贴上缘，表示已发散
    let mut c = Canvas::new(640.0, 400.0, 0.0, epochs1 as f64, 0.0, 8.0);
    c.axes("epoch", "log10(loss)");
    for (pts, color) in &curves {
        c.polyline(&clip_y(&decimate(pts, 2), ymax_log), color, 2.0);
    }
    hline(&mut c, (100.0f64 / 12.0 / 2.0).log10(), "loss 地板 ≈ 16.7");
    note(&mut c, 120.0, 7.3, "红、橙两线约 20~150 轮间冲出上缘 = 已发散");
    note(&mut c, 60.0, 5.6, "紫线前 ~15 轮曾冲至 5×10^14，随后俯冲回落");
    legend(&mut c, 2100.0, 7.5, &[
        ("#0f766e", "η = 0.0001"),
        ("#d97706", "η = 0.001"),
        ("#dc2626", "η = 0.01"),
        ("#7c3aed", "η=0.001 衰减 τ=3"),
    ]);
    c.save(&format!("{out_dir}/lr-compare.svg"));

运行输出（种子固定，每次完全一致）：

=== 实验一：学习率（未标准化数据，各跑 3000 轮） ===
η = 0.0001：没爆，但 3000 轮后 loss 仍有 99.8545（地板 ≈ 16.7）
η = 0.001：第 155 轮数值发散（此前 loss 已到 3.3e304），提前中断
η = 0.01：第 76 轮数值发散（此前 loss 已到 1.0e303），提前中断
带衰减：η₀ = 0.001 按 1/(1+t/3) 衰减，最终 loss = 104.0749

三条固定步长的曲线给出一个尴尬的局面：$\eta = 0.01$ 和 $\eta = 0.001$ 在第 76、155 轮先后数值爆炸（loss 冲到 $10^{303}$ 以上，f64 再装不下就是 inf）；$\eta = 0.0001$ 不炸，但 3000 轮后还停在 $99.85$——离地板 $16.7$ 遥遥无期。夹在三个数量级之间，没有一个 $\eta$ 是好选择。

这不是 $\eta$ 的错，是山谷的错。前面推导过，Hesse 矩阵 $X^{\top}X/n$ 的特征值从约 $10^4$（面积方向）一路跨到 $1$ 以下（房龄与截距张成的方向），差了四个数量级。梯度下降的收敛要求 $\eta$ 小于 $2/\lambda_{\max} \approx 0.0002$，否则面积方向上的每步更新反而把误差放大——这就是爆炸；而房龄、截距方向上的有效步长是 $\eta \lambda_{\min}$，小得可怜——这就是磨蹭。一个标量 $\eta$ 要同时伺候四个数量级，顾此失彼是必然的。

学习率衰减是实践中救急的标准动作：让步长随轮数缩小，前期大步冲、后期小步磨。紫线演示的就是它——$\eta_0 = 0.001$ 原本必爆，按 $\eta_t = \eta_0 / (1 + t/3)$ 衰减后，前十几轮也曾冲到 $5 \times 10^{14}$，但步长迅速降下来，硬是把自己摁回了下降通道，最终稳在 $104$。不过请注意：它只是”不爆”，并没有比固定小步长更快摸到地板——慢方向依然慢。治本的办法，是把山谷本身变圆。

实验二：特征缩放——把狭长的山谷变圆

几何直觉：损失等高线是椭圆，梯度方向垂直于等高线。在山谷又窄又长时，垂直于谷壁的梯度几乎指向谷壁两侧，下降一步就在两侧来回弹跳，前进缓慢；把特征标准化到同一尺度后，椭圆被”掰圆”，任何一点的梯度都直指圆心，每步都是有效前进。标准化就是两个公式：

$\mu = \frac{1}{n} \sum_{i=1}^{n} x_i, \qquad \sigma = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (x_i - \mu)^2}, \qquad z_i = \frac{x_i - \mu}{\sigma}$

Rust 实现就是前面见过的 mean/std_dev/standardize 三个一行函数。实验对比三条曲线：未缩放配 $\eta = 0.01$（会爆）、未缩放退而用 $\eta = 0.0001$（给足 20000 轮）、标准化后配 $\eta = 0.01$（只跑 2000 轮）：

    // ---- 实验二：特征缩放 ----
    println!("\n=== 实验二：特征缩放 ===");
    println!("同一 η = 0.01：");
    let (_, losses_raw_big) = train(&x_raw, &y, 0.01, 2000, 0.0, None);
    println!(
        "  未标准化：{}",
        if losses_raw_big.len() < 2000 {
            format!("第 {} 轮数值发散", losses_raw_big.len() + 1)
        } else {
            format!("收敛，最终 loss = {:.4}", losses_raw_big.last().unwrap())
        }
    );
    // 未缩放时能把步长压到不爆的 η = 0.0001，给它 20000 轮看多久能摸到地板
    let (_, losses_raw_small) = train(&x_raw, &y, 0.0001, 20000, 0.0, None);
    let t_raw = epochs_to_reach(&losses_raw_small, 20.0);
    let t_raw_txt = match t_raw {
        Some(t) => format!("{} 轮", t),
        None => "没做到".to_string(),
    };
    println!(
        "  未缩放（改用 η = 0.0001 跑 20000 轮）：最终 loss = {:.4}，降到 20 以下用了 {}",
        losses_raw_small.last().unwrap(),
        t_raw_txt
    );
    // 标准化 z = (x − μ)/σ（此处 300 条全用于训练，统计量取自这 300 条本身）
    let z_area = standardize(&area, mu_a, sd_a);
    let z_age = standardize(&age, mu_g, sd_g);
    let x_std = design_matrix(&z_area, &z_age);
    let (w_std, losses_std) = train(&x_std, &y, 0.01, 2000, 0.0, None);
    let t_std = epochs_to_reach(&losses_std, 20.0);
    let t_std_txt = match t_std {
        Some(t) => format!("{} 轮", t),
        None => "没做到".to_string(),
    };
    println!(
        "  标准化后（η = 0.01，2000 轮）：最终 loss = {:.4}，降到 20 以下用了 {}",
        losses_std.last().unwrap(),
        t_std_txt
    );
    // 把标准化空间学到的系数还原成原始量纲，应当接近真值 [0.8, -1.2, 50]
    let (o1, o2, o3) = restore(&w_std, mu_a, sd_a, mu_g, sd_g);
    println!("  还原为原始量纲：[面积 {:.4}, 房龄 {:.4}, 截距 {:.4}]", o1, o2, o3);

    let sc_curves = [
        log_curve(&losses_raw_big),
        decimate(&log_curve(&losses_raw_small), 20),
        decimate(&log_curve(&losses_std), 5),
    ];
    let mut c = Canvas::new(640.0, 400.0, 0.0, 20000.0, 0.0, 8.0);
    c.axes("epoch", "log10(loss)");
    c.polyline(&clip_y(&sc_curves[0], 7.9), "#dc2626", 2.0);
    c.polyline(&clip_y(&sc_curves[1], 7.9), "#0f766e", 2.0);
    c.polyline(&clip_y(&sc_curves[2], 7.9), "#2563eb", 2.0);
    hline(&mut c, (100.0f64 / 12.0 / 2.0).log10(), "loss 地板 ≈ 16.7");
    note(&mut c, 3000.0, 2.6, "蓝线在 370 轮处击穿 20（左缘的陡降）");
    note(&mut c, 300.0, 7.3, "红线 ~76 轮冲出上缘 = 已发散");
    legend(&mut c, 13300.0, 7.5, &[
        ("#dc2626", "未缩放, η=0.01"),
        ("#0f766e", "未缩放, η=0.0001"),
        ("#2563eb", "标准化, η=0.01"),
    ]);
    c.save(&format!("{out_dir}/scaling.svg"));

输出：

=== 实验二：特征缩放 ===
同一 η = 0.01：
  未标准化：第 76 轮数值发散
  未缩放（改用 η = 0.0001 跑 20000 轮）：最终 loss = 82.9877，降到 20 以下用了 没做到
  标准化后（η = 0.01，2000 轮）：最终 loss = 15.9008，降到 20 以下用了 370 轮
  还原为原始量纲：[面积 0.7962, 房龄 -1.2808, 截距 51.4568]

同一 $\eta = 0.01$：未缩放的红线 76 轮爆炸；标准化后的蓝线却只用 370 轮就把 loss 击穿 20，2000 轮后停在地板 $15.9$。而那个”安全”的小步长青线，给了十倍轮数（20000 轮）最终还在 $82.99$，“降到 20 以下”这一项直接没做到——慢了近两个数量级且依然遥遥无期。条件数把”难收敛”量化了出来：标准化前 $X^{\top}X/n$ 的特征值跨四个数量级，标准化后各特征列方差都为 1、且中心化后的特征列与截距列恰好正交（$\sum_i (x_i - \mu) \cdot 1 = 0$），Hesse 矩阵近似单位阵，山谷被掰圆，同一个 $\eta = 0.01$ 从”必然爆炸”变成”一路狂飙”。

最后那行”还原为原始量纲”是正确性自检：$[0.7962, -1.2808, 51.4568]$ 对上真值 $[0.8, -1.2, 50]$（不完全相等正是噪声的本职工作）。还原公式就是标准化的逆变换：斜率除以 $\sigma$，截距补回被减掉的 $\mu$。

一个必须养成的纪律：测试集要用训练集的 $\mu$/$\sigma$ 来变换，绝不允许用全体数据（更不允许用测试集）重新估计——那等于偷看了答案，验证误差会系统性偏低。实验二里 300 条全用于训练、统计量来自这 300 条本身，是因为这里没有留出测试集；下一节的实验三有训练/验证分割，你会看到 $\mu$、$\sigma$ 只从训练集那 200 条估计，再拿去变换验证集。

实验三：L2 正则化——给参数戴上紧箍咒

岭回归对训练器的改动小得惊人：梯度里多一项 $\lambda \boldsymbol{w}$，目标函数里多一项 $\lambda \lVert \boldsymbol{w} \rVert^2 / (2n)$——train 函数早就内置了。实验设计：前 200 条训练、后 100 条验证；$\mu$、$\sigma$ 只用训练集估计；扫 $\lambda \in \{0.1, 1, 10, 100, 1000, 10000\}$ 画图，另取 $\lambda \in \{0, 1, 10, 100\}$ 列表打印系数：

    // ---- 实验三：L2 正则化（岭回归） ----
    println!("\n=== 实验三：L2 正则化（岭回归） ===");
    let n_train = 200usize;
    let area_tr = &area[..n_train];
    let area_va = &area[n_train..];
    let age_tr = &age[..n_train];
    let age_va = &age[n_train..];
    let price_tr = &price[..n_train];
    let price_va = &price[n_train..];
    // 关键纪律：μ、σ 只用训练集估计，再拿去变换验证集
    let (mu_a_tr, mu_g_tr) = (mean(area_tr), mean(age_tr));
    let (sd_a_tr, sd_g_tr) = (std_dev(area_tr, mu_a_tr), std_dev(age_tr, mu_g_tr));
    println!(
        "训练/验证 = {}/{}；μ、σ 只用训练集估计（面积 μ={:.2} σ={:.2}，房龄 μ={:.2} σ={:.2}）",
        n_train, n - n_train, mu_a_tr, sd_a_tr, mu_g_tr, sd_g_tr
    );
    let x_tr = design_matrix(
        &standardize(area_tr, mu_a_tr, sd_a_tr),
        &standardize(age_tr, mu_g_tr, sd_g_tr),
    );
    let x_va = design_matrix(
        &standardize(area_va, mu_a_tr, sd_a_tr),
        &standardize(age_va, mu_g_tr, sd_g_tr),
    );
    let y_tr = Matrix::from_vec(price_tr.to_vec(), n_train, 1);
    let y_va = Matrix::from_vec(price_va.to_vec(), n - n_train, 1);

    // λ 扫描（画图用；λ 采用除以 n 的约定，λ = 0.1 已近似 OLS）
    let mut ridge_pts: Vec<(f64, f64, f64)> = Vec::new(); // (log10 λ, log10 训练MSE, log10 验证MSE)
    for &lam in &[0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0] {
        let (w, _) = train(&x_tr, &y_tr, 0.01, 2000, lam, None);
        ridge_pts.push((lam.log10(), mse(&x_tr, &y_tr, &w).log10(), mse(&x_va, &y_va, &w).log10()));
    }

    println!("{:<6} | {:^24} | {:^34} | {:>9} | {:>9}", "λ", "w(标准化空间)", "还原为原始量纲 [面积, 房龄, 截距]", "训练MSE", "验证MSE");
    for &lam in &[0.0, 1.0, 10.0, 100.0] {
        let (w, _) = train(&x_tr, &y_tr, 0.01, 2000, lam, None);
        let (o1, o2, o3) = restore(&w, mu_a_tr, sd_a_tr, mu_g_tr, sd_g_tr);
        println!(
            "{:<6} | [{:>7.3}, {:>7.3}, {:>7.3}] | [{:>7.3}, {:>7.3}, {:>7.3}] | {:>9.2} | {:>9.2}",
            lam,
            w.get(0, 0), w.get(1, 0), w.get(2, 0),
            o1, o2, o3,
            mse(&x_tr, &y_tr, &w), mse(&x_va, &y_va, &w)
        );
    }

    let y_min = ridge_pts
        .iter()
        .flat_map(|p| [p.1, p.2])
        .fold(f64::INFINITY, f64::min)
        - 0.3;
    let y_max = ridge_pts
        .iter()
        .flat_map(|p| [p.1, p.2])
        .fold(f64::NEG_INFINITY, f64::max)
        + 0.6;
    let mut c = Canvas::new(640.0, 400.0, -2.0, 3.0, y_min, y_max);
    c.axes("log10(λ)", "log10(MSE)");
    let train_pts: Vec<(f64, f64)> = ridge_pts.iter().map(|p| (p.0, p.1)).collect();
    let valid_pts: Vec<(f64, f64)> = ridge_pts.iter().map(|p| (p.0, p.2)).collect();
    c.polyline(&train_pts, "#2563eb", 2.0);
    c.polyline(&valid_pts, "#d6491f", 2.0);
    hline(&mut c, (400.0f64 / 12.0).log10(), "噪声地板 σ²≈33.3");
    note(&mut c, -1.85, 3.3, "两线几乎重合：本例样本充足、特征近似正交，");
    note(&mut c, -1.85, 3.1, "λ=0 已接近最优，正则化的主场在共线/高维场景");
    legend(&mut c, 2.5, 1.75, &[
        ("#2563eb", "训练 MSE"),
        ("#d6491f", "验证 MSE"),
    ]);
    c.save(&format!("{out_dir}/ridge.svg"));

    println!("\n三张图已写入 {out_dir}/（lr-compare.svg / scaling.svg / ridge.svg）");
}

输出：

=== 实验三：L2 正则化（岭回归） ===
训练/验证 = 200/100；μ、σ 只用训练集估计（面积 μ=98.06 σ=27.80，房龄 μ=13.97 σ=8.27）
λ      |         w(标准化空间)         |        还原为原始量纲 [面积, 房龄, 截距]        |     训练MSE |     验证MSE
0      | [ 22.526, -10.545, 111.585] | [  0.810,  -1.275,  49.931] |     33.00 |     29.91
1      | [ 22.407, -10.480, 111.030] | [  0.806,  -1.267,  49.686] |     33.33 |     30.10
10     | [ 21.390,  -9.927, 106.272] | [  0.770,  -1.200,  47.580] |     62.76 |     58.22
100    | [ 14.726,  -6.473,  74.390] | [  0.530,  -0.782,  33.373] |   1487.10 |   1466.16

先验证底子：$\lambda = 0$ 时还原系数 $[0.810, -1.275, 49.931]$，再次对上真值 $[0.8, -1.2, 50]$。然后看 $\lambda$ 的作用机制——系数被整体朝零收缩：$\lambda = 1$ 时几乎没动；$\lambda = 10$ 时面积系数从 $0.81$ 收到 $0.77$；$\lambda = 100$ 时收到 $0.53$，截距从 $49.9$ 压到 $33.4$，训练 MSE 从 $33$ 抬到 $1487$。这就是”紧箍咒”：$\lambda$ 越大，模型越”保守”，对训练数据越将就不了——偏差换方差。图上左端（$\lambda \le 1$）曲线贴着噪声地板 $\sigma^2 \approx 33.3$，右端随 $\lambda$ 抬升而上翘，翘得越狠说明欠拟合越重。

也请诚实地看到另一件事：本例中训练与验证曲线几乎重合，验证误差并没有出现一个”先降后升”的 U 形——因为样本相对充足、标准化后的特征近似正交（Hesse 矩阵接近单位阵，方差无从压缩），$\lambda = 0$ 已接近最优，正则化能省的方差寥寥无几、引入的偏差却实打实。岭回归真正的主场是多重共线性强（$X^{\top}X$ 病态，闭式解里 $X^{\top}X + \lambda I$ 必可逆的优势开始值钱）与高维小样本的场景，这正是清单 1.6 把”适合共线性强的数据”写进优点的原因；那里验证曲线会弯出明显的 U 底，留待回归④用 L1/L2 对比展开。

实务上 $\lambda$ 怎么选？用验证集（或交叉验证）扫出来——就像这个实验一样，取验证误差的最低点。这套”扫超参 → 验证集选点”的流程，到回归③给多项式回归选阶数时还要再用一次。

Rust 语法角：迭代器、collect 与 as 类型转换

实验二里那三个小函数是 Rust 数据处理的典型写法，此处把 mean、std_dev 两个核心再引一遍（与上文相同，方便对照），逐个拆解：

/// 均值：所有元素求和除以个数。
fn mean(xs: &[f64]) -> f64 {
    xs.iter().sum::<f64>() / xs.len() as f64
}

/// 标准差（总体标准差，除以 n）。
fn std_dev(xs: &[f64], mu: f64) -> f64 {
    let var = xs.iter().map(|x| { let d = x - mu; d * d }).sum::<f64>() / xs.len() as f64;
    var.sqrt()
}

迭代器与 sum。 xs.iter() 不复制数据，只借出一个”挨个访问元素”的迭代器；.sum::<f64>() 把它消耗的求和写出来。<f64> 是类型标注——因为空列表的和可以是任意数值类型，Rust 需要知道累加器用哪种类型（sum::<f64>() 是 turbofish 写法，写成 let s: f64 = xs.iter().sum(); 让编译器从变量类型推断也行）。对照 Python：sum(xs) 是内置函数，Rust 的 sum 是迭代器上的方法——差异的根源是 Rust 把”遍历并归约”统一抽象在迭代器 trait 上，sum/product/count/map/filter 都是这一层的方法。《Rust 程序设计语言》第 8 章向量一节讲了 Vec 与 iter() 的关系，第 13 章的迭代器一节把这套惰性管线讲透了——值得整章读，机器学习代码的一半都写在迭代器上。

map + collect 造新向量。 standardize 里的 xs.iter().map(|x| (x - mu) / sigma).collect() 对应 Python 的列表推导 [(x - mu) / sigma for x in xs]：惰性变换逐元素生成，最后由 collect 收口成具体的 Vec<f64>。std_dev 里的 map(|x| { let d = x - mu; d * d }) 演示了闭包里写多行语句——花括号包起来、最后一条表达式就是返回值，和函数体同一套规则。

as 类型转换。 xs.len() as f64 里藏着本篇最容易踩的坑：len() 返回 usize，直接除以 f64 不允许（Rust 不做隐式数值转换），必须显式 as。它对应 Python 的 float(len(xs))——但 as 是运算符不是函数调用，而且可以串在表达式中间：xs.len() as f64 读起来几乎像类型注解。精度上，usize 是 64 位整数而 f64 尾数 53 位，天文数字（超过 $2^{53}$）会舍入；数据规模这种量级的 n 转换精确无损。反向转换（f64 as usize）则会截断小数，用时要更谨慎。

后缀字面量。 代码里还有两处第 3 篇就见过、但值得再点一次的写法：300usize 给整数字面量标注类型，1u64 << 53 说明这是 64 位无符号运算。它们与 as 配合，让”哪个数是什么类型”在代码层面一目了然——这正是第 1 篇说的”让错误死在编译期”的具体样子。

运行结果汇总

三个实验在同一次 cargo run 里按序完成（完整输出如下，种子固定、可复现）：

=== 数据 ===
n = 300，真值模型 price = 0.8·面积 − 1.2·房龄 + 50 + ε（ε ~ U(-10,10)）
面积: μ = 98.09, σ = 28.22；房龄: μ = 14.14, σ = 8.09

=== 实验一：学习率（未标准化数据，各跑 3000 轮） ===
η = 0.0001：没爆，但 3000 轮后 loss 仍有 99.8545（地板 ≈ 16.7）
η = 0.001：第 155 轮数值发散（此前 loss 已到 3.3e304），提前中断
η = 0.01：第 76 轮数值发散（此前 loss 已到 1.0e303），提前中断
带衰减：η₀ = 0.001 按 1/(1+t/3) 衰减，最终 loss = 104.0749

=== 实验二：特征缩放 ===
同一 η = 0.01：
  未标准化：第 76 轮数值发散
  未缩放（改用 η = 0.0001 跑 20000 轮）：最终 loss = 82.9877，降到 20 以下用了 没做到
  标准化后（η = 0.01，2000 轮）：最终 loss = 15.9008，降到 20 以下用了 370 轮
  还原为原始量纲：[面积 0.7962, 房龄 -1.2808, 截距 51.4568]

=== 实验三：L2 正则化（岭回归） ===
训练/验证 = 200/100；μ、σ 只用训练集估计（面积 μ=98.06 σ=27.80，房龄 μ=13.97 σ=8.27）
λ      |         w(标准化空间)         |        还原为原始量纲 [面积, 房龄, 截距]        |     训练MSE |     验证MSE
0      | [ 22.526, -10.545, 111.585] | [  0.810,  -1.275,  49.931] |     33.00 |     29.91
1      | [ 22.407, -10.480, 111.030] | [  0.806,  -1.267,  49.686] |     33.33 |     30.10
10     | [ 21.390,  -9.927, 106.272] | [  0.770,  -1.200,  47.580] |     62.76 |     58.22
100    | [ 14.726,  -6.473,  74.390] | [  0.530,  -0.782,  33.373] |   1487.10 |   1466.16

三张图已写入 ../../../frontend/public/images/series/rust-ml-04-gd-tips/（lr-compare.svg / scaling.svg / ridge.svg）

复现方式：

cd series/rust-ml/rust-ml-04-gd-tips
cargo run

三张 SVG（lr-compare.svg 学习率对比、scaling.svg 缩放前后、ridge.svg 的 λ 扫描）会写入相对工程目录的 ../../../frontend/public/images/series/rust-ml-04-gd-tips/，正文以上方的方式引用。

怎么读这些图

三张图都用 log 纵轴（log10(loss) 或 log10(MSE)），先建立读感：

• 等距下落 = 指数级改善。纵轴上每掉一格代表缩小 10 倍，直线式下落就是教科书式的指数收敛；变平则说明接近极限（地板或噪声）。
• 冲出上缘 = 已发散。纵轴截断在 8（loss 图），冲顶的曲线被 clip_y 摁在上缘——贴着上缘的横线不是收敛，是”已经飞出画布”。
• 水平虚线 = 地板。loss ≈ 16.7（噪声方差的一半）、MSE ≈ 33.3（噪声方差）是任何模型都无法逾越的极限，曲线贴近它就该知足。

具体到每张图：

• lr-compare.svg：看三件事——红、橙两线的”上升沿”（那就是发散的指数上冲）；青色线的”缓坡”（3000 轮只降了不到两个数量级，这是磨蹭）；紫线的”冲顶—俯冲”（衰减自救：前十几轮冲至 $5 \times 10^{14}$，步长降下来后被摁回下降通道）。四条线合在一起回答：固定步长在未缩放数据上没有好选择，衰减只救”爆”不救”慢”。
• scaling.svg：横轴拉到 20000 轮就为了看清青线有多慢——它几乎水平地爬完全程；蓝线在左缘 370 轮处近乎垂直地击穿 20，然后就贴在地板上了。同一 $\eta = 0.01$，红爆、青爬、蓝飞，唯一的变量是特征缩不缩放。
• ridge.svg：横轴是 $\log_{10}\lambda$（$\lambda = 0$ 无法取对数，图中左端 $\lambda = 0.1$ 已近似 OLS，精确值看上面输出表）。左端两条线都贴着噪声地板（模型 fitting 到位）；向右抬升越快，说明偏差引入得越狠。训练、验证两线几乎重合是本例的诚实结论——没有过拟合可供压制，$\lambda = 0$ 已接近最优；等回归③、④遇到共线的升维特征，这里才会弯出 U 底。

优缺点与适用场景

本篇主角是岭回归（清单 1.6），原文照录：

核心思想：在最小二乘目标上加一个「系数平方和」的惩罚项 λΣβⱼ²，强迫系数变小。

数学要点：目标 min ‖y - Xβ‖² + λ‖β‖²；闭式解 β̂ = (XᵀX + λI)⁻¹Xᵀy ——注意 XᵀX + λI 一定可逆，天然解决多重共线性；λ→0 退化为 OLS；λ→∞ 系数趋近 0；贝叶斯视角：等价于给系数加高斯先验。

优点：稳定、系数都保留（只是缩小）、适合共线性强的数据。

缺点：不能做特征选择（系数不会精确为 0）。

清单没有为岭回归单列适用场景，但它的优点即场景：共线性强、又要求系数全部保留的数据（光谱、经济指标这类高度相关的特征群）。本系列后文会多次回到它：回归④把它与 L1、弹性网络摆在一起对比特征选择能力，回归⑤会从贝叶斯角度重新推出它。

下一篇的主角多项式回归（清单 1.4）也先混个脸熟，原文照录：

核心思想：把特征升维（x, x², x³, …），让线性模型拟合曲线。

优点：简单即可引入非线性。

缺点：过拟合的典型教材案例——阶数过高时曲线在样本间剧烈震荡；高阶项之间高度相关，引发数值不稳定。

关键认知：多项式回归是理解「模型复杂度 vs 泛化能力」的最佳入口，必须配合正则化或交叉验证选择阶数。

注意这条关键认知里”必须配合的正则化”——就是本篇戴上的紧箍咒；而本篇实验三末尾演示的”扫超参、验证集选点”，正是”交叉验证选择阶数”的雏形。三篇内容在这里闭环。

至于本篇一直默默服役的多元 OLS 本身，清单 1.1 给它的定位是基线模型、以及”需要解释系数的场景（经济学、社会科学）“——本例的系数天然可解释：每平米 $+0.8$、每多一年房龄 $-1.2$、基准 $+50$。

小结与下篇预告

兑现承诺的部分：一元到多元，$x_i$ 换成特征向量、Vec 换成 Matrix，损失函数 $\frac{1}{2n}\lVert X\boldsymbol{w} - \boldsymbol{y} \rVert^2$ 与梯度 $\frac{1}{n} X^{\top}(X\boldsymbol{w} - \boldsymbol{y})$ 的写法原封不动；第 3 篇内层循环里的 dw += diff * xs[i] 升维后整个消失，化作一次矩阵乘法。更重要的是此后的所有实验——学习率、衰减、标准化、$\lambda$——都只是同一个 train 函数的不同参数。

三个技巧的收束：学习率决定每一步迈多大，太大在山谷里弹跳爆炸、太小原地磨蹭，衰减能救命但救不了慢；特征缩放把狭长的山谷掰圆，标准化后 Hesse 矩阵近似单位阵，同一个学习率从必爆变成狂飙，并要牢记”测试集用训练集的 $\mu$/$\sigma$“的纪律；L2 正则化用 $\lambda$ 把系数往零点拽，$\lambda$ 的选择交给验证集扫出来的曲线。三个实验共用 300 条种子固定的房价样本、696 行零依赖的 main.rs，每一行都能在本地复现。

下一篇进入回归③：多项式回归与过拟合。我们将把唯一的特征升维成 $(x, x^2, x^3, \dots)$，看线性模型如何拟出曲线，亲眼看清单那句”曲线在样本间剧烈震荡”长什么样——再用本篇的交叉验证流程给它选阶数、用紧箍咒给它收骨头。