正规方程与梯度下降：小规模线性回归的闭式解实践

接触线性回归头几个月我一直默认所有模型参数都得靠梯度下降一步一挪地逼近。直到有一天同事看完我的训练代码后问了一句“你的设计矩阵才五百多行、十来个特征为什么不直接求逆矩阵拿解析解”这句话点醒了我——很多教程把梯度下降当成唯一解法却忽略了正规方程这种闭式解才是小规模回归问题的天然选择。这篇博文记录的是我基于波士顿房价数据集做房价预测的完整过程从数据体检、正规方程的数学推演到用数值计算库手写闭式解再到测试集指标和残差图的逐项验证。适合刚学完线性回归、想弄明白“自己写的解方程代码为什么能工作”的初学者也适合想把最小二乘原理抠细一点、搞懂算法选型边界的进阶读者。1. 上手前的账本为什么 506×14 的矩阵适合闭式解1.1 数据集规模与计算成本的粗算波士顿房价数据集包含 506 个样本、13 个原始特征加上偏置项后设计矩阵 X 的形状是 506×14。它的目标变量 MEDV 是自住房屋价格中位数单位是千美元。这个规模在机器学习里属于“小到不能再小”的级别正因如此它特别适合用来走一遍正规方程的完整流程。正规方程的计算量主要在两块构造 XᵀX 的矩阵乘法和对其求逆。前者复杂度是 O(m·n²)代入 m506、n14大约是 506×196也就是 10 万次左右的乘加运算后者用高斯消元求逆是 O(n³)大概 2700 多次操作。这个量级对现代 CPU 来说连“毫秒级”都算不上。如果换成梯度下降呢每迭代一次需要 O(m·n) 次运算大概 7000 次乘加。听起来单次更便宜但梯度下降通常要几千到几万步才能收敛而且还得反复调学习率、做特征缩放、检查收敛条件。在只有 14 个参数的模型里正规方程一次到位完全符合“杀鸡不用牛刀”的原则。1.2 正规方程相较于梯度下降的三大实际优势第一没有学习率这个超参数。梯度下降里 α 设大了震荡、设小了龟速而正规方程闭式解不存在这个问题对新手极其友好也少了一个“调参玄学”点。第二结果是确定性的。同样的数据、同样的预处理理论上能得到完全一致的系数梯度下降则受初始值、迭代次数和学习率影响两次运行结果可能有微小差异排查问题时多了一层干扰。第三无需保存中间状态。梯度下降要记录每轮迭代的损失用于画曲线、判断收敛正规方程在 14×14 的矩阵求逆上干脆利落代码量更少也更容易调试。提示这不代表正规方程永远比梯度下降好。它有一个硬伤——需要对 XᵀX 求逆复杂度随特征数 n 的三次方增长特征一旦上万内存和时间都吃不消。所以算法没有绝对优劣只有适用边界。1.3 特征尺度问题严格说不必须但影响数值稳定性很多文章说“正规方程不需要特征缩放”严格讲是对的因为闭式解在数学上不依赖特征单位。但在实际数值计算中如果某个特征量级是百万另一个是零点几XᵀX 的特征值会拉开极大差距造成病态矩阵求逆结果被舍入误差污染。波士顿数据的原始尺度还算温和但像 TAX 动辄几百、而 CHAS 只是 0/1 哑变量不统一量级时 XᵀX 的条件数会明显偏大。我下面的实现里仍然做了标准化——不是为了满足正规方程的“必须条件”而是为了让矩阵条件数更健康也让系数之间可比较。这里先埋个伏笔第五节我会量化对比标准化前后的条件数差异。2. 正规方程的推导过程从误差平方和到矩阵求逆2.1 目标函数与矩阵化的写法线性回归的假设函数可以写成 hθ(x) xᵀθ其中 x 是包含偏置项的 (n1) 维向量。把所有样本堆叠起来就是 Xθ于是损失函数写成J(θ) (1/2m) ||Xθ - y||²前面的 1/2m 是历史惯例1/2 是为了求导后消掉二次方带来的 2m 是为了把损失量纲变成“平均每个样本”。这一点很重要——因为 1/(2m) 对 θ 来说是常数正规方程的极值点不会受到它的影响。也就是说无论带不带这个系数最后解出的 θ 都一样。也许有人会问既然不影响结果为什么还要写因为如果用梯度下降做同样的任务1/m 这个平均因子直接影响步长与学习率的搭配而在闭式解里它只是推导路上一个“路过”的常数。把平方项展开就能看得更清楚||Xθ - y||² (Xθ - y)ᵀ(Xθ - y) θᵀXᵀXθ - 2θᵀXᵀy yᵀy对 θ 求梯度∇J(θ) (1/m)( XᵀXθ - Xᵀy )令梯度为零得到 XᵀXθ Xᵀy。这就是著名的正规方程。如果 XᵀX 可逆解就是θ (XᵀX)⁻¹Xᵀy2.2 几何视角最小二乘就是正交投影代数推导虽然简洁但我记牢它的辅助方式是几何图景。所有可能的预测值 Xθ 都落在 X 的列空间里。我们要找的是这个列空间中距离真实观测向量 y 最近的那个点。而“最近”的充要条件是残差向量 y - Xθ 与列空间正交。把正交性用矩阵语言写出来就是 Xᵀ(y - Xθ) 0也就是 Xᵀy XᵀXθ。这恰好就是正规方程。换句话说y 被分解成两部分位于列空间内的“可解释部分”Xθ与垂直于列空间的“残差部分”。最小二乘的直觉就是“把所有解释不了的东西甩到垂直方向上去”。用生活类比想象阳光垂直照向地面一个斜插着的杆子投出的影子就是“列空间投影”杆子顶端到影子的垂直连线距离就是残差。最小二乘要找的正是让这根垂直距离最短的投影位置。2.3 矩阵求导的两个小坑推导时最容易出错的环节有两个对二次型求导d(θᵀXᵀXθ)/dθ 2XᵀXθ。不少人会忘记系数 2导致最后方程差一个因子。对线性项求导d(θᵀXᵀy)/dθ Xᵀy。另外要注意正规方程结果唯一的前提是 XᵀX 满秩。当特征之间存在线性相关比如把类别变量做成一堆 0/1 且未去掉参考组XᵀX 就会奇异求逆直接报错这就是经典的“虚拟变量陷阱”。第六节我再展开讲。3. 波士顿房价数据的预处理身体体检决定结果上限3.1 字段清单与含义动手建模之前我习惯把每个字段的物理含义先摆出来。这个数据集里13 个特征加 1 个目标变量分别是字段含义备注CRIM城镇人均犯罪率分布严重右偏ZN大面积住宅用地占比大量样本为 0INDUS非零售商业用地占比与 NOX 强相关CHAS是否临河0/1 哑变量NOX一氧化氮浓度与 INDUS 强相关RM每套住宅平均房间数与房价正相关AGE建于 1940 年前的业主自住房比例与 DIS 负相关DIS到就业中心的加权距离与 AGE 负相关RAD高速公路可达性指数与 TAX 强相关TAX财产税率与 RAD 强相关PTRATIO学生与教师比例与房价负相关B早期统计的一项连续性指标原始定义存在伦理争议LSTAT低收入人口比例与房价强负相关MEDV房价中位数千美元目标变量关于 B 字段多说一句它的原始定义基于上世纪 70 年代的人口构成统计后来的研究者普遍认为该设计有问题这也是这个数据集近些年被多个平台下架、社区转而使用修正版的原因。我这次用的是社区整理的替代版本实际分析中也建议关注数据来源的说明。3.2 缺失、异常与分布形态数据加载进来后我会先用描述性统计过一遍。这个数据集通常没有缺失值但分布并不健康CRIM 的最大值和均值差一个数量级明显右偏ZN 有大量样本是 0说明不少城镇根本没有允许大占地住宅的区域CHAS 只有约 7% 的样本取 1属于典型的类别不平衡哑变量。对于右偏严重的列常见处理是取对数或做分位数变换。但对本篇这种“演示正规方程数学原理”的目标我选择保持原始特征只在连续特征上做标准化。原因是变换会改变特征的物理含义增加系数解释难度而且我们关心的重点不是冲刺最佳 RMSE而是理解闭式解的机理。注意如果是要参加比赛或追求上线精度右偏特征还是要处理的。经验上 CRIM 取 log 后模型对高犯罪率样本的拟合通常会好一些。3.3 相关矩阵里的共线性线索训练前看相关矩阵等于体检。我算完特征与目标变量的相关系数以及特征与特征之间的相关系数两组信息都很有价值与房价强相关的有 RM正相关约 0.7、LSTAT负相关约 -0.74、PTRATIO负相关约 -0.51、DIS正相关约 0.25。RM 越大、低收入人口比例越低房价越高符合直觉。特征之间的强相关有多组TAX 和 RAD 相关系数约 0.91NOX 和 INDUS 约 0.76AGE 和 DIS 约 -0.75。这组特征存在明显的多重共线性。共线性最直接的后果是XᵀX 接近奇异求逆时数值不稳系数的方差被大幅放大。你可能看到某个变量的系数符号和业务常识相悖或者同一个系数在标准化前后剧烈跳动。这都不是正规方程的“错”而是数据本身的结构问题我留到第六节专门处理。4. 手写实现与踩坑笔记从设计矩阵到 np.linalg.solve4.1 构造设计矩阵时最容易犯的错先给最终能跑的代码骨架去注释大概 60 行import numpy as np import pandas as pd df pd.read_csv(boston_housing.csv) X df.drop(columns[MEDV]).values y df[MEDV].values m, n X.shape # 1) 划分训练/测试注意要先划分再标准化 split int(m * 0.8) X_train, X_test X[:split], X[split:] y_train, y_test y[:split], y[split:] # 2) 标准化只 fit 训练集 mu, sigma X_train.mean(axis0), X_train.std(axis0) X_train (X_train - mu) / sigma X_test (X_test - mu) / sigma # 3) 加偏置列X_train_b 的形状是 (404, 14) X_train_b np.c_[np.ones(len(X_train)), X_train] X_test_b np.c_[np.ones(len(X_test)), X_test] # 4) 正规方程用 solve 而非显式求逆 theta np.linalg.solve(X_train_b.T X_train_b, X_train_b.T y_train) # 5) 预测与指标 y_pred X_test_b theta rmse np.sqrt(np.mean((y_test - y_pred) ** 2)) mae np.mean(np.abs(y_test - y_pred)) ss_res np.sum((y_test - y_pred) ** 2) ss_tot np.sum((y_test - y_test.mean()) ** 2) r2 1 - ss_res / ss_tot print(fRMSE{rmse:.3f} MAE{mae:.3f} R2{r2:.3f})这段代码我踩过的第一个坑就是忘记加偏置列。如果不加全 1 列正规方程算出的是一条过原点的超平面截距被强制为 0。这时 RMSE 会明显变差而且所有系数都被歪曲。检查方法很简单看训练集预测值的均值是否接近 y 的均值如果截距缺失残差往往会呈现系统性正负。提示有没有加偏置看 θ 的维度就知道——应该有 n1 个值第一个是截距项。4.2 用 np.linalg.solve 而不是显式求逆正规方程的教科书写法是 θ (XᵀX)⁻¹Xᵀy很多入门代码真的去写 np.linalg.inv(X.T X) X.T y。功能上没问题但工程习惯上我更推荐 np.linalg.solve(X.T X, X.T y)。原因有两个。其一solve 走的是 LU 分解或类似的高斯消元路线不需要真正算出逆矩阵运算量更小。对 14×14 的矩阵差别可以忽略但这是值得培养的习惯。其二显式求逆在病态矩阵上更容易放大舍入误差。inv 计算会先把逆的所有元素求出来再拿去乘中间多过一道手误差就多一次累积的机会solve 直接解线性方程组数值稳定性更好。4.3 标准化时机与数据泄漏标准化必须在划分训练/测试集之后、且只使用训练集的均值和标准差。如果不小心先对全量数据做了标准化再划分训练阶段就会看到测试集的信息属于一种常见但隐蔽的测试信息泄漏。直观理解是这样所谓测试集是在模拟“未来出现的新数据”。用全量数据的均值去标准化等于提前知道了未来数据的分布重心评估结果自然偏乐观。工程上的正规做法是把标准化器放进预处理管道里对训练数据 fit对训练和测试数据分别 transform。4.4 防止矩阵求逆直接炸掉的兜底方案如果真遇到 XᵀX 奇异比如某个特征是其他特征的线性组合np.linalg.solve 会直接抛出异常。这时候有两个选择改用伪逆np.linalg.pinv(XᵀX) 或直接对 X 求伪逆。伪逆在秩亏时仍能给出最小范数解虽然不一定唯一但至少不会中断流程。给对角线加一个小扰动用 (XᵀX λI)⁻¹ 代替 XᵀX⁻¹这就是岭回归闭式解的雏形。实际写代码时我会在求解前加一行 np.linalg.matrix_rank(X_train_b.T X_train_b)对比它和 n1 是否相等提前判断矩阵是否满秩而不是等报错再回头查。5. 测试集上的真实表现RMSE、R² 与残差图的合理解读5.1 我这次随机划分下的指标我用 80/20 划分不同随机种子下结果会在一个小区间内浮动。以我一次典型运行为例指标大致落在 RMSE 4.6 左右、MAE 3.3 左右、R² 0.78 左右单位是千美元。换算成原始价格平均误差约三千三百美元。考虑到这个数据集本身噪声不小这个水平符合经典文献的记录区间。我特意做了交叉验证把手写正规方程和机器学习库里底层走 SVD 分解的线性回归结果对比。两者系数完全一致RMSE 在小数点后三位内吻合。这个对照很有意义——它说明我没有写错公式也没有算错矩阵。5.2 残差图告诉我们什么把测试集残差画出来能发现两个明显模式。第一中低房价样本的残差集中在零附近而大于 30 千美元的高价样本残差几乎全部为负也就是预测值系统性偏低。原因是数据里高房价样本量少而且高房价区域往往由 LSTAT、RM 之外的隐藏因素决定线性模型扛不动这部分非线性。第二残差方差不恒定呈“漏斗形”拟合值越大残差波动越大。这就是异方差性。严格来说OLS 的无偏估计性质不依赖同方差但最优方差性质依赖。也就是说正规方程在这种数据上仍然给出“平方误差意义下”的最优解只是置信区间需要稳健修正不能直接拿普通标准误去解释。5.3 系数解释要小心共线性下的符号陷阱我这次得到的系数里RM 和 LSTAT 的符号符合直觉一个正、一个负。但 TAX 和 RAD 这两个高度相关变量的系数数值会随随机种子变化而有不小抖动偶尔还会出现符号翻转。这不是 bug——现实中它们俩挤占了彼此的方差正规方程虽然能把 XᵀX 求出来但逆矩阵对微小扰动极其敏感。所以阅读系数时不要轻率地说“税率每提升一个单位房价下降多少”。共线性场景下的单个系数并不代表因果效应最多只能说“在控制其他变量后该变量与房价的线性关联”。5.4 一个反直觉的检查点标准化会改变系数数值如果你分别用原始特征和标准化特征跑正规方程得到的 θ 数值会不同但换算回原始尺度的预测值应当一致。有人会把标准化前后的 θ 混为一谈以为模型变了。实际模型没变只是参数的“度量衡”变了。标准化后的系数衡量的是“该特征每变化一个标准差房价变化多少千美元”这让系数之间可以直接比大小也是我在这里做标准化的主要收益。打印 theta 时记得第一个数是截距项对应所有特征为标准差零的状态也就是特征均值的样本人群。6. 继续往里挖多重共线性、岭回归与算法选型边界6.1 条件数判断正规方程“危不危险”的体检指标正规方程的数值稳定性一个关键指标是 XᵀX 的条件数定义为最大奇异值与最小奇异值的比值σ_max / σ_min。条件数越大矩阵越“病态”求逆结果对输入的舍入误差越敏感。在我这次试验里打印 np.linalg.cond 得到的数值很清楚不做标准化时XᵀX 的条件数在千量级做完标准化后降到几百量级。虽然远没到危险区间经验上超过 10¹² 就很容易爆掉但已经能看出量级差异对数值的影响。经验法则如果条件数超过 10⁶~10⁸就要警惕求逆结果的可靠度。这时候先检查共线性再考虑是否要加正则化。提示条件数不是只看大小还要和你使用的浮点精度挂钩。双精度浮点大约能保证 15~16 位有效数字条件数每放大 10 倍求解结果就大约损失一位有效数字。6.2 岭回归闭式解同一条正规方程的“加了阻尼”版本当共线性把 XᵀX 推得很“扁”时一个自然的补救方案是对角线加一个小的常数项θ_ridge (XᵀX λI)⁻¹Xᵀy代码只改一行lam 1.0 A X_train_b.T X_train_b lam * np.eye(X_train_b.shape[1]) theta_ridge np.linalg.solve(A, X_train_b.T y_train)直觉理解XᵀX 是半正定矩阵特征值非负。当共线性导致某个特征值接近零时求逆等价于除以一个接近零的数结果自然被放大。加上 λI 等于给所有特征值抬升一个底把那些接近零的特征值拉离危险区代价是引入一点点偏差。我在同样的训练集上试了 λ1.0特征标准化后RMSE 几乎没有变化但 TAX 和 RAD 的系数抖动明显变小。这就是所谓的“偏差换方差”也是岭回归在共线性数据上常用的原因。6.3 什么时候该放弃正规方程复杂度与选型边界最后把选型边界说清楚特征数 n 很小几百以内正规方程几乎总是最优选择快、确定、无需调参。n 在几千到一万之间要留意 XᵀX 的存储是 O(n²)求逆是 O(n³)。内存和时间开始吃紧需要用 solve 而非 inv。n 上到十万、百万级绝对不要用正规方程。这时候就得转向随机梯度下降或小批量梯度下降它们的单步代价是 O(m·n)且可以流式处理数据。另外scikit-learn 的 LinearRegression 底层走的是 SVD 分解处理秩亏矩阵比直接求逆更稳健。生产环境中我通常直接用它但正因为亲手写过一遍闭式解才能判断它什么时候会失效、什么时候该切到岭回归。最后再说两点实操体会。第一我后来拿到新数据集会先跑一下 np.linalg.cond(XᵀX)只要条件数超过 10¹²就立刻意识到是共线性或量级问题而不是闷头调参。第二学正规方程最大的价值不在性能——毕竟生产中大部分人直接用封装好的库——而在于它把最小二乘的几何和代数两层含义打通了。之后你再用梯度下降、再用岭回归都会明白每一步其实都在解同一个问题在误差曲面最低点站稳脚跟。