手写线性回归：从梯度下降到特征缩放的NumPy实现

1. 为什么放着现成的库不用，非要手写线性回归

1.1 从一次调包翻车说起

刚接触机器学习那会儿，我和大多数人一样，拿到任务第一反应就是from sklearn.linear_model import LinearRegression，三行代码拟合、预测、打分，一气呵成。直到有一次，数据里混进了几个量纲差异极大的特征，模型系数跑出来大得离谱，预测结果完全没法解释。我盯着coef_数组看了半天，才意识到问题出在特征缩放上——可 sklearn 的封装太厚了，我根本不知道它内部到底做了什么、没做什么。

那次之后我就下定决心，把最基础的线性回归用原生 Python 手写一遍。不是为了炫技，而是想彻底搞清楚三件事：梯度到底是怎么下降的、学习率到底怎么影响收敛、特征缩放为什么能救命。这篇内容就是那次折腾的完整记录，从数学推导到 NumPy 实现，再到 Jupyter 里的可视化调试，全部摊开讲。适合已经会调包、但想搞懂底层原理的朋友，也适合正在做课程作业、需要自己实现算法的同学。

1.2 手写算法到底能学到什么

很多人觉得手写线性回归是"重复造轮子"，没意义。我的看法恰恰相反。调包你学到的是 API 用法，手写你学到的是机器学习的思维方式。具体来说，手写一遍你能拿到这些别人拿不走的东西：

• 损失函数的直觉：你会亲眼看到均方误差（MSE）随着迭代次数下降的曲线，而不是只看到一个score数字。
• 梯度下降的脾气：学习率设大了会震荡甚至发散，设小了收敛慢得像蜗牛，这个"手感"只有自己调过才有。
• 矩阵运算的威力：向量化实现和 for 循环实现的性能差距，可能是几十倍甚至上百倍，这个认知会改变你以后写所有代码的习惯。
• 数值稳定性问题：特征量纲不统一时，梯度下降会走"之字形"，收敛极慢，这是教科书上一笔带过、但实操中天天遇到的坑。

提示：手写算法的目的不是替代 sklearn，而是让你在用 sklearn 时，知道每个参数背后发生了什么。生产环境该调包还是调包，但心里要有底。

2. 线性回归的数学骨架，用大白话拆开讲

2.1 从一条直线到矩阵形式

线性回归最朴素的形式，就是找一条直线去拟合数据：y = wx + b。一个特征时好理解，w 是斜率，b 是截距。但真实数据往往有多个特征，比如预测房价，可能同时用到面积、房龄、楼层、距离地铁远近等。这时候就变成：

y = w1*x1 + w2*x2 + ... + wn*xn + b

为了写代码方便，我们通常把 b 吸收进权重向量，给每个样本补一列全 1 的特征，这样公式就统一成y = X @ w。这里的@是矩阵乘法。假设有 m 个样本、n 个特征，那么 X 是 m×n 的矩阵，w 是 n×1 的向量，y 是 m×1 的向量。

这个"补一列 1"的操作叫增广特征矩阵，是手写实现里第一个容易忽略的细节。我见过不少同学手写时忘了补这一列，结果模型永远过不了原点，拟合效果差一大截还找不到原因。

2.2 损失函数为什么选均方误差

衡量预测好坏，我们需要一个"打分"标准，这就是损失函数。线性回归最常用的是均方误差：

J(w) = (1/2m) * Σ(y_pred - y_true)²

前面那个1/2是为了求导时把平方的 2 约掉，纯粹是数学上的便利，不影响最优解的位置。为什么用平方而不是绝对值？因为平方函数处处可导，梯度下降好处理；而且它对大误差的惩罚更重，模型会更努力去修正那些偏离很远的点。

不过这里有个坑：均方误差对异常值非常敏感。一个离群点可能把整条拟合直线拽偏。如果你的数据里有明显异常值，要么先清洗，要么考虑换成平均绝对误差（MAE）。这个取舍在实际项目里经常要面对。

2.3 梯度下降的推导，一步步来

