尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

手写线性回归:NumPy从零实现解析解与梯度下降

手写线性回归:NumPy从零实现解析解与梯度下降 动手实现一个线性回归模型听起来像是机器学习教材里的课后作业但真正自己用 NumPy 把每个公式敲出来之后你才会发现之前调 sklearn 时那些理所当然的东西全都有了全新的理解。这篇文章就是一次完整的代码实战记录目标只有一个不借助任何机器学习框架从零手动实现经典的 Linear Regression 模型并把它用在一个可验证的模拟数据集上。这不是一篇调包侠教程也不是贴一段公式就完事的理论科普。我会把整个推导、编码、训练、评估、排错的过程拆开揉碎讲清楚每一步为什么要这么做、背后的数学直觉是什么、实际跑代码时又会踩到哪些坑。适合刚入门机器学习、想夯实基础的开发者也适合那些用惯了 sklearn 但对底层原理总觉得差点意思的朋友。文章里的所有代码都基于 NumPy可以直接复制运行边看边敲效果最好。1. 手动实现之前先搞懂线性回归在解决什么问题1.1 从直觉到数学模型到底在学什么线性回归解决的是回归问题也就是预测一个连续的数值。最经典的场景就是房价预测给出一套房子的面积、房龄、地段评分等特征模型输出一个价格。在最简单的一元情形下它就是在找一条直线让这条直线尽可能贴合训练数据点。这条直线可以写成y w * x b其中 w 是斜率权重b 是截距偏置。放到多维场景里输入特征有 n 个式子就变成y w1x1 w2x2 ... wn*xn b写成矩阵形式会更紧凑y X * theta这里 X 是形状为 (m, n) 的特征矩阵m 个样本n 个特征theta 是形状为 (n1,) 的参数向量包含 w 们和 b。为了让这个式子成立我们通常在 X 的最左边加一列全 1作为偏置项的哑特征这样 b 就可以当作一个普通的参数 w0 来统一处理。所以说白了线性回归要做的就一件事找到一组参数 theta让模型预测出来的 y 和真实标签 y 之间的差距尽可能小。1.2 为什么手动实现比调库更有价值很多初学者会有疑问sklearn 里LinearRegression一行代码就能搞定为什么要自己造轮子我的看法是调库解决的是用手动实现解决的是懂这两者之间的差距往往决定了你能走多远。手动实现至少有三个不可替代的价值。第一你会真正理解损失函数、梯度、学习率这些概念是活的而不是教科书里的名词。当你看到 loss 曲线在迭代过程中一点点下降当你发现学习率调大一点就发散、调小一点就慢得让人抓狂你才真正明白这些超参数背后的意义。第二你会具备独立调试模型的能力。工作中遇到模型不收敛、结果全是同一个值、梯度爆炸这类问题如果只会在 sklearn 里换参数基本无从下手但如果你亲手写过梯度下降你会自然地从梯度计算、特征尺度、学习率这几个方向去排查。第三这是你后续学习一切复杂模型的基石。逻辑回归、Softmax 回归、神经网络它们的核心训练过程仍然是前向计算 - 计算损失 - 反向求梯度 - 更新参数线性回归就是这套流程最简单的样本。把它啃透了后面的路会顺畅很多。1.3 评估模型好坏的标尺MSE 与 R²训练模型之前必须先定义好的标准。线性回归最常用的损失函数是均方误差Mean Squared Error, MSEMSE (1/m) * sum((y_pred - y_true)^2)之所以用平方而不是绝对值是因为平方误差处处可导而且对大误差的惩罚更大梯度下降用起来非常方便。为了让求导后系数好看很多教材会写成 J (1/(2m)) * sum((y_pred - y_true)^2)这个 1/2 在求导时会被抵消本质上不影响最优解。另一个常用指标是 R²决定系数它衡量的是模型解释了数据中多少方差R² 1 - SS_res / SS_totSS_res 是残差平方和SS_tot 是总平方和预测值等于均值时的误差。R² 的最大值是 1越接近 1 说明模型拟合效果越好如果 R² 为 0说明模型和直接用均值预测没什么区别如果 R² 为负说明模型比用均值预测还差基本是搞错了方向。我通常会两个指标一起看MSE 反映绝对误差量级R² 反映相对拟合程度组合起来才能对模型效果有个全面判断。2. 两条路线解析解与梯度下降2.1 最小二乘解析解一步到位的正规方程线性回归有一个非常漂亮的数学性质它的损失函数是凸函数也就是说存在全局唯一的最优解而且这个最优解可以直接用矩阵运算求出来不需要迭代。方法就是令损失函数对 theta 的导数等于 0解这个方程。推导过程我这里直接给结论最终得到的就是正规方程theta (X^T * X)^(-1) * X^T * y这个公式的意义很直观X^T * X 是特征间的协方差矩阵加了个偏置列X^T * y 是特征与标签的相关性两者一组合直接就能算出最优参数。整个过程没有任何超参数一次矩阵运算搞定。但解析解有个明显的局限当特征数量 n 非常大时X^T * X 是一个 n×n 的矩阵求逆的计算复杂度是 O(n³)特征上万时基本跑不动。另外如果 X^T * X 不可逆比如特征之间存在多重共线性或者样本数少于特征数求逆就会出问题虽然可以用伪逆兜底但解就不再那么稳定了。2.2 梯度下降用迭代逼近最优解梯度下降的思路完全不一样不奢求一步到位而是从一个初始参数出发沿着损失函数下降最快的方向也就是负梯度方向一步步走直到走到最低点附近。直观理解就是下山你站在山顶不知道最低点在哪但你知道脚下哪个方向最陡朝最陡的方向迈一步然后重新判断方向再迈一步如此反复最终总能到达谷底。对于 MSE 损失参数更新的公式是theta theta - learning_rate * (1/m) * X^T * (X * theta - y)其中(X * theta - y)是预测值与真实值的误差向量X^T * 误差就是梯度方向。learning_rate学习率控制每一步迈多大太大容易一步跨过头在山谷两边震荡甚至发散太小则收敛极慢半天到不了底。梯度下降的优势在于它只需要算矩阵乘法和减法每次迭代的计算量和样本数、特征数线性相关即使特征非常多也能跑。而且它是一个普适框架换一个损失函数、换一个模型结构训练流程几乎不用改。这也是为什么深度学习训练全都基于梯度下降的变体。2.3 两种方案怎么选我的实践建议说一个我自己的经验如果是纯粹的教学演示、特征数小于一万、样本量适中直接用解析解因为它精确、无超参数、代码还短非常适合用来验证梯度下降实现得对不对。但如果是真实业务场景特征成千上万、样本以万计或者后续想在这个基础上扩展成 Ridge、Lasso、逻辑回归那就直接上梯度下降把训练流程搭好后面所有模型都能复用这套骨架。对比维度解析解正规方程梯度下降计算方式一次矩阵运算迭代逼近超参数无学习率、迭代次数特征量级适合特征较少10000适合特征较多计算复杂度O(n³) 求逆O(iter * m * n)扩展性差换模型要重推好通用训练框架数值稳定性特征相关时可能不稳定依赖学习率和归一化3. 代码实战从零实现 Linear Regression3.1 准备环境与模拟数据动手之前先把环境备好。只需要 NumPy 和 Matplotlib前者做矩阵运算后者画 loss 曲线和数据点。如果你在 Jupyter Notebook 里跑直接两个 import 就行import numpy as np import matplotlib.pyplot as plt接下来生成一份我们自己知道答案的模拟数据。这样做的最大好处是模型训练完之后我们可以拿学习到的参数和真实参数直接对比一眼就能看出模型学得对不对。这里我设置真实权重 w2.5、真实偏置 b1.2并在标签中加入高斯噪声模拟真实世界的不完美np.random.seed(42) n_samples 100 X np.random.rand(n_samples, 1) * 10 # 特征范围 0~10 true_w 2.5 true_b 1.2 y true_w * X true_b np.random.randn(n_samples, 1) * 1.5这里生成 100 个样本特征 X 在 0 到 10 之间均匀分布标签 y 由 2.5x 1.2 加上标准差为 1.5 的高斯噪声构成。先画一张散点图确认数据分布plt.scatter(X, y, alpha0.7) plt.xlabel(X) plt.ylabel(y) plt.title(Synthetic Data) plt.show()从图上你应该能看到一条明显的线性趋势但点并不完全落在直线上这正是噪声带来的效果。后续的训练就是去伪存真把藏在散点背后的那条直线找回来。3.2 实现解析解版本解析解版本代码非常短核心就是把正规方程翻译成矩阵运算。注意第一步一定要给 X 加上一列全 1否则偏置项没地方放def fit_normal_equation(X, y): X_b np.c_[np.ones((X.shape[0], 1)), X] theta np.linalg.inv(X_b.T X_b) X_b.T y return theta theta_ne fit_normal_equation(X, y) print(f解析解结果: w{theta_ne[1][0]:.4f}, b{theta_ne[0][0]:.4f})代码里np.c_用于按列拼接是矩阵乘法np.linalg.inv求逆。跑完之后你应该能看到输出大概在 w2.48、b1.2 附近和真实的 w2.5、b1.2 很接近说明只要模型没学错方向解析解能非常精准地逼近真实参数。有个细节值得注意这里theta的形状是 (2, 1) 而不是 (2,)原因是 y 本身是列向量。如果你在别处看到theta是 (2,) 的写法只是向量形状不同数学上完全等价。3.3 实现梯度下降版本梯度下降代码稍微多一点但逻辑非常清晰。我习惯拆成两个函数一个算梯度一个跑训练循环。算梯度的函数如下def compute_gradient(X_b, y, theta): m len(y) error X_b theta - y gradient (1 / m) * X_b.T error return gradient训练循环里每轮迭代做四件事用当前参数预测、算误差、算梯度、更新参数。同时记录每一轮的损失值方便事后画 loss 曲线检查收敛情况def gradient_descent(X, y, learning_rate0.01, n_epochs1000): X_b np.c_[np.ones((X.shape[0], 1)), X] theta np.random.randn(X_b.shape[1], 1) * 0.1 loss_history [] m len(y) for epoch in range(n_epochs): gradient compute_gradient(X_b, y, theta) theta theta - learning_rate * gradient loss np.mean((X_b theta - y) ** 2) loss_history.append(loss) return theta, loss_history theta_gd, loss_history gradient_descent(X, y, learning_rate0.01, n_epochs500) print(f梯度下降结果: w{theta_gd[1][0]:.4f}, b{theta_gd[0][0]:.4f})跑完之后你大概率会发现一个问题梯度下降学出来的 w 和 b 跟解析解差了不少而且 loss 曲线下降得很慢500 轮迭代远没有收敛。原因就是我前面埋的坑——这里的特征 X 在 0 到 10 之间尺度偏大导致梯度在某个方向上特别陡峭学习率 0.01 在这种尺度下显得太小了。这正是我们要在下一节重点解决的问题。先做个对比验证加深印象把学习率直接调到 0.1 再跑一次看看 loss 曲线是更陡还是直接发散。我建议你亲手试一下震荡和缓慢下降两种现象都见过之后你才对学习率有真正的体感。3.4 用 R² 评估训练效果参数学出来之后需要用指标量化效果。我写了一个简单的 R² 计算函数def r2_score(y_true, y_pred): ss_res np.sum((y_true - y_pred) ** 2) ss_tot np.sum((y_true - np.mean(y_true)) ** 2) return 1 - ss_res / ss_tot y_pred_ne np.c_[np.ones((X.shape[0], 1)), X] theta_ne pred_gd np.c_[np.ones((X.shape[0], 1)), X] theta_gd print(f解析解 R²: {r2_score(y, y_pred_ne):.4f}) print(f梯度下降 R²: {r2_score(y, pred_gd):.4f})如果一切正常解析解的 R² 应该在 0.75 左右噪声占比决定了上限不可能无限接近 1而梯度下降在这个学习率下会明显偏低。这个对比就是后续优化的动力R² 不达预期别急着换模型先检查训练过程是不是真的收敛了。4. 训练过程中的关键细节与避坑要点4.1 数据归一化影响收敛的关键一步梯度下降最经典的一个坑就是特征尺度不一致。为什么特征尺度会影响收敛用一个二维的损失函数地形图来解释如果两个特征的取值范围差异很大比如一个在 0~1、另一个在 0~10000那么损失函数形成的等高线会是非常扁长的椭圆形梯度方向并不直接指向圆心而是以之字形缓慢前进。每一步都走得歪歪扭扭不仅收敛慢还容易在学习率稍大时发生震荡。解决办法是特征缩放最常用的是标准化把每个特征变成均值为 0、方差为 1 的分布def standardize(X): mean np.mean(X, axis0) std np.std(X, axis0) X_scaled (X - mean) / std return X_scaled, mean, std注意标准化之后模型的解释性会变差——学出来的权重不再是x 每增加 1 个单位y 增加 w而是x 每增加 1 个标准差y 增加 w。如果你只需要预测准确直接标准化没有任何问题但如果你要解释模型系数就需要把标准化后的权重换算回原始尺度或者干脆用解析解。在这个模拟数据例子里X 已经落在 0~10还算能忍但真实数据里年龄和收入、面积和房价这种跨数量级的特征组合非常常见不做归一化的梯度下降会让人非常痛苦。4.2 学习率怎么选从发散到收敛的实战调参学习率是梯度下降里最敏感的超参数。我踩过无数坑总结了一套简单可执行的调试策略先跑一个很短的训练比如 50 轮画出 loss 曲线观察它是上升、下降还是震荡。如果 loss 爆炸式上升说明学习率太大果断调小通常从 0.01 开始尝试不行就 0.001、0.0001。如果 loss 稳定下降但速度极慢说明学习率偏小可以适当放大。比较理想的情况是loss 前期快速下降中后期进入一个平滑的收敛区间。如果后期一直在小范围抖动要么是学习率偏大要么是数据噪声本身就很大可以接受。我习惯用 10 的整数幂做网格搜索比如 [0.1, 0.01, 0.001, 0.0001]先粗后细省时省力。4.3 Batch Gradient Descent 与 Mini-batch 的差别目前为止我们实现的是最朴素的批量梯度下降Batch Gradient Descent简称 BGD每次更新都要用全部样本计算梯度。它的优点是梯度方向稳定、数学上收敛性最好缺点是当样本量达到百万级别时一次迭代就要做一次超大矩阵运算太慢了。另一种极端是随机梯度下降Stochastic Gradient DescentSGD每次只随机挑一个样本算梯度。它的优点是更新极快、还能跳出局部极小点缺点是梯度噪声太大loss 曲线会像心电图一样抖个不停需要精心设计学习率衰减才能稳定收敛。工程上用得最多的是折中方案——Mini-batch Gradient Descent每次随机取一小批样本比如 32、64、128 个计算梯度兼顾效率和稳定性def minibatch_gradient_descent(X, y, learning_rate0.01, n_epochs200, batch_size32): X_b np.c_[np.ones((X.shape[0], 1)), X] theta np.random.randn(X_b.shape[1], 1) * 0.1 m len(y) loss_history [] for epoch in range(n_epochs): indices np.random.permutation(m) X_shuffled X_b[indices] y_shuffled y[indices] for i in range(0, m, batch_size): X_batch X_shuffled[i:i batch_size] y_batch y_shuffled[i:i batch_size] grad (1 / len(y_batch)) * X_batch.T (X_batch theta - y_batch) theta theta - learning_rate * grad loss np.mean((X_b theta - y) ** 2) loss_history.append(loss) return theta, loss_history注意一个关键细节每次 epoch 之前一定要做np.random.permutation打乱数据否则如果原始数据按某个规律排序模型可能会学到虚假的样本顺序关联。5. 常见问题排查与调试实录5.1 损失不降反升问题出在哪我在训练神经网络时经常遇到 loss 越跑越高的情况线性回归里虽然少见但一旦出现基本就两三个原因。第一是学习率过大更新步长超过了下山所需的距离直接跳到更高的山坡上第二是数据没有归一化某些特征尺度特别大导致梯度方向严重偏移第三是代码 bug比如更新公式里的正负号搞反了、梯度的m忘记除、或者 X_b 拼接的位置不对。排查方法也很简单把学习率降到极小值比如 1e-7如果 loss 还是上升那就不是学习率的问题而是代码本身有 bug如果 loss 开始下降但慢如蜗牛说明你离正确的学习率区间不远了。5.2 收敛到全 0 的权重是怎么回事这是比较隐蔽的一种情况模型训练完后预测值全是一个常数或者权重全是 0。通常原因是特征矩阵里混进了一列零或者某两个特征完全共线导致梯度在某个方向上是 0参数更新永远停在那里。我在一次用真实数据建模房价时也遇到过这个问题后来发现有两个特征本质上是同一个信息的不同单位表达相关性接近 1解析解直接不稳。解决方法是去掉冗余特征或者加 L2 正则化这就是 Ridge 回归的用武之地。5.3 训练集效果很好测试集崩了欠拟合还是过拟合线性回归本身模型复杂度很低和深度网络那种严重的过拟合不太一样但依然存在两种极端情况。欠拟合表现为训练集和测试集的 R² 都很低说明线性假设本身就不成立需要考虑加特征、加多项式项或者换更复杂的模型。过拟合则表现为训练集 R² 接近 1、测试集明显偏低如果你的样本量特别少、特征数接近样本数线性回归也照样会过拟合解决办法是加正则化、提前停止训练、或者增加样本量。5.4 和 sklearn 结果对不上别慌很多人在对比自己的实现和 sklearn 时发现权重存在微小差异第一反应是我是不是哪里写错了。其实大多数情况下不是的。sklearn 的LinearRegression默认用的是基于scipy.linalg.lstsq的最小二乘解法它在计算过程中会做矩阵分解比如 SVD数值稳定性更好此外正则化参数、是否对数据做中心化这些细节也可能导致细微差别。只要你的实现和 sklearn 的结果在合理的数值误差范围内比如 1e-8 量级就没有问题。真正确认实现正确的方法是用我前面提到的模拟数据你知道真实的 w 和 b 是多少模型学出来的 w 和 b 只要和真实值接近就说明代码逻辑没毛病。6. 进一步扩展给线性回归加上正则化6.1 岭回归在最小二乘基础上做手脚当特征很多或者特征之间存在相关性时X^T * X 接近奇异矩阵直接求逆会让权重变得非常大、对噪声极其敏感。岭回归Ridge Regression的改进是在对角线上加一个小的惩罚项theta (X^T * X alpha * I)^(-1) * X^T * y这里的 alpha 是正则化强度。加了这个项之后矩阵变得严格可逆而且权重被约束在比较小的范围模型的泛化能力通常更强。手动实现只在解析解那行代码里多一个项def fit_ridge(X, y, alpha0.1): X_b np.c_[np.ones((X.shape[0], 1)), X] n_features X_b.shape[1] theta np.linalg.inv(X_b.T X_b alpha * np.eye(n_features)) X_b.T y return theta从梯度下降的角度来看岭回归相当于在每一次参数更新时额外减去一个 alpha * theta相当于给权重泻火让它别长得太大。这也是为什么它经常出现在特征共线性严重的场景。6.2 Lasso 回归用稀疏性强制特征选择Lasso 回归用的是 L1 正则化惩罚的是权重的绝对值之和而不是平方和。它的最大特点是能把一部分权重直接压成 0从而实现特征选择。但 L1 正则化在原点处不可导梯度下降不能直接套用通常用坐标下降法Coordinate Descent去优化实现起来比岭回归要绕一些。如果想体验你可以用一个最简版的坐标下降思路每次固定其他权重只优化一个维度重复到收敛。不过对于入门来说先把岭回归吃透就足够理解正则化的价值了——L1 和 L2 最大的区别不在于实现细节而在于你希望模型稀释权重还是舍弃特征。手动实现到这里线性回归的骨架你已经完整搭建起来了生成数据、构造矩阵、设计损失、推导梯度、迭代训练、指标评估、调参排错、加正则化扩展。这套流程并不复杂但它涵盖了机器学习模型从无到有的完整生命周期。之后无论你是去学逻辑回归、决策树还是进到深度学习领域你会发现训练的思路始终是那一套定义损失、计算梯度、更新参数。把这条主线刻在脑子里比记住任何一个库的调用方式都重要。
返回列表