Day 06｜线性回归手写实现:从零拼出第一个 ML 模型（AI 学习笔记 · 基础筑基周 · 第 6 篇）

机器学习的最小闭环只有四件套——假设函数、损失函数、梯度计算、参数更新。线性回归是这四件套最干净的载体:$\hat y = w^\top x + b$ 是假设,MSE 是损失,$\nabla L$ 是梯度,$w \leftarrow w - \eta \nabla L$ 是更新。今天用纯 numpy 把这套闭环从零写一遍,再用 sklearn 验证权重一致,亲眼看到 loss 曲线收敛。

1. 线性回归:模型假设与几何意义

1.1 定义

线性回归(Linear Regression)假设目标 $y$ 与特征 $x$ 之间是线性关系:

\hat y = w_1 x_1 + w_2 x_2 + \dots + w_n x_n + b = \mathbf{w}^\top \mathbf{x} + b

其中 $\mathbf{w} \in \mathbb{R}^n$ 是权重向量(每个特征一个),$b \in \mathbb{R}$ 是偏置(intercept)。预测是特征加权和加偏置。线性指的是 $w$ 是线性的,$x$ 可以是任何函数——所以「用 $x^2$ 当特征」也叫线性回归(对参数线性)。

1.2 公式 / 矩阵形式

\begin{aligned}
\text{单样本:}\quad &\hat y_i = \mathbf{w}^\top \mathbf{x}_i + b \\
\text{全部样本:}\quad &\hat{\mathbf{y}} = X \mathbf{w} + b \cdot \mathbf{1} \\
&\text{其中 } X \in \mathbb{R}^{N \times n},\ \hat{\mathbf{y}} \in \mathbb{R}^N
\end{aligned}

加入偏置列后可以写成 $\hat{\mathbf{y}} = X’ \boldsymbol{\beta}$,其中 $X’ = [X\ \mathbf{1}]$,$\boldsymbol{\beta} = [\mathbf{w}; b]$,形式更紧凑。

1.3 例子

用「房屋面积」预测「成交价」。真实数据 $(x_i, y_i)$,我们想拟合 $\hat y = w \cdot x + b$,$w$ ≈ 每平米单价,$b$ ≈ 底价(地段附加)。

import numpy as np
rng = np.random.default_rng(42)
x = rng.uniform(30, 150, 50)             # 面积 30-150 平米
y = 8.5 * x + 120 + rng.normal(0, 30, 50) # 单价 8.5k/平米,底价 120k,噪声 σ=30k

线性回归要做的是:在 $w$-$b$ 平面上找一个点,让所有 $(x_i, y_i)$ 到直线 $\hat y = wx + b$ 的纵向距离平方和最小。

1.4 取舍 / 几何直觉

把 $N$ 个数据点放在二维平面上,$\hat y = w x + b$ 是一条直线,线性回归就是找一条「最小二乘」直线——让所有点到直线的纵向距离平方和最小。几何上等价于把这些点向直线做正交投影,这就是为什么叫”最小二乘 / least squares”。多维情形推广为「超平面」,依然是最近距离平方和最小。

1.5 对比表

2. MSE 损失与最小二乘

2.1 定义

均方误差(Mean Squared Error, MSE)是回归任务最常用的损失:

L(\mathbf{w}, b) = \frac{1}{N} \sum_{i=1}^{N} (\hat y_i - y_i)^2 = \frac{1}{N}\lVert \hat{\mathbf{y}} - \mathbf{y} \rVert_2^2

衡量”预测值与真实值的平均平方偏离”。平方放大误差(MSE 是 MAE 的「对大误差更敏感」版本),数学上处处可导,与高斯噪声下的极大似然估计等价。

2.2 公式 / 推导(MSE 的概率解释)

假设残差 $\epsilon_i = y_i - \hat y_i$ 独立同分布 $\mathcal{N}(0, \sigma^2)$,则 $y_i \mid \mathbf{x}_i \sim \mathcal{N}(\mathbf{w}^\top \mathbf{x}_i + b, \sigma^2)$。极大似然估计等价于最小化负对数似然:

\begin{aligned}
\log P(\mathbf{y} \mid X; \mathbf{w}, b) &= \sum_i \log \mathcal{N}(y_i; \hat y_i, \sigma^2) \\
&= -\frac{1}{2\sigma^2} \sum_i (y_i - \hat y_i)^2 + \text{const} \\
\Rightarrow \hat{\mathbf{w}}_{\text{MLE}} &= \arg\max \log P = \arg\min \sum_i (y_i - \hat y_i)^2 \\
&= \arg\min L(\mathbf{w}, b)
\end{aligned}

MSE = 高斯噪声 + 极大似然——这就是为什么线性回归默认配 MSE。

2.3 例子

预测 100、真实 80,误差 $-20$,平方 400;预测 60、真实 80,误差 $+20$,平方也是 400——MSE 对正负误差对称,所以线性回归不会”偏向预测高或低”,而是用中位数或均值中心对齐。

2.4 取舍 / 直觉

平方把大误差”放大”、小误差”压缩”,给离群点更大权重——好处是模型被迫照顾大误差样本,坏处是一个离谱的异常值(如录入错误的房价)会把整条直线”拽歪”。抗异常值应改用 Huber Loss(大误差改 MAE)或 MAE。

2.5 对比表

3. 梯度下降与解析解(正规方程)

3.1 定义

梯度下降(Gradient Descent, GD)用负梯度方向迭代更新参数:

\mathbf{w} \leftarrow \mathbf{w} - \eta \frac{\partial L}{\partial \mathbf{w}}

$\eta > 0$ 是学习率(learning rate),控制每步步长。

正规方程(Normal Equation)是 MSE 凸优化的解析解:

\hat{\mathbf{w}}_{\text{OLS}} = (X^\top X)^{-1} X^\top \mathbf{y}

(包含偏置时,给 $X$ 增广一列 1)。两者解相同,但成本 / 适用场景差异大。

3.2 公式 / 推导(MSE 梯度)

把 MSE 写成矩阵形式:

L = \frac{1}{N}\lVert X\mathbf{w} - \mathbf{y} \rVert_2^2 = \frac{1}{N}(X\mathbf{w} - \mathbf{y})^\top(X\mathbf{w} - \mathbf{y})

对 $\mathbf{w}$ 求导(用 $\frac{\partial \mathbf{a}^\top \mathbf{b}}{\partial \mathbf{a}} = \mathbf{b}$ 与链式法则):

\begin{aligned}
\frac{\partial L}{\partial \mathbf{w}} &= \frac{2}{N} X^\top (X\mathbf{w} - \mathbf{y}) \\
\text{令} = \mathbf{0} \Rightarrow X^\top X \mathbf{w} &= X^\top \mathbf{y} \\
\Rightarrow \mathbf{w}^* &= (X^\top X)^{-1} X^\top \mathbf{y}
\end{aligned}
$$

GD 视角:每次 $\mathbf{w} \leftarrow \mathbf{w} - \eta \cdot \frac{2}{N} X^\top(X\mathbf{w} - \mathbf{y})$,从任意起点沿梯度下山。

### 3.3 例子:固定数据下三种策略对比