有了损失函数，接下来就是找让它最小的 w。梯度下降的思路很直白：算出当前点的梯度（也就是最陡的上升方向），然后往反方向走一小步，反复迭代。

对 J(w) 求关于 w 的偏导，推导过程如下：

J(w) = (1/2m) * (Xw - y)ᵀ(Xw - y) ∂J/∂w = (1/m) * Xᵀ(Xw - y)

这个结果非常优雅。Xw - y是预测值和真实值的差（残差），Xᵀ乘上残差就得到了梯度。用代码写出来就是一行：

gradient = (1 / m) * X.T @ (X @ w - y)

然后更新权重：

w = w - learning_rate * gradient

就这么简单。但简单不代表好调，学习率的选择、迭代次数、是否收敛，全是学问，后面实操部分细讲。

3. 在 Jupyter 里从零搭一个线性回归

3.1 环境准备与数据构造

我习惯用 Jupyter Notebook 做这类实验，因为可以一段一段跑、随时看中间结果，调试体验比脚本好太多。需要的库就三个：NumPy 做矩阵运算，Matplotlib 画图，Pandas 可选（用来加载真实数据时方便）。

import numpy as np import matplotlib.pyplot as plt # 固定随机种子，保证每次跑结果一致，方便调试 np.random.seed(42) # 构造 100 个样本，2 个特征 m = 100 n_features = 2 # 真实权重，我们最后要看看能不能拟合回来 true_w = np.array([[3.0], [-2.0]]) true_b = 1.5 X_raw = np.random.randn(m, n_features) * 2 y = X_raw @ true_w + true_b + np.random.randn(m, 1) * 0.5 # 加点噪声

这里我故意让两个特征的量纲接近，先跑通流程。等会儿会专门演示量纲差异大时会发生什么。加噪声是为了模拟真实情况，纯线性无噪声的数据拟合起来太"完美"，反而看不出算法的真实表现。

3.2 特征缩放，这一步千万别省

在把数据喂给梯度下降之前，标准化是必须做的。原因很实在：如果特征 A 的取值在 0 到 1 之间，特征 B 的取值在 0 到 10000 之间，那么损失函数的等高线会变成一个又扁又长的椭圆，梯度下降会沿着"之字形"慢慢挪，收敛极慢。

标准化的做法是每个特征减去均值、除以标准差：

def standardize(X): mu = X.mean(axis=0) sigma = X.std(axis=0) return (X - mu) / sigma, mu, sigma X_scaled, mu, sigma = standardize(X_raw)

注意，测试数据必须用训练数据的均值和标准差来缩放，不能各自算各自的。这是新手最容易犯的错误之一，会导致训练和推理时的数据分布不一致，模型表现莫名其妙地差。

3.3 增广矩阵与权重初始化

前面说过要补一列 1，代码里这样写：

def add_bias_column(X): return np.hstack([np.ones((X.shape[0], 1)), X]) X = add_bias_column(X_scaled)

权重初始化一般用全零或小的随机数。线性回归的损失函数是凸的，理论上从哪初始化都能收敛到全局最优，所以全零完全没问题：

w = np.zeros((X.shape[1], 1))

3.4 核心训练循环，逐行注释

把前面的推导组装起来，就是完整的训练函数：

def gradient_descent(X, y, w, learning_rate=0.1, num_iters=1000): m = X.shape[0] loss_history = [] for i in range(num_iters): # 1. 前向计算预测值 y_pred = X @ w # 2. 计算残差 residual = y_pred - y # 3. 计算梯度 gradient = (1 / m) * X.T @ residual # 4. 更新权重 w = w - learning_rate * gradient # 5. 记录损失，用于画收敛曲线 loss = (1 / (2 * m)) * np.sum(residual ** 2) loss_history.append(loss) return w, loss_history

跑起来：

w_final, losses = gradient_descent(X, y, w, learning_rate=0.1, num_iters=1000) print("拟合权重:", w_final.ravel()) print("真实权重:", np.r_[true_b, true_w.ravel()])

