★ FEATURED ARTICLE 我真见过太多人包括我自己早期调sklearn.linear_model.LinearRegression调了大半年模型是跑通了面试官一句“你给我讲讲最小二乘法怎么算的”直接愣住。这个标题看着简单但“底层代码实现”这六个字才是重点——它不是教你调包是让你把线性回归从数学推到代码、再从代码跑出结果的整条链路走通。这篇文章适合所有会一点numpy、但从来没自己手写过模型原理的人也适合已经会用sklearn、想补上底层盲区的同学。1. 为什么还要自己写线性回归的底层实现1.1 调库一年我发现自己回答不了这个问题先交代一下背景。有段时间我在做房价预测的模拟项目X数据集里有面积、房龄、周边配套评分几个特征当时我的做法简单粗暴from sklearn.linear_model import LinearRegressionfit一下predict一下完事。结果某次内部技术交流某导师问我“如果数据集只有一列特征你告诉我这条拟合直线的截距和斜率是怎么算出来的”我当时只能回答“最小二乘法”但具体的矩阵推导、甚至代码里怎么求逆矩阵一个字都说不出来。这个经历让我意识到一个很现实的问题调包能力不是模型能力。你当然可以一辈子用sklearn但一旦你要做在线学习数据一条条进来模型边跑边更新、要自定义损失函数、要把模型嵌入到不能装sklearn的嵌入式环境你不会底层实现就只能干瞪眼。所以这篇文章我不讲任何高级内容就盯着一件事把简单线性回归从“输入特征”到“输出预测”的每一步用numpy老老实实写出来。1.2 一条线性回归代码里的三个核心部件任何线性回归实现拆开看就三样东西假设函数Hypothesis、损失函数Loss、优化算法Optimizer。这个框架非常重要因为后面学逻辑回归、甚至学神经网络你都会发现它们是同一个骨架。先定义问题。假设我们有一个特征x和一个目标值y数据长这样(x₁, y₁), (x₂, y₂), ..., (xₙ, yₙ)。线性回归要找到一条直线y_pred w * x b这里的w是权重斜率b是偏置截距。写成向量形式就是y_pred θᵀ · X其中θ [b, w]X [1, x]。为什么要塞一个1进去因为这样b就可以看作一个特别的权重跟x统一处理后续矩阵运算方便得多。损失函数我们选均方误差MSEJ(θ) (1/m) · Σᵢ (y_predᵢ - yᵢ)²也就是预测值和真实值差的平方再取平均。为什么是平方不是绝对值因为平方函数处处可导而且惩罚大误差的力度更狠这正好符合大多数回归场景的直觉——差2个单位比差1个单位严重四倍。有了假设函数和损失函数剩下的问题就是怎么求出一组θ让J(θ)最小这就是优化算法的任务。经典做法有两条路正规方程Normal Equation和梯度下降Gradient Descent。1.3 手写实现到底难在哪说实话线性回归的底层代码不算难核心代码加起来不超过二十行。但难在三处第一矩阵维度。很多人推导都会一写代码就报ValueError: shapes (3,1) and (3,1) not aligned本质是对矩阵相乘的两个维度关系没建立直觉。第二数值稳定性。矩阵求逆不是万能的特征之间高度相关时逆矩阵会非常不稳定算出来的参数可能大得离谱。第三调参与验证。梯度下降的学习率、迭代次数、归一化策略每一项都依赖经验。你写完了怎么确认它对拿sklearn的结果对拍当然是最直接的办法但如果手写的结果和sklearn差一点你要知道差在哪。这篇文章的后半部分我会把这三块都过一遍。先来最简单的方案——正规方程。2. 直接求解最小二乘法的正规方程实现2.1 从平方误差到参数求解正规方程的推导很优雅。我们有损失函数J(θ) (1/m) · ||Xθ - y||²写成矩阵形式后对θ求导并令导数为0∂J/∂θ (2/m) · Xᵀ(Xθ - y) 0 XᵀXθ Xᵀy θ (XᵀX)⁻¹ · Xᵀy这就是正规方程。它告诉我们一个很漂亮的结论线性回归的损失函数是凸函数只有一个全局最小值直接一步到位求出来就行不需要像后面那些复杂模型一样一步步迭代。我举个例子。假设我们有五个数据点x [1, 2, 3, 4, 5]y [3, 5, 7, 9, 11]。肉眼可见y 2x 1那你用正规方程算出来的结果应该无限接近b ≈ 1, w ≈ 2。为什么是“无限接近”因为如果用真实采集的数据通常会有噪声干扰不会刚好一条直线。我这里用无噪声数据只是为了感受正规方程的精确性——它对凸函数是最小化问题的精确解。2.2 NumPy版正规方程完整代码直接上代码这是完整可运行的版本。import numpy as np def normalize(X): 标准化为了后续梯度下降做准备正规方程可不做 mean np.mean(X, axis0) std np.std(X, axis0) std[std 0] 1 # 防止除零 return (X - mean) / std, mean, std def normal_equation(X, y): 正规方程求解线性回归参数 X: shape (m, n)m是样本数n是特征数 y: shape (m,)目标值 返回 theta: shape (n,) # 添加偏置项在X最左边加一列全1 X_b np.c_[np.ones((X.shape[0], 1)), X] # 核心公式: theta (X^T X)^(-1) X^T y theta np.linalg.inv(X_b.T X_b) X_b.T y return theta # 构造数据y 3x 5 噪声 np.random.seed(42) X np.random.rand(100, 1) * 10 # 特征范围0~10 true_w 3.0 true_b 5.0 y true_w * X.flatten() true_b np.random.randn(100) * 1.5 theta normal_equation(X, y) print(截距 b:, theta[0]) print(斜率 w:, theta[1]) # 对比 sklearn from sklearn.linear_model import LinearRegression model LinearRegression() model.fit(X, y) print(sklearn 截距:, model.intercept_) print(sklearn 斜率:, model.coef_[0])我在本地跑出来的结果大致是b ≈ 4.95w ≈ 3.02和sklearn的结果在10⁻⁴级别完全一致。这对接完成说明实现正确。这里有个细节值得展开np.c_这个操作。它把一列全1拼到原特征矩阵的左边很多人不理解为什么非要多这一列。从几何上讲b和w本来就是两种角色——b不依赖任何输入特征对应的是直线在y轴上的截距w则和特征x相乘。数值计算上如果不加这列1θ里就没有位置放b了预测时Xθ少了常数项直线就变成一条必须过原点的线。以后你写任何线性模型逻辑回归、Softmax回归都要记得这一步。2.3 正规方程在什么情况下会“翻车”正规方程很直接但它有几个天生的短板我实际踩过之后印象特别深。第一个是矩阵求逆的计算复杂度。(XᵀX)是一个n×n矩阵n是特征数不是样本数。求逆的复杂度大致是O(n³)。如果你有两万维特征比如文本的TF-IDF向量这个矩阵求逆在普通笔记本上会算到怀疑人生——别问我怎么知道的问就是曾经跑过。第二个是矩阵不可逆问题。当特征之间存在多重共线性比如同时放入“摄氏温度”和“华氏温度”XᵀX的行列式趋近0逆矩阵数值上会变得极其不稳定算出来的θ动辄几百上千毫无意义。np.linalg.inv甚至可能直接报LinAlgError: Singular matrix。第三个是样本量小于特征量的情况。比如100个样本、500个特征XᵀX是500×500但秩最多只有100必然奇异。这种情况sklearn也救不了得靠正则化岭回归或者降到维度方法。所以正规方程适合的特征量级我的经验是最多几千维再多就建议换梯度下降或者np.linalg.pinv伪逆过渡一下。接下来就聊梯度下降——这也是“底层实现”四个字含金量最高的部分。3. 梯度下降另一个视角的底层实现3.1 梯度下降的本质是“沿着下坡路走”正规方程是直接算极值点梯度下降则是模拟一个人在山坡上往下走的过程。你要下到山谷最深处损失最小但你看不到全局地图只能感觉脚下的坡度然后每一步都往最陡的下坡方向迈一步反复走直到坡度接近0。这个“坡度”在数学上就是损失函数对参数的偏导数梯度。对线性回归的MSE损失梯度公式可以手推出来∂J/∂w (2/m) · Σᵢ (y_predᵢ - yᵢ) · xᵢ ∂J/∂b (2/m) · Σᵢ (y_predᵢ - yᵢ)矩阵化之后更简洁grad (2/m) · Xᵀ · (Xθ - y)这里X是加了偏置列的设计矩阵θ是[b, w]。每次迭代就做一次更新θ θ - learning_rate * grad学习率learning_rate就是“迈步的大小”。步子太大可能越过山谷直接冲上对面的山步子太小则要走到天荒地老。3.2 批量梯度下降的完整实现与收敛记录下面是一个完整实现包括数据生成、归一化、训练循环和损失记录方便你直接观察收敛过程。import numpy as np def prepare_data(X, y): 添加偏置列并保持为二维矩阵方便后续矩阵乘法 X_b np.c_[np.ones((X.shape[0], 1)), X] y y.reshape(-1, 1) return X_b, y def compute_loss(X_b, y, theta): m len(y) y_pred X_b theta return (1 / m) * np.sum((y_pred - y) ** 2) def batch_gradient_descent(X_b, y, theta_init, learning_rate0.01, n_iterations1000): m len(y) theta theta_init.copy() loss_history [] for iteration in range(n_iterations): y_pred X_b theta grad (2 / m) * (X_b.T (y_pred - y)) theta theta - learning_rate * grad loss compute_loss(X_b, y, theta) loss_history.append(loss) if iteration % 200 0: print(fIter {iteration:4d} | loss: {loss:.6f}) return theta, loss_history # 生成带噪声的线性数据 np.random.seed(42) X np.random.rand(200, 1) * 5 y 2.5 * X.flatten() 1.2 np.random.randn(200) * 0.8 # 数据标准化梯度下降的必需品 X_scaled, mean, std (X - X.mean()) / X.std(), X.mean(), X.std() X_b, y prepare_data(X_scaled, y) theta_init np.zeros((X_b.shape[1], 1)) theta, loss_history batch_gradient_descent(X_b, y, theta_init, learning_rate0.1, n_iterations1000) # 把标准化后的参数还原到原始尺度 w_original theta[1][0] / std b_original theta[0][0] - w_original * mean print(原始尺度下的 b:, b_original) print(原始尺度下的 w:, w_original) # 对拍验证 from sklearn.linear_model import LinearRegression model LinearRegression() model.fit(X, y) print(sklearn b:, model.intercept_) print(sklearn w:, model.coef_[0])训练过程里一个有意思的点是loss_history前几十轮下降非常快之后慢慢趋于平稳最后无限逼近正规方程的解。我第一次跑的时候打印250轮损失从85.7一路降到1.8直观感受到了什么叫“凸函数的梯度下降一定收敛”。这里要强调一个关键操作——数据标准化。我先把X减均值除以标准差训练完后再把参数还原。这么做有三个原因第一如果X的量纲很大比如房价几百万而y很小梯度的数值会非常大训练直接发散第二标准化后所有特征尺度一致学习率这单个超参数才能对全部维度都有合适的行为第三坐标轴等比缩放后损失函数的等高线从狭长椭圆变成近似正圆梯度方向更直接指向圆心收敛更快。3.3 学习率、归一化和收敛判据学习率是梯度下降里最需要“手感和经验”的超参数。我在模拟项目里测试过一组对照记录下来非常直观学习率迭代500次后的行为0.001损失缓慢下降500次未完全收敛0.01收敛正常曲线平滑0.1收敛较快曲线稳定0.5前几步骤剧烈震荡后期勉强收敛1.0损失爆炸性增长直接发散这里还有个小技巧如何写一个“自适应收敛”的训练循环不要写死n_iterations1000而是设置一个阈值当两次迭代的损失变化小于1e-6就提前停。工程实现时这个比硬编码迭代数靠谱得多下面给个片段def gradient_descent_with_early_stop(X_b, y, learning_rate0.1, tol1e-6, max_iter5000): m len(y) theta np.zeros((X_b.shape[1], 1)) loss_prev compute_loss(X_b, y, theta) for i in range(max_iter): grad (2 / m) * (X_b.T (X_b theta - y)) theta theta - learning_rate * grad loss compute_loss(X_b, y, theta) if abs(loss_prev - loss) tol: print(f收敛于第 {i1} 次迭代) break loss_prev loss return theta实际经验是如果损失曲线长时间“平着走”先别急着加迭代次数优先检查梯度是否计算错误、数据是否归一化。这两个检查能解决我遇到过的九成训练不收敛问题。4. 从一维到多维向量化实现与特征工程4.1 用一个矩阵运算统一所有维度前面所有代码都只用了一个特征x。但真实场景没有这么善良——我处理过的数据分析任务里最少的也有三四个特征多的几十上百个。好消息是上面代码几乎一行不用改就能从一元线性回归扩展到多元线性回归。原因在于numpy的矩阵运算天然支持多维。你只需要把X从(m, 1)变成(m, n)剩下的X_b np.c_[np.ones((m, 1)), X]、X_b theta、X_b.T (X_b theta - y)全部照常工作。维度变化如下X_b(m, n1)theta(n1, 1)X_b theta(m, 1)就是所有样本的预测值X_b.T (X_b theta - y)(n1, 1)对应每个特征的梯度分量所以向量化不只是“写得整齐一点”它直接决定了你代码能不能跑起来。如果你用三层嵌套for循环去更新参数当数据量到几万条的时候就非常痛苦换成矩阵乘法底层是高度优化的BLAS线性代数库快几十倍不止。这也是为什么我建议所有入门者从一开始就养成矩阵思维把样本当成一个整体别一个样本一个样本地循环。4.2 真实场景里非数值特征怎么处理多元线性回归还有一个在底层实现教程里几乎不讲的隐形门槛特征预处理。我自己在第一个真实数据集上就栽过跟头。当时特征里有个“城区类别”取值是[A区, B区, C区]我直接把字符串映射成[0, 1, 2]就扔进模型结果出来的拟合结果很差因为这等于给类别强加了一个“B区是A区的两倍”的荒谬假设。正确的做法是做独热编码One-Hot Encoding每个类别变成一个独立的0/1列。用numpy实现起来也很直接def one_hot_encode(labels): unique np.unique(labels) matrix np.zeros((len(labels), len(unique))) for i, label in enumerate(labels): matrix[i, np.where(unique label)[0][0]] 1 return matrix如果你可以用pandaspd.get_dummies()一行搞定。但自己写一遍能帮你理解独热编码的本质它是一个合法的特征变换把名义变量映射到多维欧氏空间每列之间相互正交避免人为引入顺序关系。除此之外还有两个常见操作。一个是特征缩放除了标准化z-score还有MinMaxScaler把特征压到[0, 1]区间适合特征分布比较均匀的情况。另一个是特征组合比如房间面积和楼层单独作为特征时效果一般但面积 × 楼层往往有更强的解释性。这些在sklearn里有专门的PolynomialFeatures底层原理就是numpy的列拼接和幂运算理解了本质后你完全可以在数据预处理环节自己构造。4.3 向量化实现的性能收益我在一个模拟数据集上做过简单对比m10000n50同样跑1000次梯度下降未优化的循环实现耗时约38秒向量化矩阵实现不到0.3秒。差距超过百倍。这个差距在算法竞赛、实时预测场景里就是能不能用的区别。向量化实现的额外好处是代码更接近数学公式可读性和可调试性都更好。你看上面grad (2/m) * (X_b.T (X_b theta - y))直接对应梯度公式想出错都难。而for循环版本里面一堆下标出错概率成倍增长。5. 验证与排错如何确认自己写对了5.1 用梯度检查定位实现错误手写模型最大的恐惧是“不知道错在哪”因为训练能跑通并不代表参数对不对。一个可能的情况是你梯度算错了但恰好损失还在下降拟合效果“看起来还行”实际性能已经离最优差很远。这时候最有效的工具就是数值梯度检查Gradient Checking。原理很简单导数的定义是极限那我们用很小的eps去近似它。对参数θᵢ分别计算θᵢeps和θᵢ-eps的损失然后numeric_grad_i ≈ (J(θ eps·eᵢ) - J(θ - eps·eᵢ)) / (2·eps)如果这个数值梯度和我们解析算出来的梯度在每个维度上都接近相对误差小于1e-5量级那梯度计算基本没错。放一段实现def gradient_check(X_b, y, theta, eps1e-6): m len(y) analytic_grad (2 / m) * (X_b.T (X_b theta - y)) numeric_grad np.zeros_like(theta) for i in range(len(theta)): theta_plus theta.copy() theta_minus theta.copy() theta_plus[i][0] eps theta_minus[i][0] - eps loss_plus compute_loss(X_b, y, theta_plus) loss_minus compute_loss(X_b, y, theta_minus) numeric_grad[i][0] (loss_plus - loss_minus) / (2 * eps) diff np.linalg.norm(analytic_grad - numeric_grad) / (np.linalg.norm(analytic_grad) np.linalg.norm(numeric_grad)) print(相对误差:, diff) return diff 1e-6梯度检查这个习惯我从线性回归一直保留到后面的神经网络实现每次新写一个模型、不能确认反向传播是否正确的时候跑一次梯度检查就能定位问题比盲目调参高效得多。5.2 回归指标怎么选模型写完了总得有一个判断好坏的量化标准。线性回归最经典的是决定系数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)R² 0.85意味着模型解释了85%的方差波动。但要注意R²对特征数量敏感特征越多R²天然越高哪怕那些特征毫无意义。所以跨模型比较时我会额外看调整后的R²。另外实际业务中我也会习惯性看RMSE均方根误差因为这个指标的单位和y一致可以直接说“预测房价平均偏差多少钱”老板听得懂比R²直观得多。5.3 我踩过的三个典型坑第一个坑是忘记归一化直接跑梯度下降。我用未归一化的数据算出来的损失曲线像心电图数字飙到1e12当时第一反应是学习率太大调到1e-9还是抖人麻了。后来才反应过来是特征尺度问题——一个特征取值几万一个特征取值零点几梯度方向被大尺度特征完全绑架。第二个坑是把标准化参数直接当最终结果输出。我第一次写模型时在归一化的数据上训练得到w0.87就直接拿去给别人解释“面积每平米房价涨0.87万”完全忘了这是在标准化尺度下的数字。正确做法是先还原到原始尺度再给结论。第三个坑是用训练集评分代替泛化性能。我犯过最愚蠢的错误在训练集上R²达到0.98开心得不行结果换个新数据直接掉到0.6。后来才意识到线性回归的基础假设里有一条“训练集和测试集独立同分布”数据分布一旦偏移模型立刻现原形。所以每次训练我都会把数据切成训练/验证两部分验证部分的分数才是真正能相信的。写在最后这一套手写线性回归的流程走下来最大的收获不是会写那十几行代码而是建立了一个“模型是怎么被训练出来的”完整心智模型。以后再遇到逻辑回归你只需要把假设函数换成sigmoid、损失函数换成交叉熵剩下的梯度下降、归一化、验证流程全部可以照搬。这篇文章里提到的正规方程、批量梯度下降、梯度检查你面试的时候能亲手推导一遍和只会说概念完全是两种状态。最后再分享一个小技巧如果你拿到的数据集特征量很大正规方程慢得跑不动可以先试试在维度方向上做PCA降维再跑正规方程速度能快很多。遇到多重共线性也别慌给XᵀX的对角线加上一个小λ这就是岭回归你的底层实现立刻升级了一版。希望这篇分享能帮你把线性回归这块地基真正打扎实。
阅读完成 · 觉得有帮助？
