从零手写线性回归：基于NumPy与Jupyter的梯度下降实战

1. 为什么我不调库，偏要手写线性回归

1.1 从一次面试翻车说起

前几年我面过一家做数据服务的公司，技术面聊得挺顺，最后面试官扔过来一道题：“不用sklearn，手写一个线性回归，要求支持多元输入，用梯度下降训练，写在一张纸上。”我当时脑子一空，公式记得，但矩阵维度怎么对齐、学习率怎么设、收敛怎么判断，全卡壳了。那次之后我回去老老实实用NumPy把线性回归从头到尾撸了一遍，跑在Jupyter里，一步步打印损失、画拟合线，才算真正把这块知识焊死在脑子里。

这篇东西就是那次复盘的完整记录。核心关键词就三个：机器学习、Jupyter、Python线性回归，而且重点是“自己写的算法”，不是from sklearn.linear_model import LinearRegression一行了事。它解决的是一个很具体的问题：当你需要真正理解梯度下降在干什么、为什么学习率会影响收敛、正规方程和迭代法各自适合什么场景时，调库是给不了你答案的。适合谁看？刚入门机器学习、被公式劝退过、或者像我一样在面试里栽过跟头的人。哪怕你只会一点点Python和一点点矩阵乘法，跟着走也能跑通。

1.2 手写和调库，到底差在哪

很多人觉得，线性回归这么简单的东西，调库不就行了？我一开始也这么想。但实际做下来你会发现，调库你只拿到了一个结果，手写你拿到的是整个过程的控制权。举个最直接的例子：数据里有量纲差异极大的特征，比如一个特征是面积（几千到几百），另一个是房间数（1到5），调库时你不管它，模型可能也能跑，但手写梯度下降时如果不做标准化，损失会直接炸成NaN，你被迫去理解“为什么特征缩放对梯度下降是刚需”。这种被迫理解，才是手写的价值。

再一个，手写让你能改。比如你想把损失从MSE换成MAE，想加L2正则，想换成随机梯度下降，调库你得翻文档找参数，手写你直接改几行代码。Jupyter在这里的角色就是实验台，一个cell一个cell地跑，改一个参数立刻看损失曲线变化，这种即时反馈对建立直觉太重要了。所以我的定位很明确：这不是一个“替代sklearn”的方案，而是一个教学向、可拆解、可魔改的实现，生产环境该调库还是调库。

2. 开工前的环境与数据准备

2.1 Jupyter环境怎么搭最省心

我试过几种组合，最后稳定在Anaconda + Jupyter Notebook这套上。原因很简单，Anaconda把NumPy、Matplotlib、Pandas这些科学计算的包一次性装齐了，不用一个个pip，省得版本打架。装完之后命令行敲jupyter notebook，浏览器自动弹出来，新建一个.ipynb文件就能开写。

如果你不想装Anaconda这么重的东西，用pip install jupyter numpy matplotlib也行，但要注意Python版本别太老，3.8以上比较稳。我踩过的坑是：某些系统里pip装的NumPy和系统自带的冲突，导致import numpy报奇怪的错，这时候用虚拟环境（python -m venv myenv）隔离一下就好了。

提示：Jupyter里建议第一个cell就写%matplotlib inline，这样画图直接嵌在notebook里，不用弹窗，看损失曲线特别方便。

2.2 造一份自己能看懂的数据

真实数据往往脏且难解释，学习阶段我建议先自己造数据，因为你知道真实规律是什么，能验证算法学出来的对不对。我用的是最经典的“面积-房价”单变量场景，但为了后面演示多元，我多加了一个“房间数”特征。

import numpy as np import matplotlib.pyplot as plt np.random.seed(42) m = 100 # 样本数 area = np.random.uniform(50, 150, m) # 面积 rooms = np.random.randint(1, 5, m) # 房间数 # 真实规律：房价 = 0.8*面积 + 5*房间数 + 20 + 噪声 price = 0.8 * area + 5 * rooms + 20 + np.random.normal(0, 5, m) X = np.column_stack((area, rooms)) y = price.reshape(-1, 1)