实测下来，拟合出的权重和真实值非常接近，误差在小数点后两位。这说明算法本身是对的。

4. 调试与可视化，让收敛过程看得见

4.1 画损失曲线判断收敛

光看最终权重不够，得看损失是怎么降下来的。一条健康的收敛曲线应该是前期快速下降、后期逐渐平缓：

plt.plot(losses) plt.xlabel("迭代次数") plt.ylabel("损失值") plt.title("损失收敛曲线") plt.show()

如果曲线还在明显下降就停了，说明迭代次数不够；如果曲线上下震荡，说明学习率太大；如果曲线几乎是一条平线，说明学习率太小或者数据没标准化。

4.2 学习率对比实验

我强烈建议做一次学习率对比，直观感受它的影响。下面这段代码同时跑三个学习率：

for lr in [0.001, 0.01, 0.1, 1.0]: _, loss_hist = gradient_descent(X, y, np.zeros((X.shape[1], 1)), learning_rate=lr, num_iters=200) plt.plot(loss_hist, label=f"lr={lr}") plt.legend() plt.xlabel("迭代次数") plt.ylabel("损失值") plt.show()

实测结果很有意思：lr=0.001时曲线下降极慢，200 次迭代还没收敛；lr=0.1时又快又稳；lr=1.0时直接震荡甚至发散，损失值飙到很大。这个实验做一遍，你对学习率的理解会比看十篇文章都深。

4.3 量纲差异的"之字形"演示

为了验证标准化的必要性，我特意构造了一组量纲差异巨大的数据：一个特征范围 0-1，另一个范围 0-10000。不标准化直接跑梯度下降，损失曲线会呈现出明显的震荡，收敛需要的迭代次数是标准化后的几十倍。这个对比实验做完，你就再也不会跳过标准化这一步了。

注意：标准化后记得保存均值和标准差，预测新数据时要用同一套参数转换，否则结果全错。

5. 手写实现里那些教科书不讲的坑

5.1 常见问题速查表

5.2 向量化 vs 循环，性能差多少

我做过一个对比：用 for 循环逐样本计算梯度，和用矩阵运算一次性算完，在 10 万样本上性能差了将近 50 倍。原因是 NumPy 底层用了高度优化的线性代数库，而 Python 的循环开销极大。所以手写算法时，能用矩阵运算就别用循环，这是铁律。

5.3 数值稳定性的小技巧

计算损失时，如果样本量很大，np.sum(residual ** 2)可能溢出。稳妥的做法是先除以 m 再求和，或者用np.mean代替。另外，判断收敛时不要用"损失等于 0"，而应该用"损失变化量小于某个阈值"，因为浮点数永远不可能精确等于 0。

6. 从手写版到生产版的过渡思路

6.1 什么时候该用 sklearn

手写版适合学习和理解，但生产环境我依然推荐 sklearn。原因很简单：sklearn 的LinearRegression用的是正规方程（闭式解），不需要调学习率、不需要迭代，在小到中等规模数据上又快又稳。手写梯度下降的优势在于数据量极大、内存装不下时可以做小批量训练，这是另一个话题了。

6.2 把学到的原理反哺调包

手写一遍之后，你再看 sklearn 的参数就会有全新的理解。比如fit_intercept参数，你现在知道它对应的是"要不要补那一列 1"；normalize参数（旧版本）对应的是"要不要标准化"。这种"知道自己在调什么"的感觉，是单纯调包永远给不了的。

6.3 后续可以扩展的方向

这套手写框架很容易扩展。加 L2 正则化，只需要在损失函数里加一项λ * ||w||²，梯度里加λ * w；换成逻辑回归，把预测值过一个 sigmoid 函数，损失换成交叉熵即可。骨架是一样的，理解了线性回归，这些变体都是顺水推舟。

我在实际教学中发现，凡是认真手写过一遍线性回归的同学，后面学神经网络时对反向传播的理解都明显更快，因为梯度下降的直觉已经建立起来了。这个基础打扎实，后面走多远都不慌。