★ FEATURED ARTICLE 波士顿房价预测这个项目入门机器学习的朋友基本都绕不过去。而我更想说的是越是那种“代码实现看起来只有几行”的项目越值得把原理抠明白。拿正规方程Normal Equation来解线性回归很多时候代码就是矩阵乘法和求逆跑完一片岁月静好但真到了数据多一点、特征乱一点的时候各种奇怪问题才冒出来。这篇文章我准备把这个经典项目完完整整拆开从数学推导到数据预处理从最小实现到踩坑排查全部按实际做项目的思路来讲目标是看完你不仅能跑通波士顿房价预测还能真正理解每一行代码背后的逻辑。需要说明的是波士顿房价数据集在新版的 scikit-learn 中已经被移除了所以这里我直接用 UCI 上的原始数据文件来做顺便也模拟一下真实工程里“拿到原始数据自己清洗”的过程。1. 为什么波士顿房价预测还在用正规方程方案选型拆解1.1 正规方程是什么一个能直接算出来的闭式解正规方程是线性回归模型参数的一种解析式解法。所谓“解析解”意思是不需要像梯度下降那样一点一点去逼近而是通过数学公式直接把参数算出来。线性回归要拟合的函数长这样h_θ(x) θ_0 θ_1x_1 θ_2x_2 ... θ_nx_n写成矩阵形式就是 h_θ(x) Xθ其中 X 是样本矩阵每一行是一个样本每一列是一个特征θ 是我们要找的模型参数向量。正规方程的核心结果就一个公式θ (X^T X)^(-1) X^T y看着很短对吧实际代码也就一行。但这一行是整个线性回归里最有信息量的一行因为它直接把“最小化误差”这个优化问题变成了一个纯粹的线性代数问题。我用这个项目去讲正规方程主要是因为它非常适合展示这个公式背后的逻辑。506个样本、13个特征数据量不大不小既不会因为矩阵太大导致半天算不出来也不至于样本太少看不出统计规律。而且数据本身存在多重共线性这对讲解“什么时候正规方程会失效”提供了绝佳素材。1.2 和梯度下降怎么选一张表看清边界很多新手第一次接触线性回归时会困惑为什么有的资料讲梯度下降有的资料讲正规方程到底该学哪个我的观点是两者都要懂但理解优先级不同。梯度下降是给“大场面”准备的正规方程则能让你快速理解线性回归的本质。维度正规方程梯度下降参数求解方式直接计算解析解迭代逼近需要调整超参数不需要没有学习率需要调学习率需要特征缩放不是必需的强烈建议计算复杂度O(n^3)与特征数量强相关O(k·m·n)与迭代次数相关适用场景特征数量相对较少特征数量大、样本量大这里有一个关键判断标准特征数量 n 是否超过了 1 万左右。如果 n 在几千以内正规方程通常很快因为矩阵求逆的耗时还算可控。但当特征数上万甚至更高时O(n^3) 的复杂度会让人非常难受训练一次可能需要几十秒甚至几分钟这时候用梯度下降会舒服很多。波士顿房价这个项目特征一共就13个即便做了多项式扩展也就几十个用正规方程必然是首选。这也是我拿它来演示正规方程的原因之一场景合适数学背景清晰不会因为算力问题喧宾夺主。2. 正规方程的数学原理从损失函数到矩阵推导2.1 最小二乘目标在做什么很多人一看到公式就开始背却不理解为什么要求出那个 θ。我们先用一个生活化的例子建立直觉。假设你在估算一套房子的价格手上有几个特征面积、房间数、地段评分。你给出的预测是这些特征的线性组合。预测值和真实值之间必然有误差。最小二乘法的思路非常朴素——找到一组系数让所有样本的预测值和真实值之差的平方和最小。这个“差的平方和”就是损失函数J(θ) (1/2) ∑ (h(x^i) - y^i)²前面那个 1/2 纯粹是为了后面求导方便凑出来的不影响最小值所在的位置。写成矩阵形式J(θ) (1/2)(Xθ - y)^T (Xθ - y)为什么要用平方而不是直接用绝对值两个原因一是平方函数处处可导方便用微积分找极值二是平方会放大较大误差的权重让模型更“重视”那些偏差大的样本。当然代价是对异常值敏感这个在后面的章节我会单独提到。2.2 一次求导得出模型参数既然 J(θ) 是关于 θ 的函数那么它的最小值出现在梯度为零的地方。对 J(θ) 求关于 θ 的偏导并令其等于 0∇J(θ) X^T(Xθ - y) 0移项X^T Xθ X^T y如果 X^T X 是可逆矩阵那么θ (X^T X)^(-1) X^T y这就是正规方程的完整推导过程。整个过程没有任何迭代没有学习率就是一个冷冰冰的矩阵运算。所有“学习”的因素都隐藏在 X^T X 的求逆里面。这里想补一句初学者容易忽略的点X 矩阵里通常要加一列全 1 的向量对应的是截距项 θ_0。如果不加你的模型就被强制穿过原点预测能力会大打折扣。这一点我后面在代码里也会再强调一次。2.3 求逆的替代方案解线性方程更稳正规方程的公式里有个显式的矩阵求逆(X^T X)^(-1)。但真要写代码时我一般不直接调 np.linalg.inv而是用 np.linalg.solve 来解线性方程组。原因有两个一是数值稳定性。直接求逆矩阵的算法如伴随矩阵法或初等行变换在矩阵条件数很大时误差会被放大而 solve 方法底层用的是 LU 分解数值表现要稳健得多。二是性能。求解线性方程组 Aθ b 的计算量通常比显式求逆再乘向量要小而且更不容易因为舍入误差导致结果漂移。实际工程中我用得最多的其实是 np.linalg.lstsq它走的是 SVD奇异值分解即使 X^T X 不可逆也能给出最小范数解属于“兜底稳妥”方案。但如果要配合教学让大家理解正规方程的原始形态用 solve 或 pinv 都合适。3. 波士顿房价数据集实战准备字段、探索与预处理3.1 数据集字段与含义波士顿房价数据集一共 506 条样本每条样本包含 13 个特征和 1 个目标变量。这个数据集来自美国波士顿地区目标变量 MEDV 是区域内自住房的中位数价格单位是千美元。字段名含义CRIM城镇人均犯罪率ZN占地面积超过2.5万平方英尺的住宅用地比例INDUS城镇非零售商业用地比例CHAS是否临查尔斯河1为临河0为不临河NOX一氧化氮浓度RM平均每个住宅的房间数AGE1940年之前建成的自住单元比例DIS到波士顿就业中心的加权距离RAD到高速公路的可达性指数TAX每1万美元的不动产税率PTRATIO城镇学生与教师比例B一个按城镇人口计算的比例相关指标LSTAT低收入人群占比百分数MEDV自住房中位数价格单位千美元目标变量从工程角度看这个数据集的特色在于特征量级差异极大。比如 TAX 的数值普遍在 200 到 700 之间而有些比例类特征可能只在小数点后两位这种量级跨度对正规方程的数值稳定性是一个考验。3.2 数据探索需要关注什么写代码之前我习惯先做三件事看维度、看缺失、看分布。看维度确认数据形状是 (506, 14)避免读文件时列错位。看缺失虽然这个数据集没有缺失值但真实数据里这一步绝不能省。看分布对 MEDV 做一次直方图观察能发现它基本符合正态分布的形态但右尾略长。这提醒我在评估模型时不能单看平均误差还要结合误差分布来分析。还有一个容易被忽视的点特征之间的共线性。比如 TAX 和 RAD 之间的相关性往往很高因为它们都跟区域可达性有关。这种共线性不会让正规方程完全算不出来但会让参数估计的方差变大导致结果不稳定。观察相关矩阵热力图能帮你提前意识到这个问题。3.3 预处理该做哪些针对正规方程预处理的原则和梯度下降不一样需要单独梳理一下。特征缩放不是必须的。因为正规方程没有迭代过程不需要通过缩放来加速收敛。但量级差异过大会影响 X^T X 的条件数造成求解结果波动。所以我的做法是对数值型特征做标准化减均值除标准差特别是后边还要做多项式扩展的时候标准化能避免高次项把数值范围撑爆。标准化有一个细节——必须先拆分训练集和测试集再在训练集上计算均值和标准差最后用训练集的统计量去变换测试集。如果先对整个数据集做标准化再划分测试集的信息就提前“泄漏”进了训练过程评估结果会虚高。虚拟变量的处理通常只需要关注 CHAS它本身就是 0/1 二值变量不需要额外 one-hot。其他特征都是连续型直接进入模型即可。4. 代码实现从零手写正规方程完成房价预测4.1 数据加载与格式整理因为新版 scikit-learn 移除了load_boston我直接从 UCI 仓库读取原始文件。这个文件是纯文本格式每行 14 个数值用空格分开。import numpy as np import pandas as pd # 从UCI读取波士顿房价原始数据 df pd.read_csv( https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data, headerNone, delim_whitespaceTrue ) # 为每一列补充字段名 df.columns [ CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, B, LSTAT, MEDV ] print(df.shape) print(df.head())如果网络不稳定读不出来可以先把 housing.data 下载到本地然后改成read_csv(housing.data, ...)。这里提醒一句原始文件没有表头headerNone一定不能漏否则第一行数据会被当成列名。4.2 核心代码矩阵运算一行求出模型参数拿到数据后先把特征矩阵和目标变量拆开然后给 X 添加一列全 1 的截距项最后调用正规方程公式求解。X df.drop(MEDV, axis1).values y df[MEDV].values # 添加截距项 X_b np.c_[np.ones((X.shape[0], 1)), X] # 正规方程求解 theta np.linalg.solve(X_b.T X_b, X_b.T y) print(theta)如果你更想贴近教科书里的原始公式可以这样写theta np.linalg.inv(X_b.T X_b) X_b.T y两种写法结果一致。但实际项目中我更推荐solve一方面数值更稳另一方面当矩阵接近奇异时solve会直接报错反而帮你提早发现数据问题。预测和评估也很直接y_pred X_b theta # 计算RMSE与R2 rmse np.sqrt(np.mean((y_pred - y) ** 2)) r2 1 - np.sum((y - y_pred) ** 2) / np.sum((y - np.mean(y)) ** 2) print(RMSE:, rmse) print(R2:, r2)在全量数据上跑完RMSE 大约在 4.7 左右R2 大约 0.73 到 0.75。也就是说平均预测误差约为 4700 美元模型解释了目标变量约七成多的方差。对于一个没有任何特征工程的普通线性模型这个结果属于正常水平。4.3 模型评估与结果解读全量拟合看起来不错但真实的建模流程会要求拆出测试集。我用一个简单的 80/20 划分来演示。np.random.seed(42) shuffle_idx np.random.permutation(len(df)) train_idx shuffle_idx[:int(len(df) * 0.8)] test_idx shuffle_idx[int(len(df) * 0.8):] X_train, X_test X[train_idx], X[test_idx] y_train, y_test y[train_idx], y[test_idx] X_b_train np.c_[np.ones((len(X_train), 1)), X_train] X_b_test np.c_[np.ones((len(X_test), 1)), X_test] theta_train np.linalg.solve(X_b_train.T X_b_train, X_b_train.T y_train) y_test_pred X_b_test theta_train rmse_test np.sqrt(np.mean((y_test_pred - y_test) ** 2)) r2_test 1 - np.sum((y_test - y_test_pred) ** 2) / np.sum((y_test - np.mean(y_test)) ** 2) print(Test RMSE:, rmse_test) print(Test R2:, r2_test)因为样本量只有 506划分方式不同测试集指标波动会比较明显。我自己实测下来R2 在 0.65 到 0.85 之间浮动都是正常的这恰好说明小样本下模型评估的不稳定性。所以跑这个项目时不要过分迷信单次划分的结果多换几个随机种子看看分布才是正确的打开方式。还有一个很容易踩的细节划分数据时一定要保证特征和标签按同一索引对齐。如果直接 shuffle 了特征矩阵却忘了 shuffle 目标变量模型学到的就是错位的数据结果会极其离谱。我在第一次写这段代码时就栽过这个跟头。4.4 高级扩展多项式特征进一步提升效果波士顿房价数据集中有不少特征与房价的关系并不是线性的比如 LSTAT 和房价之间近似反比例关系。这时候可以考虑做多项式扩展让模型能够拟合曲线。from sklearn.preprocessing import PolynomialFeatures poly PolynomialFeatures(degree2, include_biasFalse) X_poly poly.fit_transform(X) X_b_poly np.c_[np.ones((X_poly.shape[0], 1)), X_poly] theta_poly np.linalg.solve(X_b_poly.T X_b_poly, X_b_poly.T y) y_pred_poly X_b_poly theta_poly rmse_poly np.sqrt(np.mean((y_pred_poly - y) ** 2)) r2_poly 1 - np.sum((y - y_pred_poly) ** 2) / np.sum((y - np.mean(y)) ** 2) print(Poly RMSE:, rmse_poly) print(Poly R2:, r2_poly)把 13 个特征扩展成 2 次多项式后特征数量会膨胀到一百多维。正规方程的求逆复杂度是 O(n^3)但一百多维的矩阵求逆依然非常快这就是小特征集下用正规方程做多项式回归很舒服的原因。实测 RMSE 能降到 3.5 左右R2 可以超过 0.85。但这里要提醒一下多项式特征 正规方程的组合很容易过拟合。特征全部参与拟合模型可以做到训练集误差很低但测试集表现可能反而变差。解决办法是在损失函数中加入正则化项也就是下一章要讲的岭回归闭式解。5. 踩坑记录正规方程实现中的常见问题速查5.1 矩阵奇异时怎么办正规方程计算的硬性前提是 X^T X 可逆。可逆性被破坏通常有两个原因一是特征之间存在完全线性相关比如你把列“面积_平米”和列“面积_亩”同时放进模型两列只差一个常数倍数X^T X 的行列式为零。二是样本数小于特征数。比如做了 3 次多项式扩展之后特征数超过样本数这时 X^T X 必然是奇异的。遇到这种情况怎么办先检查特征是否有多余的线性组合尽量删除冗余列。如果删完还是奇异就用伪逆theta np.linalg.pinv(X_b.T X_b) X_b.T y更稳妥的方案是直接上 SVD 求解theta, residuals, rank, singular np.linalg.lstsq(X_b, y, rcondNone)这两种方式都能在矩阵奇异时给出一个可用解。但请注意你的模型此时可能已经过拟合了伪逆给出的参数范数通常也很大要结合正则化一起用。5.2 特征量级差异对结果影响大吗正规方程不像梯度下降那样“怕”特征量级因为不涉及梯度更新的步长问题。但量级差异会通过条件数影响数值稳定性。举个例子TAX 的数值范围是几百而有些比例类特征在小数点左右浮动两者差距一千倍以上。计算 X^T X 时TAX 那一列产生的数值会占据主导导致矩阵接近奇异。此时求逆的结果可能在数值上不稳定稍微换一批数据参数就大幅波动。我的习惯是对所有连续特征做标准化。这不只是为了稳定数值还有一个附带好处标准化之后模型参数的大小可以粗略反映特征重要性这对解释模型很有帮助。标准化的代码非常简单from sklearn.preprocessing import StandardScaler scaler StandardScaler() X_scaled scaler.fit_transform(X)要提醒的是多项式扩展如果接在标准化之后做高次项会自然保持合理的量级。如果顺序反了先做多项式再标准化虽然也能用但中间会多出很多不必要的极大值计算。5.3 多特征下如何防止过拟合岭回归的闭式解正规方程天然会把训练集误差压到很低尤其特征数多、样本数少的时候过拟合几乎是必然事件。解决办法是在损失函数里加一个 L2 正则项J(θ) (1/2)||Xθ - y||² (λ/2)||θ||²求导后得到的闭式解变成θ (X^T X λI)^(-1) X^T y注意I 是单位矩阵但通常不对截距项施加惩罚所以实际实现中会把 I 的第一个对角线元素改成 0。l2_lambda 1.0 I np.eye(X_b.shape[1]) I[0, 0] 0 # 不惩罚截距项 theta_ridge np.linalg.solve(X_b.T X_b l2_lambda * I, X_b.T y)这就是岭回归的解析解。它和普通正规方程只差了一个 λI但作用非常大加了 λ 之后即使 X^T X 接近奇异矩阵 X^T X λI 也一定是可逆的同时参数会被压缩过拟合现象明显缓解。在波士顿房价项目里做多项式扩展之后我建议试一下 λ 从 0.01 到 10 的几组取值你会发现测试集误差呈先下降后上升的趋势。这个“最小点”对应的 λ 就是你要找的当然更严谨的做法是交叉验证来选。尾声一点经验之谈把这个项目从头到尾做一遍你会发现自己收获的远不止“会调一个线性回归模型”这么简单。正规方程的可贵之处在于它用最直白的方式揭示了线性回归的本质拟合一条线让误差的平方和最小。你亲手推过公式、亲手写过求逆代码之后再去看 sklearn 里封装好的 LinearRegression理解完全是两个层级。最后分享一个我自己的小习惯每次实现完正规方程我都会顺手打印一下 X^T X 的条件数或者直接看 lstsq 返回的奇异值。这个数能非常直接地告诉你当前数据是不是病态的。波士顿房价这个数据集上条件数通常在千级到万级属于“还能接受”的范围但如果哪天你在自己的项目里发现条件数到了亿级别犹豫先去处理特征共线性再回来看模型结果。这个项目后续还有很多可以玩的地方换成学习曲线分析过拟合、加入交叉验证调 λ、尝试特征组合筛选。但主线是不变的——矩阵运算、解析解、数值稳定性这些底层逻辑才是真正能迁移到任何模型上的能力。
阅读完成 · 觉得有帮助？