这里np.random.seed(42)是为了让每次跑出来的随机数一样，方便复现。column_stack把两个一维特征拼成(100, 2)的矩阵，y转成(100, 1)的列向量，这个维度对齐后面会反复强调，是新手最容易翻车的地方。

2.3 特征标准化，不做这一步后面全白搭

刚才说了，面积是几十到几百，房间数是个位数，量纲差了两个数量级。梯度下降对这种情况极其敏感，因为损失函数在不同方向上的“坡度”差异巨大，学习率稍微大一点，在面积那个方向上就震荡发散了。所以标准化是必须的，我用的是Z-score标准化：

mu = X.mean(axis=0) sigma = X.std(axis=0) X_norm = (X - mu) / sigma

axis=0表示按列（每个特征）算均值和标准差。标准化之后每个特征均值0、方差1，梯度下降就好走多了。记住把mu和sigma存下来，预测新数据时要用同样的参数去标准化，不然结果全错。这个细节调库时是自动帮你做的，手写就得自己操心。

3. 线性回归的数学底子，用人话讲一遍

3.1 假设函数就是个加权求和

线性回归的假设函数，说白了就是“每个特征乘一个权重，加起来再加个偏置”：

$$h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n$$

写成矩阵形式就是 $h = X\theta$，其中$X$是$(m, n+1)$的矩阵（多出来的一列全是1，对应偏置$\theta_0$），$\theta$是$(n+1, 1)$的参数向量。为什么要加一列1？因为偏置项可以看成“特征恒为1”对应的权重，这样统一用矩阵乘法就能表达，代码写起来干净。我一开始老忘了加这列1，结果偏置永远学不出来，拟合线总是过原点，排查半天才发现。

3.2 损失函数：衡量“错得有多离谱”

我们用均方误差（MSE）作为损失：

$$J(\theta) = \frac{1}{2m}\sum_{i=1}^{m}(h_\theta(x^{(i)}) - y^{(i)})^2$$

那个$\frac{1}{2}$是故意的，因为求导时平方的2会和它约掉，公式更简洁。$m$是样本数，除以$m$是为了让损失和样本量无关，方便比较。MSE的几何意义就是所有点到拟合直线的垂直距离平方的平均，越小说明拟合越好。这个函数是凸函数，只有一个全局最小值，所以梯度下降一定能找到最优解，不用担心局部最优，这也是线性回归好教好学的根本原因。

3.3 梯度下降：蒙眼下山，靠坡度找方向

梯度下降的直觉特别简单：你站在山坡上，想走到谷底，但看不见路，只能感受脚下哪个方向最陡，然后朝那个方向迈一小步，重复直到走不动。数学上就是：

$$\theta_j := \theta_j - \alpha \frac{\partial}{\partial \theta_j} J(\theta)$$

$\alpha$是学习率，也就是“步子大小”。对MSE求偏导，结果是：

$$\frac{\partial}{\partial \theta_j} J(\theta) = \frac{1}{m}\sum_{i=1}^{m}(h_\theta(x^{(i)}) - y^{(i)}) \cdot x_j^{(i)}$$

矩阵形式就是 $\frac{1}{m}X^T(X\theta - y)$。这个式子我建议你手推一遍，推完你就明白为什么代码里是X.T.dot(error) / m这个写法了，而不是死记。

4. 手写核心算法：从零到能跑

4.1 给X补上那列1

def add_bias_column(X): m = X.shape[0] return np.column_stack((np.ones(m), X))

np.ones(m)生成全1列，拼在特征矩阵最左边。这样$\theta_0$就对应这一列，后面统一处理，不用单独写偏置逻辑。

4.2 损失函数实现

def compute_loss(X, y, theta): m = len(y) error = X.dot(theta) - y return (1 / (2 * m)) * np.sum(error ** 2)

