最小二乘法原理推导与实战：从正规方程到异常值处理

1. 最小二乘法到底在解决什么问题第一次接触最小二乘法很多人会被那一堆求和符号和矩阵求导吓退。但如果你把它放回它诞生的场景里其实特别朴素我们手里有一堆看起来杂乱的数据点想找一条“最合适”的线或者曲线去描述它们之间的关系。这条线不要求穿过每一个点那不现实数据本身有噪声它只要求整体上离所有点都“最近”。这个“最近”怎么定义就是最小二乘法的灵魂。你当然可以用“每个点到线的距离之和最小”来定义但绝对值在数学上不好求导处理起来很麻烦。所以高斯和勒让德当年选了一个更聪明的办法把每个点的误差先平方再求和让这个平方和最小。平方有两个好处一是把正负误差都变成正的不会互相抵消二是平方函数光滑可导求极值特别方便。这就是“最小二乘”这个名字的由来——最小化误差的平方和。它解决的问题非常广。你看到的任何一条趋势线、任何一个拟合出来的经验公式、任何一个从实验数据反推物理常数的过程背后大概率都是最小二乘法在工作。从天文观测里推算行星轨道到工程上标定传感器再到机器学习里最基础的线性回归全是它的地盘。适合谁来学只要你在跟数据打交道不管是做实验、写代码、搞金融分析还是做控制算法这一关都绕不过去。我见过太多人直接调numpy.polyfit或者sklearn.LinearRegression就把活干了结果数据里有个明显的异常点拟合出来的线偏得离谱还不知道问题出在哪。根子就在于没搞懂最小二乘法的原理和它的软肋。这篇文章我就按自己的理解把原理、推导、实操和踩过的坑一次讲透。2. 从几何直觉到代数表达最小二乘法的原理拆解2.1 用生活化的例子理解“最小误差平方和”假设你在做一个物理实验测弹簧的伸长量和所挂钩码重量的关系。理论上胡克定律告诉你 F kx是一条过原点的直线。但你测了五组数据画在坐标纸上发现点并不在一条直线上有的偏上有的偏下。现在你要从这五组数据里估出劲度系数 k。最直觉的做法是拿一把透明直尺在图上比划让尺子尽量贴近所有点然后读斜率。你脑子里其实就在做最小二乘法只不过你是用眼睛在最小化“点到尺子的垂直距离”。最小二乘法把这个过程数学化了它定义每个点的误差是“实际 y 值减去直线预测的 y 值”然后让所有误差的平方加起来最小。为什么不是让误差绝对值之和最小因为绝对值函数在零点不可导你没法用求导等于零这套流畅的工具去解。而平方和是个二次函数二次函数的极值有现成的公式求导令其为零就能得到解析解。这就是数学上的“好算”压倒了“直觉上更合理”。实际用下来平方和对大误差的惩罚更重所以拟合结果对异常值也更敏感这一点后面会专门讲。2.2 一元线性回归的公式推导全过程我们把场景简化到最经典的一元情况有一组观测数据 (x_i, y_i)i 从 1 到 n假设它们近似满足 y a b x。注意我这里用 a 表示截距b 表示斜率跟有些教材的记号可能相反但含义一样。每个点的残差记作 e_i y_i - (a b x_i)。最小二乘的目标就是让S(a, b) Σ e_i² Σ [y_i - (a b x_i)]²取最小值。这里 Σ 都是从 i1 到 n 求和。因为 S 是 a 和 b 的二次函数分别对 a 和 b 求偏导令偏导等于零就能得到两个方程。先对 a 求偏导 ∂S/∂a -2 Σ [y_i - (a b x_i)] 0化简得到 Σ y_i - n a - b Σ x_i 0 也就是 n a b Σ x_i Σ y_i …… (1)再对 b 求偏导 ∂S/∂b -2 Σ x_i [y_i - (a b x_i)] 0化简得到 Σ x_i y_i - a Σ x_i - b Σ x_i² 0 也就是 a Σ x_i b Σ x_i² Σ x_i y_i …… (2)这两个方程就是著名的正规方程组。两个未知数 a、b两个方程直接解就行。从 (1) 式解出 a a (Σ y_i - b Σ x_i) / n ȳ - b x̄其中 ȳ 是 y 的均值x̄ 是 x 的均值。这个结果很漂亮拟合直线一定穿过样本点的重心 (x̄, ȳ)。这是个很重要的几何性质后面检查拟合结果时可以用。把 a 代回 (2) 式经过一通整理中间过程就是移项和合并同类项可以得到 b 的表达式b [Σ (x_i - x̄)(y_i - ȳ)] / [Σ (x_i - x̄)²]这个公式特别好记分子是 x 和 y 的协方差未除以 n分母是 x 的方差未除以 n。所以斜率本质上就是“x 和 y 一起变化的程度”除以“x 自己变化的程度”。我建议你亲手把这个推导在纸上写一遍不要只看。写的过程中你会对“为什么最小二乘一定有解析解”有肌肉记忆。很多人卡在求和符号的运算上其实就是移项和展开的基本功慢一点都能推出来。2.3 从一元到多元矩阵形式的统一表达实际工作中一元的情况很少更多是多个自变量。比如房价预测影响因素有面积、房龄、楼层、地段等等。这时候写成矩阵形式会清爽很多。设设计矩阵 X 是 n 行 p 列n 个样本p 个特征通常第一列全是 1 用来吸收截距参数向量 β 是 p 维列向量观测值 y 是 n 维列向量。模型写作y ≈ X β残差向量 e y - X β。目标函数是S(β) eᵀe (y - X β)ᵀ (y - X β)展开 S(β) yᵀy - yᵀXβ - βᵀXᵀy βᵀXᵀXβ注意 yᵀXβ 和 βᵀXᵀy 都是标量且互为转置所以相等。于是 S(β) yᵀy - 2 βᵀXᵀy βᵀXᵀXβ对 β 求梯度矩阵求导的规则这里直接用不展开证明 ∂S/∂β -2 Xᵀy 2 XᵀXβ令其为零向量 XᵀXβ Xᵀy这就是矩阵形式的正规方程。如果 XᵀX 可逆解就是β (XᵀX)⁻¹ Xᵀy这个公式是整个最小二乘法的核心。你看到的所有线性回归的解析解都是它。理解了这个你就理解了为什么线性回归训练那么快——它根本不需要迭代一步矩阵运算就出结果。注意XᵀX 可逆的前提是特征之间不能有完全的线性相关也就是不能有某一列是其他列的线性组合。实际数据里如果两个特征高度相关比如“面积”和“房间数”XᵀX 会接近奇异求逆数值不稳定这时候就要用正则化或者剔除特征。3. 实操落地从零手写最小二乘拟合3.1 用 Python 手算一遍正规方程光看公式没用我带你用代码走一遍。先造一组带噪声的数据然后用正规方程解出参数再跟库函数的结果对比。import numpy as np # 造数据真实关系 y 3 2x加一点高斯噪声 np.random.seed(42) n 50 x np.linspace(0, 10, n) y_true 3 2 * x noise np.random.normal(0, 1.5, n) y y_true noise # 构造设计矩阵第一列全1用于截距 X np.column_stack([np.ones(n), x]) # 正规方程求解 beta (X^T X)^-1 X^T y XTX X.T X XTy X.T y beta np.linalg.inv(XTX) XTy print(手算结果截距 %.4f, 斜率 %.4f % (beta[0], beta[1])) # 用 numpy 内置的最小二乘对比 beta_np, residuals, rank, sv np.linalg.lstsq(X, y, rcondNone) print(lstsq结果截距 %.4f, 斜率 %.4f % (beta_np[0], beta_np[1]))跑下来你会发现两者几乎一模一样差异只在浮点误差级别。真实参数是 3 和 2拟合出来大概在 3.1 和 1.98 附近取决于噪声。这就是最小二乘的威力即使数据有噪声它也能把真实参数估得八九不离十。这里有个细节值得说我用了np.linalg.inv直接求逆这在教学演示里没问题但实际工程中不要显式求逆。原因有两个一是求逆计算量大二是数值稳定性差。正确做法是用np.linalg.solve(XTX, XTy)或者直接用lstsq它们内部用 QR 分解或 SVD更稳。3.2 参数求解的数值稳定性处理上面提到不要显式求逆这里展开说一下为什么。XᵀX 的条件数是 X 的条件数的平方。如果 X 本身条件数就大比如特征尺度差异巨大平方之后条件数会爆炸求逆结果误差极大。举个例子假设你有两个特征一个是“面积”单位平方米数值在 50 到 200 之间另一个是“房龄”单位年数值在 0 到 30 之间。这两个特征尺度差了一个数量级XᵀX 的条件数就会比较大。解决办法是标准化把每个特征减去均值再除以标准差让它们都在同一量级。标准化之后求出来的系数是相对于标准化数据的要还原回原始尺度需要做相应变换。我自己的习惯是只要特征超过两个一律先标准化再拟合。多花两行代码省去后面调参的无数麻烦。另外如果特征数量 p 接近甚至超过样本数 nXᵀX 直接不可逆这时候必须上正则化也就是岭回归相当于在 XᵀX 对角线上加一个小的正数 λI强行让它可逆。3.3 拟合效果的评估与可视化拟合完了不能只看系数得评估效果。最常用的指标是决定系数 R²定义是R² 1 - SS_res / SS_tot其中 SS_res Σ (y_i - ŷ_i)² 是残差平方和SS_tot Σ (y_i - ȳ)² 是总平方和。R² 越接近 1 说明拟合越好等于 1 表示完美拟合等于 0 表示跟直接用均值预测一样烂。注意 R² 可以是负数说明模型比瞎猜均值还差。y_pred X beta ss_res np.sum((y - y_pred) ** 2) ss_tot np.sum((y - np.mean(y)) ** 2) r2 1 - ss_res / ss_tot print(R² %.4f % r2)可视化也很重要。把原始散点和拟合直线画在一起一眼就能看出有没有系统性的偏差。如果残差图呈现出喇叭形或者弯曲的规律说明线性模型不合适可能需要对 y 做变换或者引入 x 的平方项。实操心得残差图比 R² 更能说明问题。R² 高不代表模型对可能只是数据本身方差大。我习惯把残差对预测值画散点图理想情况是随机分布的一条水平带如果出现任何规律都说明模型漏掉了某些结构。4. 常见问题与排查技巧实录4.1 异常值为什么能把拟合带偏最小二乘对异常值极其敏感这是它平方误差定义的直接后果。一个偏离很远的点它的误差平方会非常大为了最小化总和拟合线会向这个点“妥协”整体偏移。我做过一个测试在 50 个正常点里插入一个 y 值偏了 10 倍的点拟合斜率直接变了 30%。而如果改用最小绝对偏差L1 损失那个异常点的影响就小得多。所以如果你的数据里可能有录入错误或者传感器偶发故障一定要先做异常值检测。常用办法是看标准化残差绝对值大于 3 的点标记为可疑然后人工核查。也可以用箱线图看 y 的分布超出 1.5 倍四分位距的点先拎出来看看。不要无脑删但一定要查。4.2 多重共线性导致系数符号反常多元最小二乘里一个经典坑是多重共线性。两个或多个特征高度相关时XᵀX 接近奇异解出来的系数会非常大而且符号可能跟常识相反。比如你用“总价”和“单价”同时预测“面积”这两个特征几乎线性相关拟合出来的系数会乱跳。判断方法是算方差膨胀因子 VIF如果某个特征的 VIF 大于 10就说明共线性严重。解决办法有几种删掉相关性高的特征之一、用主成分分析降维、或者上岭回归。我一般优先考虑删特征因为解释起来最直接。4.3 非线性关系硬套线性模型的后果最小二乘本身不要求模型一定是线性的它要求的是“参数线性”。也就是说 y a b x² 也是线性模型因为对 a 和 b 是线性的。但如果你硬用 y a b x 去拟合一个明显弯曲的关系残差就会有系统性规律。排查方法还是看残差图。如果残差呈现 U 形或者倒 U 形说明漏了二次项。解决办法很简单加一列 x² 进设计矩阵重新拟合。这就是多项式回归的本质依然是最小二乘只是特征变了。下面这张表是我整理的最小二乘常见问题速查问题现象可能原因排查方法解决手段拟合线明显偏向某几个点存在异常值标准化残差绝对值大于3核查后剔除或改用稳健回归系数符号与常识相反多重共线性计算VIF删特征、降维或岭回归残差图呈U形模型欠拟合漏了非线性项残差对预测值散点图加入x²等高阶项R²很高但预测新数据很差过拟合交叉验证减少特征或加正则化求逆报奇异矩阵错误特征完全线性相关或p大于n检查特征相关性删冗余特征或用伪逆4.4 数值计算中的精度陷阱最后说一个容易被忽略的点数值精度。当 x 的数值很大时比如年份 2020、2021 这种Σ x² 会非常大而 Σ x 和 n 也不小正规方程里两个方程的数量级差异巨大直接求解会有精度损失。解决办法是中心化把 x 减去均值再拟合。这样 Σ (x - x̄) 0正规方程的第一个方程直接给出 a ȳ第二个方程只涉及 (x - x̄)数值小很多精度大幅提升。这个技巧在多项式拟合里尤其重要因为 x 的高次幂会让数值爆炸。我见过有人用原始年份做五次多项式拟合结果完全跑飞中心化之后立刻正常。提示中心化不改变斜率只改变截距的含义。拟合完把截距换算回原始尺度即可换算关系是 a_原始 a_中心化 - b * x̄。5. 最小二乘法的扩展与个人经验5.1 加权最小二乘当数据可信度不同时普通最小二乘假设每个样本的误差方差相同但实际中往往不是。比如你测某个量高浓度区域的测量误差比低浓度区域大这时候应该给可信度高的点更大的权重。加权最小二乘的目标函数变成 S Σ w_i e_i²w_i 是权重。推导过程跟普通最小二乘几乎一样只是每个求和项多了个权重。矩阵形式的解变成β (XᵀWX)⁻¹ XᵀWy其中 W 是对角矩阵对角线元素是 w_i。权重怎么定如果知道每个测量的方差 σ_i²通常取 w_i 1/σ_i²这是理论上最优的。不知道的话可以根据经验给比如测量次数多的点权重大。5.2 递推最小二乘数据流式到达时的处理有时候数据不是一次性给你的而是一个一个来。比如在线系统里每来一个新样本就要更新模型参数。这时候如果每次都用全部历史数据重新算一遍计算量会越来越大。递推最小二乘解决了这个问题。它维护当前的参数估计和协方差矩阵每来一个新样本用一组简单的更新公式修正参数不需要保留历史数据。核心公式是K P φ / (λ φᵀ P φ) β_new β_old K (y - φᵀ β_old) P_new (P - K φᵀ P) / λ其中 φ 是新样本的特征向量P 是协方差矩阵λ 是遗忘因子通常取 0.95 到 1 之间用来让旧数据的影响逐渐衰减。这套公式在自适应滤波和控制领域用得非常多我当年做传感器标定的时候就是靠它实现在线校准的。5.3 我踩过的几个坑和对应建议第一个坑是忘记加截距列。有次我直接用 X 拟合忘了第一列全 1结果拟合线强行过原点偏差巨大。后来养成习惯构造设计矩阵后先打印形状和第一列确认。第二个坑是特征尺度差异太大导致收敛慢。虽然最小二乘有解析解不需要迭代但如果你用梯度下降去优化尺度不统一会让收敛慢得让人崩溃。即使是用解析解标准化也能提升数值稳定性。所以不管怎样标准化都是好习惯。第三个坑是把相关当因果。最小二乘只能告诉你变量之间有线性关联不能证明因果。我见过有人用最小二乘拟合出“冰淇淋销量和溺水人数高度相关”然后得出吃冰淇淋导致溺水的荒谬结论。背后其实是气温这个共同因素在起作用。做数据分析时因果推断需要额外的实验设计或工具变量方法不能只靠拟合。5.4 后续可以深入的方向如果你已经把基础最小二乘玩熟了可以往这几个方向走。一是正则化方法岭回归和 Lasso解决过拟合和共线性问题。二是广义最小二乘处理误差项存在相关性的情况。三是非线性最小二乘比如高斯-牛顿法和列文伯格-马夸尔特法用于拟合非线性模型。四是贝叶斯线性回归把参数看成随机变量给出后验分布而不是点估计能更好地量化不确定性。这些方向每一个都够写一篇长文但根都在最小二乘这里。把今天讲的原理和推导吃透后面学这些会顺很多。我自己是从手推正规方程开始到用代码实现再到踩坑排查一步步走过来的。最深的体会是公式要亲手推代码要亲手写坑要亲手踩这三样缺一不可。看别人的推导觉得懂了自己一写就卡壳这是常态。慢下来把每一步的来龙去脉搞清楚比快速刷完十篇教程有用得多。