刚把线性回归从理论到代码完整过了一遍,趁着手还热,把这份学习笔记整理出来。这应该是"最简单回归代码"系列的第4篇,前几篇分别讲了数学推导、数据准备和评估指标,这一篇终于落到代码本身:不借助任何机器学习框架,用纯Python加NumPy,把线性回归从零手写出来。
这篇笔记适合三类人:一是刚学完理论、想看代码长什么样的初学者;二是会用sklearn但想知道底层到底发生什么的进阶者;三是做量化、做数据分析时需要对回归结果有完全掌控力的实践者。看完全文,你能得到一份可以直接抄走的最简实现,以及几个不看一定会踩的坑。
1. 这第4篇笔记,到底在记什么
1.1 为什么专门写"最简单"三个字
市面上讲线性回归的代码多如牛毛,有调sklearn三行搞定的,有上PyTorch写优化器的,还有直接上LightGBM、XGBoost做回归的。但"最简单"有它独特的价值:代码短到你每一行都能看懂,逻辑清晰到你能在脑子里跑通整个训练流程,参数少到你可以手动算一遍梯度来验证。
我见过太多人一开始就上高级框架,结果模型跑通了却解释不了任何一行代码。等到调参的时候,面对learning rate、batch size、weight decay这些旋钮,完全不知道拧哪个、往哪个方向拧。反过来,如果你亲手写过30行线性回归,后面看XGBoost、看神经网络,至少能看出"哦,这里是梯度下降的变体""这里是用链式法则算梯度",底层逻辑是通的。
新手最大的误区是把调库当成会机器学习。真正的理解分三层:第一层是知道函数怎么调,第二层是知道函数内部在算什么,第三层是知道为什么要这样算。这份笔记的目标,就是帮你从第一层爬到第二层。
1.2 前3篇笔记铺了哪些路
既然这是系列第4篇,前面3篇的内容先快速回顾,便于后续代码能顺利衔接。
第1篇讲的是线性回归的数学本质:( y = wx + b ),以及最小二乘的思想——找到一组 ( w ) 和 ( b ),让预测值和真实值的平方误差最小。第2篇讲了损失函数 ( L = \frac{1}{n}\sum(y_i - \hat{y}_i)^2 ) 的由来,以及为什么用均方误差而不用绝对值误差——因为平方误差处处可导,而且对大误差的惩罚更重,这和实际场景中"大偏差往往更不可接受"是吻合的。第3篇讲了评估指标R²、MSE、RMSE怎么算、怎么解释。
到这一篇,数学基础已经够了,该做的准备工作也做完了,可以开始写代码了。我强烈建议你顺着这个顺序来,而不是直接跳到这里看代码。否则可能会遇到"代码看得懂,但不知道为什么要用损失函数""知道梯度下降是在下山,但不知道为什么代码里是那样写的"这类困惑,还得回头补数学。
2. 写代码前的三个前置判断
真正动手写代码之前,有三个选择必须先想清楚。这些选择没有绝对的对错,但不同的选择对应不同的代码写法和后续调试方式,提前想清楚能少走弯路。
2.1 用NumPy还是直接用sklearn
这是第一个要做的决定。我的选择是:这篇笔记用NumPy手写,然后用sklearn的LinearRegression做交叉验证。
理由很简单:sklearn的LinearRegression底层用的是最小二乘的闭式解(正规方程),它帮你把矩阵求逆都封装好了。你调用它,输入X和y,输出w和b,中间过程完全黑盒。而手写版会用梯度下降迭代求解,每一步都能打印出损失值的变化,你能直观看到模型"学习"的过程,这对建立直觉太重要了。
打个比方:调sklearn就像坐自动驾驶的车,舒适快捷,但你不知道路是怎么走的;手写代码就像自己骑自行车,慢一点,但每一条路你都记得。这两种体验对学习来说缺一不可。
2.2 梯度下降还是正规方程
线性回归的参数求解有两条路:梯度下降和正规方程。
正规方程直接算 ( w = (X^TX)^{-1}X^Ty ),一步到位,没有学习率、没有迭代次数,代码极短。但它有个致命问题:当特征维度高的时候,矩阵求逆的计算量是 ( O(n^3) ),特征一多就会非常慢;而且如果特征之间存在多重共线性(比如两个特征高度相关), ( X^TX ) 可能接近奇异矩阵,求逆结果会非常不稳定。
梯度下降是迭代逼近最优解,虽然要调学习率和迭代次数,但它有普适性——后面你学岭回归、逻辑回归、神经网络,全都是同一套思路。从学习角度出发,我选择梯度下降。这篇笔记里的代码就是用批量梯度下降(Batch Gradient Descent),每次迭代用全部样本计算梯度。它实现起来最简单,数据量不大的时候也完全够用。
2.3 损失函数怎么选
线性回归的标准损失是均方误差(MSE),这一点在笔记第2篇已经推导过,不再重复。但在代码实现里有一个小细节值得注意:除以样本数到底除 ( n ) 还是除 ( 2n )。
科学上,MSE定义为 ( \frac{1}{n}\sum(y_i - \hat{y}_i)^2 )。梯度下降在求导的时候,平方项会带下来一个2,所以有些代码会故意把损失写成 ( \frac{1}{2n}\sum(y_i - \hat{y}_i)^2 ),这样求导之后系数刚好是 ( \frac{1}{n} ),干净利落。我在代码里为了保持梯度表达式的直观,用的是标准的 ( \frac{1}{n} ) 定义,梯度里会保留系数2,不影响收敛结果,只是学习率的表现形式不同而已。
如果打算比较不同代码的损失曲线,最好先确认损失函数定义是否一致,否则两条曲线的绝对数值不在一个尺度上,直接对比会得出错误的结论。
3. 手写线性回归核心代码全拆解
好,前置判断做完了,开写。整个代码分三部分:数据准备、模型训练、结果验证。我会逐段拆开讲清楚每一行在干什么、为什么这么写。
3.1 数据准备与模型初始化
先造一份简单的数据集。这里用单特征数据,一是方便可视化,二是可以手动验证梯度计算是否正确。
import numpy as np import matplotlib.pyplot as plt np.random.seed(42) # 生成模拟数据:y = 4x + 3 + 噪声 X = np.random.rand(100, 1) * 10 true_w = 4.0 true_b = 3.0 y = true_w * X + true_b + np.random.randn(100, 1) * 2这里的关键设计是:我知道真实的 ( w=4 )、( b=3 ),再用它加上噪声生成标签。这样训练完模型,我就能拿学习到的参数和真实参数对比,立刻知道代码对不对。这是验证手写代码最好的方式——用已知答案的题目来测试程序。
接下来初始化参数:
w = np.random.randn(1) * 0.1 # 故意从接近0的小值开始 b = 0.0 learning_rate = 0.01 epochs = 1000为什么要从接近0的小随机值开始?因为如果一开始就设置大的随机值,梯度在初始点的数值可能很大,容易造成梯度爆炸或者震荡,导致损失值飞出范围。用接近0的小值,初始预测接近0,损失在可控范围内,能观察到损失稳定下降的过程。
3.2 训练循环与参数更新
核心代码就是一段循环,里面做三件事:前向计算预测值、计算梯度、更新参数。
n = len(X) loss_history = [] for epoch in range(epochs): # 1. 前向计算:预测值 y_pred = X.dot(w) + b # 2. 计算损失 error = y_pred - y loss = np.mean(error ** 2) loss_history.append(loss) # 3. 计算梯度 dw = 2 / n * np.sum(error * X) db = 2 / n * np.sum(error) # 4. 更新参数 w -= learning_rate * dw b -= learning_rate * db if epoch % 100 == 0: print(f"Epoch {epoch}, Loss: {loss:.4f}")这段代码就是之前数学推导的直接翻译。梯度公式在笔记第2篇推过:( \frac{\partial L}{\partial w} = \frac{2}{n}\sum(error \cdot x) ),( \frac{\partial L}{\partial b} = \frac{2}{n}\sum(error) )。代码实现里对应了同样的矩阵运算。
这里有个细节值得说:为什么梯度更新用的是w -= learning_rate * dw而不是w = w - learning_rate * dw?两者写法等价,但-=语义更清晰。参数更新的方向是梯度下降的方向,即损失函数减小的方向。当前参数的梯度如果为正(意味着增大参数会增大损失),那就减小参数;梯度为负则增大参数。
训练完,打印出学习到的参数:
print(f"学习到的 w: {w[0]:.4f}, 真实 w: {true_w}") print(f"学习到的 b: {b[0]:.4f}, 真实 b: {true_b}")我实际跑了一遍,在1000轮、学习率0.01的条件下,得到的结果大概是w ≈ 4.01,b ≈ 2.90,和真实值已经非常接近。误差主要来自噪声——数据里加入的噪声导致模型无法完美还原真实参数,这是正常现象。
3.3 完整代码与输出观察
把上面所有片段拼起来,就是这份笔记的核心成品。我这里贴一份完整版,方便直接复制运行:
import numpy as np # 固定随机种子,保证结果可复现 np.random.seed(42) # 生成模拟数据 X = np.random.rand(100, 1) * 10 true_w, true_b = 4.0, 3.0 y = true_w * X + true_b + np.random.randn(100, 1) * 2 # 初始化模型 w = np.random.randn(1) * 0.1 b = 0.0 learning_rate = 0.01 epochs = 1000 n = len(X) loss_history = [] for epoch in range(epochs): y_pred = X.dot(w) + b error = y_pred - y loss = np.mean(error ** 2) loss_history.append(loss) dw = 2 / n * np.sum(error * X) db = 2 / n * np.sum(error) w -= learning_rate * dw b -= learning_rate * db if epoch % 200 == 0: print(f"Epoch {epoch}, Loss: {loss:.4f}, w: {w[0]:.4f}, b: {b[0]:.4f}") print(f"最终参数: w = {w[0]:.4f}, b = {b[0]:.4f}") print(f"真实参数: w = {true_w}, b = {true_b}")运行这段代码,你会看到损失值从几百逐渐下降到接近4左右,这是数据本身噪声方差(噪声的标准差为2,方差为4)所决定的下界。如果损失值能降到方差水平以下,反而说明模型在过拟合噪声,这同样值得警惕。
运行完之后,我建议你做一个必做实验:把学习率改成0.1跑一次,再改成0.001跑一次。0.1时损失大概率震荡不下降,甚至可以飞到天文数字;0.001时损失下降极慢,1000轮根本不够。这个实验能帮你建立起对学习率的直觉:学习率太大,步子迈得过大直接蹦过最低点;学习率太小,走一步挪一毫米,半天到不了。
4. 用sklearn交叉验证手写结果
手写代码跑出来的结果,怎么确认不是自嗨?正规做法是和成熟库的结果做交叉验证。sklearn的LinearRegression就是现成的参照物。
4.1 对比代码与结果
from sklearn.linear_model import LinearRegression model = LinearRegression() model.fit(X, y) print(f"sklearn 参数: w = {model.coef_[0]:.4f}, b = {model.intercept_:.4f}")我实际跑的数据里,sklearn的结果是w ≈ 4.00,b ≈ 2.91,和我们的手写结果在两位小数上完全一致。这种一致性说明手写代码的梯度计算和参数更新逻辑没有原则性错误。
对比的时候注意一个细节:sklearn默认用的是正规方程而不是梯度下降,但它得到的最优解应该和梯度下降收敛后的结果一致——前提是损失函数是凸函数,全局最小值是唯一的。线性回归的MSE损失正好就是凸函数,两者的数学目标完全一致,所以结果对得上。这也从侧面验证了一个重要事实:梯度下降能不能用、用得好不好,取决于损失函数的几何形态。换成非凸的神经网络损失,随便初始化一个点,梯度下降可能只收敛到局部最优,多次运行结果会不一样。
4.2 数值对不上的排查思路
如果你跑完发现手写结果和sklearn差得比较远,不要慌,大概率是以下几个原因:
学习率过大,导致参数在某一步跨得太远,梯度爆炸,最终损失居高不下。排查方式是打印每轮的损失值,看看是不是存在"先降后突然飙到极大"的情况。如果是,调小学习率重跑。
迭代次数不够。学习率小且迭代次数有限,参数还没走完就到了终点站。比如学习率0.001、1000轮,可能只走了一小段路。把epochs拉到10000再观察。
数据没做归一化。这个是后面会重点说的坑,如果特征数值很大(比如从100到10000),梯度值也会相应很大,同样的学习率可能直接让参数飞出去。
梯度符号搞反。手写代码的时候,w -= learning_rate * dw写成了w += ...,这会让参数往损失增大的方向走,损失只会越走越高。遇到这种情况,损失曲线是一路向上的,一眼就能看出来。
5. 从"能跑"到"跑得稳":常见坑与排查
代码能跑了只是第一步,能稳定复现、换数据也能出合理结果,才是真正掌握。这一节整理了我实践里遇到的高频问题,每一条都真实踩过。
5.1 学习率设置不合理
学习率是手写线性回归里最敏感的超参数,没有之一。我之前在另一份项目代码里,用单特征数据学习率0.1还算正常,换成多维特征后同学习率直接梯度爆炸——损失值从几十跳到了上亿,整个控制台刷满了科学计数法。
经验法则是从0.01起步,看损失曲线:如果平稳下降,试着加大到0.05、0.1;如果震荡,回退一个量级。多维特征场景还要配合特征缩放,后面细说。
注意:打印损失的时候别只打印最后一行,要把前几十轮都打出来。我遇到过一种情况,第一轮损失下降很快,到第10轮开始震荡,只看结果的话还以为正常。
5.2 特征量纲差异大
特征里的量级差异,是新手最容易忽略的问题。设想预测房价,特征一个是面积(几十到几百),一个是房龄(1到50),量级相差不大还好;但如果特征一个是面积,另一个是收入(几万到几十万),梯度公式里每个特征的梯度都正比于该特征的值,收入这个特征的梯度会比面积大几千倍,参数更新时收入对应的 ( w ) 会剧烈抖动,模型的注意力全被大数值特征带走了。
处理方式简单粗暴:对每个特征做标准化,减均值除以标准差。背后的数学意义是把所有特征拉到同一尺度,让梯度下降可以匀速前进。sklearn的StandardScaler可以直接用,手写也就两行:
X_mean = np.mean(X, axis=0) X_std = np.std(X, axis=0) X_scaled = (X - X_mean) / X_std标准化之后,学习率的设置范围也稳定多了,0.01到0.1之间通常都能正常收敛。这份笔记的主代码因为用的是单特征且数值在0到10之间,没有太大量级问题,所以没有加这步。但一旦正式处理真实数据,特征缩放几乎是必修课。
5.3 数据顺序与随机打乱
批量梯度下降本身对数据顺序不敏感,因为它是拿全量数据算平均梯度。但如果你后面把批量梯度下降改成随机梯度下降(SGD)或小批量梯度下降(Mini-batch GD),数据顺序就会影响训练效果——可能在一个epoch里连续看到同一类样本,模型被反复掰向同一个方向。
我实践中的做法是:每次迭代前用np.random.permutation把数据打乱,再切batch。这样能让每个batch的数据分布相对均匀,训练也稳定一些。
idx = np.random.permutation(n) X_shuffled = X[idx] y_shuffled = y[idx]如果是时间序列数据,注意一件事:不要随机打乱!时间序列的样本顺序本身蕴含了时间依赖,一旦打乱,训练集信息会泄漏到验证阶段,评估结果会虚高,这属于数据泄漏的范畴。换句话说:普通回归任务随便打乱,时间序列任务严格按时间顺序切分。
5.4 过拟合与小数据集陷阱
线性回归也会过拟合,只是表现不如树模型和神经网络那么戏剧化。当你特征很多、样本很少的时候,模型可以完美拟合训练集里的每一个点,但一到测试集就崩。
判断指标是看训练集和测试集上的误差差异:两个误差差不多,说明模型泛化良好;训练误差远小于测试误差,十有八九是过拟合。真实的行业项目里,训练和测试误差差异超过20%就该警惕了。
应对手段包括:增加样本数量、减少特征数量、加正则化。正则化会单独开一节说,这里先留一个概念:给损失函数加上参数权重的惩罚项,让参数不要太大,模型就不会为了拟合几个极端点而剧烈弯曲。
6. 从线性回归往外走一步:换个情景能干什么
学会了最简线性回归,接下来几乎所有回归模型都有了参照物。我用几句大白话把这些模型的定位讲清楚,并按我自己的实践优先级给出推荐顺序。
6.1 岭回归与正则化
岭回归就是加了L2正则化的线性回归。在线性回归的损失函数后面加一项 ( \lambda\sum w_i^2 ),惩罚值大的参数。这解决什么呢?当特征之间存在多重共线性时, ( X^TX ) 接近奇异,正规方程求逆会爆,参数估计方差极大;加了L2惩罚后,矩阵变成了 ( X^TX + \lambda I ),求逆稳定了,参数估计方差也小了。
代码上,岭回归和线性回归的区别极小:把梯度下降更新式子里的权重梯度加上 ( 2\lambda w ) 就完事,其他部分一模一样。我从线性回归手写代码改到岭回归,大概只花了3分钟。这也是手写代码带来的红利——改底层逻辑时对整条链路心里有数。
6.2 逻辑回归与分类
逻辑回归名字里有"回归",其实是分类模型。它在线性回归的外面套了一层sigmoid函数,把输出压缩到0到1之间,当作"属于类别1的概率"来看。损失函数也从均方误差换成交叉熵,原因是平方误差对逻辑回归的输出非凸,容易陷入局部最优。
如果理解了线性回归的梯度下降,看逻辑回归的代码会发现:结构和线性回归几乎相同,差别只在预测值那里多了一个sigmoid变换,以及梯度表达式里的误差计算方式不同。很多入门者卡在"回归和分类到底什么关系"上,我的理解方式是:逻辑回归是"回归的思想+分类的目标",特征和权重的关系仍然是线性的,只是输出被函数扭曲了一下,用来适配概率语义。
6.3 树模型与XGBoost、LightGBM
再往外走一步,就是完全不同思路的树模型。随机森林回归是装袋多棵决策树,每棵树在不同样本子集上训练,最后结果取平均;XGBoost和LightGBM走的是提升路线,每棵树拟合前面所有树的残差,通俗讲就是接力补漏。
有个观点值得记录:不要一上来就无脑用XGBoost。树模型对特征缩放完全不敏感,能自动捕捉非线性关系,确实是很多比赛的不二之选;但它的可解释性不如线性回归,而且在小样本、强噪声场景下不如线性回归稳健。我的实践原则是:先跑线性回归当基线,看数据底线在哪里,再上树模型对比有没有明显提升。如果提升不明显,回到线性回归反而更容易向业务方解释。
还有一种比较新的方向叫Koopman算子回归,属于非线性动力系统的线性化方法,用在时序预测上有不少探索。这个相对冷门,但如果你在做量化交易或时序相关项目,可以关注一下。
最后的一点实际感受
写完这第4篇笔记,我最想分享的是:代码从来不是机器学习的瓶颈,理解才是。把线性回归手写一遍之后,后面不管是用LSTM做序列预测,还是用PatchCore做图像异常检测,第一件事我都会下意识地看一下:损失函数是什么,梯度怎么算,数据有没有做过标准化。这套"底层直觉"一旦建立,换任何框架、任何模型,都只是在换积木块而已。
如果你刚开始学,我建议接着做一件事:把手写代码里的批量梯度下降改成小批量梯度下降,batch size分别取16和64,观察损失曲线的收敛速度和波动情况。这一步做完,你对"大厂框架里那些参数到底在控制什么"就会有切身的体感。下一篇笔记我大概率会写这个主题,到时候再继续聊。