X.dot(theta)就是$X\theta$，得到预测值，减y得到误差，平方求和除以2m。这个函数我建议单独写出来，训练过程中每隔几十轮打印一次，能直观看到损失在降，心里有底。

4.3 梯度下降主循环

def gradient_descent(X, y, theta, alpha, iterations): m = len(y) loss_history = [] for i in range(iterations): error = X.dot(theta) - y gradient = (1 / m) * X.T.dot(error) theta = theta - alpha * gradient loss = compute_loss(X, y, theta) loss_history.append(loss) if i % 100 == 0: print(f"迭代 {i}, 损失 {loss:.4f}") return theta, loss_history

这段是整个算法的核心。X.T.dot(error)就是$X^T(X\theta - y)$，除以m得到平均梯度，然后theta减去alpha * gradient完成一次更新。loss_history把每轮损失存下来，后面画曲线用。我特意每100轮打印一次，是因为如果学习率设大了，前几轮损失就会飙到天文数字或者NaN，早点发现早点调。

4.4 跑起来看看

X_b = add_bias_column(X_norm) theta = np.zeros((X_b.shape[1], 1)) theta, loss_history = gradient_descent(X_b, y, theta, alpha=0.1, iterations=1000) print("学到的参数：", theta.flatten())

初始theta全设0，这是最常用的初始化，因为线性回归是凸的，从哪起步都能到同一个谷底。alpha=0.1、迭代1000次是我试出来的一个比较稳的组合。跑完你会看到损失从几百降到个位数，参数也接近我们造数据时的真实值（标准化后参数含义变了，但趋势对）。

5. 训练过程的可视化与调参

5.1 损失曲线：判断收敛最直接的工具

plt.plot(loss_history) plt.xlabel("迭代次数") plt.ylabel("损失") plt.title("损失下降曲线") plt.show()

这条曲线如果是一条平滑下降然后趋于水平的线，说明训练正常。如果曲线震荡甚至往上翘，八成是学习率太大；如果下降得特别慢像蜗牛，那是学习率太小。我一般会跑几个不同的alpha对比，比如0.001、0.01、0.1、0.5，把曲线画在一张图上，一眼就能看出哪个合适。这个实验在Jupyter里做特别爽，改个数字重跑一个cell就行。

5.2 学习率怎么选：一个实用的试错法

学习率没有万能值，但有个经验做法：从0.001开始，每次乘10试到1.0，看损失曲线。通常0.01到0.1之间是比较舒服的区间。如果0.001都发散，那说明数据没标准化或者有异常值；如果1.0还能稳定下降，那可以再往上试。我实测下来，标准化之后的数据，0.1基本是个不错的起点。

注意：学习率太大时损失会变成inf或nan，这时候别慌，把alpha调小重跑就行。Jupyter的好处就是可以只重跑出问题的那几个cell，不用从头来。

5.3 迭代次数：看损失曲线决定

迭代次数不是越多越好，够收敛就行。判断标准就是损失曲线变平了，再迭代损失也不怎么降了。我一般先设1000次看曲线，如果还没平就加到5000。但要注意，如果学习率合适，通常几百次就收敛了；如果几千次还在明显下降，那可能是学习率偏小，与其加迭代次数不如调大学习率。

6. 从单变量到多元：维度对齐的坑

6.1 多元只是矩阵变宽了

单变量时$X$是$(m, 2)$（一列1加一列特征），多元时就是$(m, n+1)$。算法代码一行都不用改，因为矩阵乘法天然支持任意列数。这就是矩阵形式的优雅之处。我上面造的数据有两个特征，所以X_b是(100, 3)，theta是(3, 1)，跑出来的三个参数分别对应偏置、面积权重、房间数权重。

6.2 维度不匹配的报错怎么读