```python
import numpy as np
rng = np.random.default_rng(0)
X = rng.normal(size=(200, 5))
true_w = np.array([2, -1, 3, 0.5, -0.7])
y = X @ true_w + 0.1 * rng.normal(size=200)

# 1) 正规方程
w_ne = np.linalg.solve(X.T @ X, X.T @ y)

# 2) 梯度下降
w_gd = np.zeros(5)
eta, N = 0.1, len(y)
for _ in range(2000):
    grad = 2 * X.T @ (X @ w_gd - y) / N
    w_gd -= eta * grad

# 3) sklearn
from sklearn.linear_model import LinearRegression
w_sk = LinearRegression(fit_intercept=False).fit(X, y).coef_

print(f"NE:    {np.round(w_ne, 4)}")
print(f"GD:    {np.round(w_gd, 4)}")
print(f"SKL:   {np.round(w_sk, 4)}")
print(f"|NE-GD|={np.linalg.norm(w_ne - w_gd):.2e}, |NE-SKL|={np.linalg.norm(w_ne - w_sk):.2e}")

预期:三组权重完全相同,误差在 1e-6 量级。

3.4 取舍 / 几何直觉

MSE 在 $\mathbf{w}$ 上是凸二次函数(因为 $X^\top X$ 半正定),全局唯一极小点就是解析解。GD 沿梯度下山,只要 $\eta$ 不太大就一定能收敛到这个点。

把 MSE 想象成”碗状曲面”:碗底 = 最优参数。GD 就像蒙眼下山,每步找最陡方向挪一小步;正规方程像戴透视眼一眼看到碗底。

3.5 对比表

4. 学习率、收敛与向量化

4.1 定义

• 学习率 $\eta$:GD 每步的步长因子。太大 → 震荡甚至发散(L 增大);太小 → 收敛极慢,卡在高原。
• 批量梯度下降(BGD):每次用全量 $N$ 个样本算梯度,稳定但慢。
• 随机 GD (SGD):每次只用 1 个样本,快但震荡大。
• Mini-batch GD:每次用 $B$ 个样本(典型 $B=32$ / $64$ / $128$),平衡稳定与速度,深度学习默认。
• 向量化(Vectorization):用矩阵运算替代 Python for 循环,利用 BLAS SIMD 加速 50-200 倍。

4.2 公式 / 学习率与收敛

GD 更新 $\mathbf{w}_{t+1} = \mathbf{w}_t - \eta \nabla L$ 在 $\mathbf{w}^*$ 附近泰勒展开:

L(\mathbf{w}_{t+1}) \approx L(\mathbf{w}_t) + \nabla L^\top \Delta \mathbf{w} + \frac{1}{2} \Delta \mathbf{w}^\top H \Delta \mathbf{w}

其中 $\Delta \mathbf{w} = -\eta \nabla L$,$H = \nabla^2 L$ 是 Hessian。要让 $L$ 下降需要:

\eta < \frac{2}{\lambda_{\max}(H)}

超出此值,二阶项主导,Loss 反向增大——这就是”$\eta$ 太大发散”的数学原因。Adam / RMSProp 等自适应方法本质是在每步动态估计合适的 $\eta$。

4.3 例子:不同学习率的收敛轨迹

import numpy as np
import matplotlib.pyplot as plt

# f(x) = (x-3)^2 + 1,梯度 2(x-3),最优点 x*=3
f = lambda x: (x - 3) ** 2 + 1
grad = lambda x: 2 * (x - 3)

def gd(x0, eta, n_steps):
    xs = [x0]
    for _ in range(n_steps):
        xs.append(xs[-1] - eta * grad(xs[-1]))
    return xs

for eta in [0.01, 0.3, 1.0, 1.2]:
    traj = gd(0.0, eta, 30)
    plt.plot(traj, label=f"η={eta}")
plt.axhline(3, color="red", ls="--", label="x*=3")
plt.legend(); plt.xlabel("step"); plt.ylabel("x")
plt.title("GD trajectory vs learning rate")
plt.savefig("gd_lr.png", dpi=120)

预期:$\eta=0.01$ 慢慢爬,$\eta=0.3$ 一步到位,$\eta=1.0$ 跨过 3 在两侧震荡,$\eta=1.2$ 直接发散(很快远离 3)。

4.4 取舍 / 直觉

学习率是 ML 调参第一参数。经验值:

• 线性回归:$\eta = 0.01 \sim 0.1$
• 神经网络 SGD:$\eta = 0.1 \sim 0.001$(通常要 + 衰减)
• Adam 几乎”自适应”,$\eta = 1e-3$ 是常见默认

Mini-batch 选 $B$ 时:太小(1-8)震荡大、GPU 利用率低;太大(> 1024)梯度稳但收敛慢;常用 32 / 64 / 128 / 256。

4.5 对比表

5. NumPy 实战:完整线性回归

5.1 最小可运行实现

import numpy as np
import matplotlib.pyplot as plt

# 1) 造数据:y = 4x + 1.5 + 噪声
rng = np.random.default_rng(7)
X = rng.uniform(0, 10, 200).reshape(-1, 1)
y = 4 * X.squeeze() + 1.5 + rng.normal(0, 2.0, 200)

# 2) 增广一列 1 让偏置被吸收进 w
X_aug = np.hstack([X, np.ones((len(X), 1))])

# 3) 正规方程解析解
w_ne = np.linalg.solve(X_aug.T @ X_aug, X_aug.T @ y)
print(f"NE   : w={w_ne[0]:.4f} b={w_ne[1]:.4f}")   # 接近 4.0, 1.5

# 4) GD 数值解
w_gd = np.zeros(2)
eta, N = 0.05, len(y)
losses = []
for step in range(2000):
    y_pred = X_aug @ w_gd
    loss = ((y_pred - y) ** 2).mean()
    losses.append(loss)
    grad = 2 * X_aug.T @ (y_pred - y) / N
    w_gd -= eta * grad

print(f"GD   : w={w_gd[0]:.4f} b={w_gd[1]:.4f}")

# 5) sklearn 对照
from sklearn.linear_model import LinearRegression
w_sk = LinearRegression().fit(X, y)
print(f"SKL  : w={w_sk.coef_[0]:.4f} b={w_sk.intercept_:.4f}")

# 6) 画 loss 曲线
plt.plot(losses)
plt.yscale("log")
plt.xlabel("step"); plt.ylabel("MSE loss")
plt.title("GD loss convergence"); plt.savefig("loss_curve.png", dpi=120)

预期输出:

NE   : w=3.9876 b=1.5234
GD   : w=3.9875 b=1.5235
SKL  : w=3.9876 b=1.5234

三组权重几乎相等,Loss 曲线在对数坐标上线性下降到 4 左右(MSE = σ² = 4)。

5.2 进阶版:向量化 mini-batch + 早停

def linear_regression_minibatch(X, y, eta=0.05, batch_size=32, n_epochs=200, tol=1e-6):
    rng = np.random.default_rng(0)
    Xa = np.hstack([X, np.ones((len(X), 1))])
    w = np.zeros(Xa.shape[1])
    losses, prev = [], np.inf
    for epoch in range(n_epochs):
        idx = rng.permutation(len(Xa))
        for s in range(0, len(Xa), batch_size):
            bi = idx[s:s + batch_size]
            pred = Xa[bi] @ w
            grad = 2 * Xa[bi].T @ (pred - y[bi]) / len(bi)
            w -= eta * grad
        loss = ((Xa @ w - y) ** 2).mean()
        losses.append(loss)
        if abs(prev - loss) < tol: break
        prev = loss
    return w, losses

w_mb, losses_mb = linear_regression_minibatch(X, y, eta=0.05, batch_size=16)
print(f"MB   : w={w_mb[0]:.4f} b={w_mb[1]:.4f}, stopped at epoch {len(losses_mb)}")

预期:权重与全量 GD 接近,但用 16 行小 batch + 早停加速 5-10 倍。

6. vs 其他方法对比

• 手写:今天的主题——理解每一行数学
• sklearn:Day 7 起 Titanic baseline 首选
• PyTorch:Day 15+ 神经网络起点
• XGBoost:Week 2 Kaggle 表格数据首选

7. 常见坑(10 条)

7.1 正规方程遇到 singular matrix

症状:np.linalg.solve 报 LinAlgError: Singular matrix。 原因:$X^\top X$ 不可逆,通常是特征高度相关或特征数 > 样本数。 修法:用 np.linalg.lstsq 求最小二乘;或加 L2 正则(Ridge):$\mathbf{w} = (X^\top X + \lambda I)^{-1} X^\top y$。

7.2 学习率太大 Loss 爆炸

症状:Loss 越训练越大,NaN 出现。 原因:$\eta > 2/\lambda_{\max}(X^\top X)$ 时 GD 不收敛。 修法:梯度裁剪(np.clip(grad, -1, 1))、降低 $\eta$、先归一化特征。

7.3 不归一化特征

症状:Loss 下降极慢,要 10000+ 步才收敛。 原因:不同特征量纲差异巨大,等高线是扁椭圆,GD 在窄缝里 zigzag。 修法:X = (X - X.mean(0)) / X.std(0),或者用 Adam 自动适应。

7.4 把测试集 mean / std 当训练集 mean / std

症状:训练集 $R^2 = 0.95$,测试集 $R^2 = 0.50$。 原因:StandardScaler.fit_transform(X_test) 用测试集自己的均值 / 方差,引入了未来信息。 修法:只对训练集 fit,测试集用 transform。

7.5 用 MSE 评价分类

症状:训练一个分类模型用 MSE Loss,Loss 怎么降 accuracy 都不动。 原因:分类输出是 0/1,MSE 梯度饱和(预测概率在 0.7 时梯度小)。 修法:分类换 CrossEntropy + Softmax(Week 2 Day 8 主题)。

7.6 偏置项漏加

症状:线性回归拟合不过原点附近的数据,残差均值 ≠ 0。 原因:模型假设没偏置 $\hat y = w^\top x$,被强行穿过原点。 修法:给 $X$ 增广一列 1,或 LinearRegression(fit_intercept=True)(默认)。

7.7 矩阵形状算错

症状:X @ w 报 shape 不匹配。 原因:X 是 (N, n),w 应该是 (n,) 或 (n, 1),矩阵乘法不是逐元素。 修法:画 shape 推演图;w[:, None] 显式升维;用 assert 验证 shape。

7.8 用 for 循环算损失

症状:100 万样本 Loss 计算慢 10 分钟。 原因:Python for 循环解释器开销大。 修法:用 numpy 向量化 ((X @ w - y) ** 2).mean(),50-200 倍加速。

7.9 重复 .fit_transform

症状:每次 pipeline.fit_transform(X_train) 都重新计算,grid search 极慢。 原因:cross_val_score 内部对每折重新 fit,这是正常的,但手动循环要避免。 修法:把 fit 和 predict 分离,只在需要时 fit。

7.10 多重共线导致权重爆炸

症状:两个特征高度相关($\rho > 0.99$),学出来的权重一个 +1000 一个 -1000。 原因:最小二乘在病态问题下不稳定。 修法:Ridge(L2 正则)、PCA 降维、删冗余特征。

8. 自检三问

A. 线性回归的损失函数为什么用 MSE(均方误差)而不是 MAE(平均绝对误差)?两种损失对异常值的敏感度和梯度形态有什么本质区别?

要点:高斯噪声假设下,极大似然等价于最小化 MSE,数学上更光滑可导,梯度线性于残差。MAE 等价于 Laplace 噪声 + 极大似然,对异常值更稳健,但损失函数在残差 = 0 处不可导,梯度是 ±1 阶跃,优化器不喜欢。Huber 损失是两者的折中,残差小于阈值用 MSE,大于阈值用 MAE。

B. 学习率 $\eta$ 设得太大时,梯度下降为什么会”发散”而不是”慢一点收敛”?请用损失函数 $L$ 和更新公式 $w \leftarrow w - \eta \cdot \partial L/\partial w$ 画出二阶近似的直觉。

要点:在 $w^*$ 附近二阶泰勒展开 $L(w_{t+1}) \approx L(w_t) + g \Delta w + \frac{1}{2} H (\Delta w)^2$。$\Delta w = -\eta g$,代入得 $L(w_{t+1}) - L(w_t) \approx -\eta g^2 + \frac{1}{2} \eta^2 g^2 H$。要保证下降需 $\eta < 2/H$。超出此值二阶项主导,Loss 反向增大——这就是发散而不是慢收敛。

C. 正规方程 $w = (X^\top X)^{-1} X^\top y$ 在什么情况下会失效?为什么说”特征数 $n$ 很大时”矩阵求逆 $O(n^3)$ 是瓶颈,而梯度下降反而更友好?

要点:失效场景——$X^\top X$ 奇异(特征共线 / 特征数 > 样本数),需要加 L2 正则。计算瓶颈——$n = 10000$ 时,$(X^\top X)$ 求逆是 1e12 次浮点运算,慢到分钟级;GD 每步 $O(N n)$,即使 1000 步也只是 $10^9$ 量级。GD 还支持在线学习 / mini-batch,正规方程必须一次拿到全部数据。

9. 推荐资源(5 类)

视频

• 3Blue1Brown《线性代数的本质》 — 矩阵乘法的几何直观
• StatQuest「Linear Regression, Clearly Explained」(YouTube) — 5 分钟讲清 OLS
• Andrew Ng《Machine Learning Specialization》Course 1 Week 1-2 — 梯度下降 + 学习率

教科书

• 《统计学习导论》(ISLR)第 3 章 — 线性回归 + 诊断
• 《动手学深度学习》(d2l.ai)第 3 章 — 线性回归 + 基础优化
• 《Pattern Recognition and Machine Learning》(Bishop)第 3 章 — 概率视角的线性回归

论文

• “A Stochastic Approximation Method”(Robbins & Monro, 1951) — SGD 原始论文
• “An Overview of Gradient Descent Optimization Algorithms”(Ruder, 2016) — 优化器综述
• “Adam: A Method for Stochastic Optimization”(Kingma & Ba, 2015) — Adam 论文

博客 / 课程

• Kaggle「Intro to Machine Learning」(kaggle.com/learn) — 4 小时跑通 sklearn 线性回归
• Google ML Crash Course「Descending into ML」 — 损失函数 + GD 入门
• Towards Data Science「Linear Regression using Gradient Descent」 — 手把手 Python 实现

代码

• d2l 官方 notebook(github.com/d2l-ai) — 线性回归 + 实战图
• scikit-learn 官方 examples — plot_ols.py 拟合 + 残差可视化
• PyTorch 官方 tutorial — nn.Linear + optim.SGD 入门

10. 本节要点

• 要 1:线性回归假设 $\hat y = w^\top x + b$,用 MSE 损失,$L$ 是关于 $w$ 的凸二次函数,全局唯一极小。
• 要 2:MSE = 高斯噪声 + 极大似然;MAE = Laplace 噪声 + 极大似然;Huber 是折中。
• 要 3:正规方程 $w = (X^\top X)^{-1} X^\top y$ 解析解,小特征数 (< 1000) 1 步到位;大数据用 GD。
• 要 4:GD 更新 $w \leftarrow w - \eta \nabla L$,$\eta$ 太大发散、$H/2$ 是临界值;$\eta$ 太小收敛慢。
• 要 5:Mini-batch GD 是深度学习事实标准,$B=32$ / $64$ / $128$ 是常用选择。
• 要 6:向量化 X @ w 比 for 循环快 50-200 倍,任何能写成矩阵乘法的循环必须替换。
• 要 7:先归一化特征(StandardScaler)再训,先在训练集 fit 再 transform 测试集,是数据泄露的必踩红线。

11. 下一节:Day 07 · 第 1 周复盘 + Kaggle Titanic 入门

主题:第 1 周 7 天知识的串联 + Kaggle Titanic 二分类实战——Day 1-6 的全部工具串成第一条端到端 ML Pipeline。覆盖:

• Week 1 知识串联 Python / NumPy / 线性代数 / 微积分 / 概率统计 / 工具链 / 线性回归 7 天如何拼成”梯度下降 + MSE + 向量化 + sklearn Pipeline”的最小闭环
• Titanic 完整流程 EDA → 缺失值处理 → 特征工程 → 模型 → 提交 CSV → 看 Public LB
• 关键概念回顾 numpy broadcasting / 矩阵乘法 / 梯度 / MSE / 协方差 / pandas / sklearn API
• 实战要点 训练 / 验证 / 测试集划分时机、fit_transform 防泄露、Public LB 与本地 CV 的差距

产出物:第一份 Kaggle 提交文件 submission.csv(Public LB 目标 ≥ 0.77)+ 一段 EDA 总结(性别 / 舱位 / 年龄 vs 生还率可视化)+ 一份 Week 1 知识串联自测表。