如果你正在跟的教程或课程恰好走到“1.5”这一小节那你多半会在目录里看到这一行字最小二乘法求二元线性回归。一元回归的时候画个散点图、套一下斜率和截距公式还能硬扛一旦自变量变成两个公式里开始出现一堆下标和矩阵很多人就是从这里开始掉队的。这篇文章就专门把这一小节讲透先说明最小二乘法到底在优化什么再推出正规方程然后用一组五条数据完整手算三个系数最后用 Python 验证同一个结果。适合正在学统计、数据分析或机器学习入门的人也适合工作中要用回归模型、但一直没把矩阵化求解啃下来的同学。放心全程不需要复杂数学底子跟着算一遍就通了。1. 最小二乘法到底在求什么1.1 先别急着代公式弄清楚“二元线性回归”指什么“二元线性回归”这个说法在不同教材里其实有点歧义。一种理解是只有 x 和 y 两个变量拟合一条直线 y kx b另一种理解是包含两个自变量 x1、x2去预测一个因变量 y也就是最常见的学习时长 完成作业数预测期末成绩这种场景。我按第二种最普遍的情况展开y β0 β1·x1 β2·x2 ε其中 β0 是截距β1、β2 是偏回归系数ε 是误差项。很多新手第一次看到“偏回归系数”这五个字就懵了。其实用大白话讲β1 的意思是在 x2 保持不变的前提下x1 每增加 1 个单位y 平均会增加多少个单位。注意“在 x2 保持不变”这几个字这是偏回归系数和简单回归系数最核心的差别。你可以把它类比成体检报告年龄不变的情况下运动量每增加一小时血压平均下降多少。一旦把“控制变量”这个思路建立起来后面解读系数才不会出错。还要强调一点线性回归里的“线性”指的是参数 β0、β1、β2 是一次方的并不是说 x1、x2 本身必须是直线关系。即便你把 x1² 或 log(x2) 作为一个新变量塞进方程模型框架依然线性。这一节只是提醒真正的重点是接下来要解决的最核心问题三张系数表怎么定下来1.2 为什么偏偏是“最小二乘”你完全可以想象拟合一条直线或平面方法有很多。有人会说“让所有点尽量靠近”有人会说“让点和线之间垂直距离的总和最小”。但绝大多数教材最终选的是“残差平方和最小”这是有讲究的。残差就是真实值 y 与预测值 ŷ 的差。最小二乘法的目标是找一组 β让所有样本的残差平方和最小即 Q Σ(y_i - ŷ_i)² → min。平方有两个好处第一正负残差不会被相互抵消如果只求和不加平方一个 10 和一个 -10 就可能撑平为零完全失去衡量意义第二平方会放大那些偏差大的点让优化过程更在意离群样本。那为什么不干脆用所有残差的绝对值之和绝对值函数在 0 处不可导数学处理非常痛苦更重要的是最小二乘在“误差项服从正态分布”这个假设下和极大似然估计是等价的。换句话说如果你相信预测误差主要是 y 方向上的随机噪声那么平方误差目标就是最自然的那一个。这也是为什么到第 5 章我会专门提醒你留意异常值——平方项对异常点实在太敏感了。从几何上看最小二乘的目标函数是一个凸二次函数图像是一个碗状曲面只有一个全局最低点。所以它不像神经网络那样担心“陷入局部最优”后面用正规方程或者梯度下降去找最终结果理论上是一样的。搞清楚这个“为什么选平方”再看公式就不会觉得它是天上掉下来的规则。1.3 一个重要转折一元回归画直线二元回归画平面一元回归 y β0 β1x图像是在二维坐标平面上画一条直线你只需要调斜率和截距两个参数。二元回归 y β0 β1x1 β2x2图像就变成了三维空间里的一个平面三个参数共同决定这个平面怎么倾斜、怎么上下移动。可别小看这个变化。一元回归求你只需要一个公式斜率等于 x 与 y 的协方差除以 x 的方差。二元回归就麻烦多了因为 x1 和 x2 之间可能存在相关性你不能单独把 β1 抠出来单算必须三个方程联立求解。这也是为什么矩阵会出现在这一节里——矩阵本质上就是用一套统一的符号去解这种带多个变量的方程组。从思维习惯上你也得从“二维平面里的点”切换到“三维空间里的点”。学完之后你会发现这正是多元回归的核心场景用多个维度的信息去解释一个结果。以后再往三维、四维扩展数学框架完全一致只是不能画图罢了。2. 正规方程一次性求到最优解的数学底牌2.1 把 N 个样本写成一个矩阵方程先约定记号设计矩阵 X 是一个 n 行 p 列含截距项的矩阵每一行是一个样本第一列全部是 1代表截距项 β0后面分别是 x1、x2。β 是参数向量 [β0, β1, β2]ᵀy 是对应的因变量列向量。于是整个回归模型可以简写为 y Xβ ε。很多同学手算出错十有八九是忘记在设计矩阵里加那一列 1。这一列不是可选项它对应的是截距。用代码实现时如果你已经给 X 加了常数项列那么模型里就不能再让框架自动加截距反过来你直接用现成库的 fit_interceptTrue就不要再重复加常数列。这个细节后文会专门说。损失函数同样可以写成矩阵形式L(β) (y - Xβ)ᵀ(y - Xβ)。这里 (y - Xβ)ᵀ 是残差向量的转置两个向量相乘的结果正好是所有残差的平方和。接下来是整套线性回归里最漂亮的一段推导。先把损失函数展开L yᵀy - 2βᵀXᵀy βᵀXᵀXβ对 β 求导并令导数等于零利用 XᵀX 是对称矩阵这个性质可以得到-2Xᵀy 2XᵀXβ 0整理后就是著名的正规方程XᵀXβ Xᵀy如果 XᵀX 可逆最优解直接写成 β (XᵀX)⁻¹Xᵀy。这就是最小二乘法的解析解全程没有任何迭代。2.2 求解步骤与可解条件用正规方程求解操作上就三步第一算 XᵀX第二算 Xᵀy第三对 XᵀX 求逆再和 Xᵀy 乘起来。你如果把这三步背下来二元回归的系数问题就算解决一大半了。但有一个前提容易被忽略XᵀX 必须满秩。也就是说自变量之间不能出现完全线性相关比如 x2 恰好是 2 倍 x1那 XᵀX 的行列式就是 0没有逆矩阵。样本量也不能少于参数个数否则方程解不出来一个唯一值。实际数据里完全共线不常见但“接近共线”很常见。判断方法很简单先看 x1 和 x2 的相关系数超过 0.9 就要非常警惕同时可以算一下行列式行列式越小说明方程越病态系数对数据变化越敏感。真遇到不可逆或数值不稳定有几种常见替代方案用伪逆 np.linalg.pinv、改用最小二乘求解器 np.linalg.lstsq、删除冗余变量或者往 XᵀX 对角线上加一个小 λ 变成岭回归。这些方案不是今天的重点但值得记住名字后面用得上。2.3 正规方程 vs 梯度下降什么时候用什么很多讲机器学习的教程讲到这里会突然抛出一个词“梯度下降”然后开始调学习率。这里给你一个清晰的分工判断。正规方程是解析方法一步求解不需要学习率不需要迭代。它的主要成本是要做矩阵乘法和求逆。对于二元回归这种 p 只有 3 的场景正规方程是绝对的首选快、准、省心。梯度下降是迭代方法每次沿损失函数梯度的反方向走一小步适合样本量、特征量大到矩阵求逆非常吃力的场景也适合那些根本没有解析解的问题。表格对比会更直观一点方法思路优点缺点适合场景正规方程直接解线性方程组一步到位、无需调参求逆开销大特征多时慢样本不大、特征不多时首选梯度下降按梯度逐步逼近可扩展到大样本/大特征需要调学习率、迭代次数深度学习、大规模线性模型说白了二元线性回归这个规模根本不值得上梯度下降。你先把正规方程吃透将来遇到上百万样本再回头啃梯度下降也不迟。它们解的目标函数是同一个只是一个用代数一步到位一个用数值慢慢逼近。3. 手算五步彻底看懂系数是怎么算出来的3.1 一份小样本案例理论知识讲再多不亲手算一遍等于白看。我用一个非常典型的场景x1 是每周学习小时数x2 是每周练习完成份数y 是期末成绩。数据是我随手构造的示例一共 5 个样本足够展示完整的矩阵运算过程。样本x1学习小时x2练习份数y期末成绩1216223270343764538356489手工计算前先把后面矩阵里会反复出现的中间量整理出来。n5Σx120Σx213Σy380Σx1²90Σx2²39Σx1x259Σx1y1587Σx2y1035。这些数字看着多实际上就是求和一行行加出来即可。整理成表的好处是后面不管手算还是核对代码都能一眼看出问题在哪。3.2 完整推演从 X 矩阵到三个系数步骤 1构建设计矩阵。把原始数据拆成 X 和 yX 的每一行加一个常数 1X [[1, 2, 1], [1, 3, 2], [1, 4, 3], [1, 5, 3], [1, 6, 4]] y [62, 70, 76, 83, 89]ᵀ这一步在代码里就是 np.column_stack 两行的事手算时最容易漏掉的正是最左边这一列 1。步骤 2计算 XᵀX。这是一个 3×3 矩阵每个元素都可以从中间量直接填进去XᵀX [[5, 20, 13], [20, 90, 59], [13, 59, 39]]左上角是样本数 5第一行第二列是 Σx120第一行第三列是 Σx213中间 Σx1²90第二行第三列 Σx1x259右下角 Σx2²39。注意这个矩阵一定是对称的如果抄出来不对称那必然是前面的中间量算错了。步骤 3计算 Xᵀy得到一个 3×1 的列向量Xᵀy [380, 1587, 1035]ᵀ这三个数分别对应 Σy、Σx1y、Σx2y。步骤 4求 XᵀX 的逆矩阵。手算逆矩阵时先求行列式再算伴随矩阵最后除以行列式。我先算行列式det 5×(90×39 − 59×59) − 20×(20×39 − 59×13) 13×(20×59 − 90×13) 15计算行列式是最容易出错的一步代数余子式的正负号稍不注意就会翻车。如果你只是出于学习目的完整推一遍即可实际项目里没人手算这一步都是用软件完成。这里直接给出伴随矩阵的结果(XᵀX)⁻¹ (1/15) × [[29, -13, 10], [-13, 26, -35], [10, -35, 50]]我用一个简单方法验证它是不是正确的逆矩阵拿 XᵀX 的第一行 [5,20,13] 去乘以逆矩阵的第一列得到 (5×29 20×(−13) 13×10) / 15 145 − 260 130 15除以 15 正好等于 1说明对角线位置正确。用同样的方法检查非对角线位置结果都为 0。这个验证过程非常值得你在手算后做一遍。步骤 5把逆矩阵和 Xᵀy 相乘得到参数向量。按行展开β0 (29×380 − 13×1587 10×1035) / 15 739/15 ≈ 49.27β1 (−13×380 26×1587 − 35×1035) / 15 97/15 ≈ 6.47β2 (10×380 − 35×1587 50×1035) / 15 5/15 ≈ 0.33于是回归方程写出来就是y 49.27 6.47·x1 0.33·x2这个结果从业务角度很好解释在其他条件不变时每周学习时间每多 1 小时期末成绩平均多 6.47 分练习份数的边际贡献则只有 0.33 分左右因为在这个小样本里 x1 和 x2 高度相关x2 的独有信息很少。这个现象就是多重共线性的雏形到第 5 章再展开。3.3 残差、R² 与“这个方程到底好不好”得到系数之后一定要回代算一遍残差否则你根本不知道模型拟合得怎么样。把每个样本的 x1、x2 代入回归方程样本yŷ残差 y − ŷ16262.53−0.5327069.330.6737676.13−0.1348382.600.4058989.40−0.40残差有正有负求和接近 0这是最小二乘性质的直接体现。把所有残差平方后加起来得到 RSS 16/15 ≈ 1.07。再看 y 的波动样本均值为 76SST Σ(y − ȳ)² (62−76)² (70−76)² (76−76)² (83−76)² (89−76)² 450。决定系数 R² 1 − RSS/SST ≈ 0.9976。这个数字出现在示例数据里其实有点“过于漂亮”原因是样本量只有 5 个拟合自由度又多了一个模型对这种小样本天然容易做得很高。所以我必须提醒你训练集上的 R² 高只能说明拟合得很贴近这批数据不代表真实预测一定准。实践中的判断标准应该是交叉验证、测试集误差和残差图而不是单看一个 R²。4. 用 Python 五分钟跑通全流程附完整代码4.1 手写正规方程十几行实现理解了手算过程代码其实就是在复述刚才那五个步骤。用 numpy 可以非常简洁地实现import numpy as np X np.array([ [2, 1], [3, 2], [4, 3], [5, 3], [6, 4], ], dtypefloat) y np.array([62, 70, 76, 83, 89], dtypefloat) # 在设计矩阵中加一列常数项 1 X_design np.column_stack([np.ones(len(X)), X]) # 正规方程beta (XᵀX)⁻¹ Xᵀy beta np.linalg.inv(X_design.T X_design) X_design.T y print(beta) # 预期输出: [49.26666667, 6.46666667, 0.33333333]这里是矩阵乘法运算符X.T是转置np.column_stack负责拼接常数项列。运行结果和手算完全一致这本身就是对手算过程的最好验证。不过要提醒一句这种直接调用np.linalg.inv的写法教学意义很好工程上却不够稳妥。数值计算里解线性方程组远比直接求逆来得稳健。更推荐的写法是用np.linalg.lstsq它对接近病态的矩阵也能给出相对稳定结果beta_lstsq, residuals, rank, s np.linalg.lstsq(X_design, y, rcondNone) print(beta_lstsq) # 输出同样是 [49.26666667, 6.46666667, 0.33333333]如果你用lstsq发现结果和inv有细微差别不用慌那是浮点精度差异不是公式错了。4.2 用现成库交叉验证自己实现一次是学习用现成库是效率。opencv、sklearn 这类工具包已经封装好了线性回归模型代码几乎不需要思考from sklearn.linear_model import LinearRegression model LinearRegression().fit(X, y) print(model.intercept_, model.coef_) # 预期输出49.26666666666657, [6.46666666667, 0.33333333333]注意 sklearn 的LinearRegression默认fit_interceptTrue它会自动帮你处理截距项所以传入的X不要再加常数列。如果你自己加了常数列同时又让它去拟合截距等于给模型喂了两个常数项结果就会乱掉。这个坑我见过不少人踩。把三种方式放在一起对照方法β0β1β2手算49.26676.46670.3333numpy 正规方程49.26676.46670.3333sklearn49.26676.46670.3333三者结果一致说明思路没有跑偏。sklearn 底层实际上用的也不是裸求逆而是封装了更稳定的最小二乘求解器但最终解在数学上等价。4.3 可视化检查拟合效果参数算出来只是一个开始模型靠不靠谱还要看图。最常用的是残差图横轴是预测值 ŷ纵轴是残差 y − ŷ。如果模型假设成立残差点应该在 0 线附近随机散布没有明显的漏斗形或弯曲趋势。import matplotlib.pyplot as plt y_pred X_design beta residual y - y_pred plt.scatter(y_pred, residual) plt.axhline(0, colorgray, linestyle--) plt.xlabel(Predicted Value) plt.ylabel(Residual) plt.show()对于二元回归你还可以画三维散点和拟合平面直观看到数据点是不是大致落在一个斜平面上from mpl_toolkits.mplot3d import Axes3D import numpy as np x1s np.linspace(1.5, 6.5, 30) x2s np.linspace(0.5, 4.5, 30) X1g, X2g np.meshgrid(x1s, x2s) Yg beta[0] beta[1] * X1g beta[2] * X2g fig plt.figure() ax fig.add_subplot(111, projection3d) ax.scatter(X[:, 0], X[:, 1], y, colorred, labelsamples) ax.plot_surface(X1g, X2g, Yg, alpha0.4) ax.set_xlabel(x1) ax.set_ylabel(x2) ax.set_zlabel(y) plt.legend() plt.show()这类图最大的作用是帮你建立空间感拟合的本质就是找一个斜平面让所有点到平面的 y 方向距离平方和最小。可视化只是辅助手段别因为图形漂亮就跳过残差分析毕竟回归假设检验的主要依据在残差图上。5. 常见问题与避坑指南5.1 问题速查表学回归最容易遇到的几个症状基本都能在下面这张表里找到答案现象可能原因排查方法解决办法XᵀX 不可逆自变量完全线性相关、样本太少计算行列式、查看相关系数矩阵删冗余变量、用伪逆、加正则化系数符号与直觉相反多重共线性、异常值拖拽画散点图、相关性热图标准化、检查异常点、删除高度相关变量R² 接近 1 但预测很差过拟合样本量太少做交叉验证、划分测试集减少特征、收集更多数据正规方程算出巨大系数特征尺度差异大查看各特征的标准差标准化后重新建模手算与代码结果不一致设计矩阵忘加常数项列打印 X_design 第一列检查补常数列或关闭 fit_intercept这五类问题几乎覆盖了初学者在二元回归作业里遇到的 90% 的报错现场。重点提醒看到“奇异矩阵”“LinAlgError”不要慌先去看是不是两个变量几乎完全线性相关再看是不是样本量比特征数还少这两条是解不出唯一解的根本原因。5.2 三个容易被忽略的细节第一个细节是特征缩放。正规方程在数学上并不要求标准化但在实际计算中如果 x1 的量级是几千x2 的量级只有零点几XᵀX 里的数值会相差几个数量级求逆的精度会明显下降。做法很简单把每个特征减去均值、除以标准差让它们落到同一个量纲。标准化以后 β 的含义会从“原始单位每增加一单位”变成“标准差每增加一个单位”解读时要注意。第二个细节是异常值的影响。前面说过平方项会把残差放大所以最小二乘天然对离群点敏感。一个极端离群点可能把整个拟合平面拽向它导致所有正常样本的误差变大。数据进来先画箱线图或散点图看到明显离群点先判断是记录错误还是真实分布再决定剔不剔绝对不能无脑删。第三个细节是多重共线性。本文这组示例数据x1 和 x2 的相关系数已经接近 0.97属于高度相关。虽然还能求出系数但 XᵀX 的行列式只有 15说明它离病态很近。回到业务解释你会发现 x2 的系数很小且不稳定这是因为 x2 提供的增量信息太少。真正做项目时先算 X 的相关性矩阵发现高度相关的变量就只保留其中一个或者后续用岭回归这类带惩罚的模型。5.3 手算与代码不一致时的排查顺序如果你严格按照前面的步骤走了一遍发现手算结果和 Python 输出对不上我建议按下面这个顺序排查。第一步也是最高频的错误先看设计矩阵左边那列 1 有没有加进去。没有这一列截距项会被丢给回归方程自己硬凑结果一定不一样。第二步核对矩阵乘法顺序正规方程里必须是 XᵀX不是 XXᵀ第三步检查逆矩阵的代数余子式符号手算时最容易在第二行第三列这类位置丢负号。还有一个很隐蔽的坑是完成手算后又用 sklearn 验算却忘了它默认自动加截距。如果你的 X 已经手动加过常数列再用fit_interceptTrue等价于让模型有两个“第一列全为 1”的特征结果当然不同。最简单的解决办法是统一约定自己造常数列就用 numpy 跑用 sklearn就传原始 X。以我接触到的很多初学者反馈来看这个流程第一次完整走一遍确实耗时但矩阵乘法和求逆符号只要错一遍你就会永远记住 XᵀXβ Xᵀy 这个式子。后续不管使用哪种库看到它内部其实在做同一件事你对模型输出的信任感都不一样。这就是我把手算放在代码前面讲的原因理论这套底牌越早摸清后面学岭回归、Lasso、逻辑回归时就越不吃力。