新手最常见的报错是shapes (100,3) and (2,1) not aligned，意思是矩阵乘法维度对不上。矩阵乘法要求前一个的列数等于后一个的行数。(100,3)乘(3,1)得到(100,1)，没问题；但如果theta初始化成(2,1)，就会报这个错。排查方法很简单：在每个关键步骤后打印shape，比如print(X_b.shape, theta.shape)，一眼就能看出谁不对。

6.3 正规方程：另一条路，适合小数据

除了梯度下降，线性回归还有解析解，叫正规方程：

$$\theta = (X^T X)^{-1} X^T y$$

代码就一行：

theta_ne = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

它的好处是不用调学习率、不用迭代，一步到位。坏处是要算矩阵求逆，复杂度是$O(n^3)$，特征一多就慢得离谱，而且$X^T X$可能不可逆（特征间高度相关时）。所以我的建议是：特征少于几千个、数据量不大时可以用正规方程验证梯度下降的结果；特征多、数据大时老老实实梯度下降。我一般两个都跑一遍，结果接近就说明梯度下降收敛得不错。

7. 常见问题与排查实录

7.1 损失变成NaN怎么办

这是最高频的问题，原因基本就三个：学习率太大、数据没标准化、数据里有异常值或缺失值。排查顺序我一般是：先确认标准化做了没，再看学习率是不是超过1了，最后检查数据里有没有inf或nan。用np.isnan(X).sum()和np.isinf(X).sum()能快速定位。

7.2 损失不降反升

如果损失曲线往上走，几乎可以肯定是学习率太大，步子迈太大直接跨过了谷底，越走越远。把alpha除以10再试。还有一种可能是梯度算错了，比如忘了除以m，导致梯度被放大m倍，等效于学习率大了m倍。

7.3 参数学出来和真实值差很远

先别怀疑算法，检查两件事：一是标准化参数有没有在预测时一致地应用，二是迭代次数够不够。我遇到过一次，训练时标准化了，验证时忘了用同样的mu和sigma，结果预测全偏。另外，如果特征之间高度相关（比如面积和房间数强相关），参数可能不唯一，但预测值还是对的，这种情况不用纠结单个参数的值。

7.4 一个容易被忽略的坑：y的形状

y必须是(m, 1)的列向量，不能是(m,)的一维数组。因为X.dot(theta)出来是(m,1)，如果y是(m,)，相减时NumPy会广播成(m,m)，结果完全错但不报错，特别隐蔽。我踩过一次，损失看着在降但就是不对，查了半天才发现是y的形状问题。养成习惯：y = y.reshape(-1, 1)。

8. 我个人的几条实操心得

第一，先跑通再优化。别一上来就想着加正则、换优化器，先用最简单的梯度下降把单变量跑通，看到损失下降、拟合线画出来，有了正反馈再往下走。我见过太多人卡在环境配置和维度报错上就放弃了。

第二，Jupyter里多打印中间结果。theta、损失、梯度范数，想看什么就打印什么。手写算法的最大优势就是透明，别浪费这个优势。我习惯在每个函数里加一两行print，跑通之后再删掉。

第三，用正规方程做交叉验证。梯度下降跑完，用正规方程算一遍，两者结果接近就说明你的梯度下降实现是对的。这是最省事的自检方法，比盯着公式看半天强。

第四，标准化参数一定要存。训练完把mu和sigma保存下来（np.save或者直接pickle），预测新数据时用同一套参数处理，否则预测结果毫无意义。这个坑我在实际项目里见别人踩过，线上预测和线下训练不一致，查了一整天才发现是标准化没对齐。

这套代码我后来又扩展过：加了L2正则（就是在损失里加$\lambda\sum\theta_j^2$，梯度里加$2\lambda\theta$），加了随机梯度下降（每次只用一个样本更新），还试过用PolynomialFeatures造多项式特征做非线性拟合。核心的梯度下降框架一直没变，改的都是损失和梯度那几行。这也是我坚持手写的理由——你真正拥有它，才能随意改造它。