ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

最小二乘法从原理到实战:手推正规方程与代码实现

最小二乘法从原理到实战:手推正规方程与代码实现 最小二乘法这个名字很多人在大学课堂上第一次听到时觉得它不过就是个拟合工具工作几年之后才发现它几乎是所有数据类项目的底层地基。我做数据建模和算法落地这些年从传感器标定、广告归因到推荐系统的特征拟合绕来绕去总会回到这条“让误差平方和最小”的主线上。它解决的问题很朴素手上有一堆带噪声的观测点想找一条或一个最靠谱的规律去描述它们同时还要能说清楚这条规律到底靠不靠谱。这篇文章我打算把最小二乘法从直觉、原理、公式推导到代码实现完整走一遍适合刚接触拟合的初学者也适合想重新把数学底子捡回来的工程师。看完你至少能做到三件事自己推导出正规方程、手写一个可用的最小二乘拟合函数、以及在遇到病态矩阵时知道该往哪个方向排查。1. 最小二乘法到底在解决什么问题1.1 从一个生活场景说起假设你开了个奶茶店想研究“气温”和“冰饮销量”之间的关系。你连续记录了十天每天的最高气温和当天冰饮卖出的杯数。把这些点画在坐标纸上你会发现它们大致沿着一条斜向上的直线分布但又不完全在一条线上——因为销量还受周末、促销、隔壁新店开业等因素影响。这时候你想画一条“最能代表整体趋势”的直线问题就来了能穿过所有点的直线几乎不存在那到底哪条线才算“最好”最小二乘法给出的答案非常直接对每一条候选直线计算每个真实点到这条直线的垂直偏差残差把这些偏差平方之后全部加起来谁的平方和最小谁就是最优的那条线。为什么是平方而不是绝对值这个问题我后面会专门讲因为它牵扯到最小二乘法最核心的数学动机。用一句话概括最小二乘法是一种通过最小化残差平方和来寻找最优拟合参数的数学方法。它既能做直线拟合也能做多项式拟合、多元回归甚至是复杂的非线性问题的线性化求解。1.2 为什么是“二乘”而不是“一次”或“三次”“二乘”在中文里就是平方的意思。选择平方而不是绝对值背后有三个非常实际的考量这也是我在实际项目中反复验证过的可导性平方函数处处可导而绝对值函数在零点不可导。可导意味着我们可以用求导等于零的方式直接解出解析解这是最小二乘法能推出漂亮闭式公式的根本原因。对大误差更敏感平方会放大较大的偏差抑制较小的偏差。换句话说它更“讨厌”离谱的离群点倾向于让整体拟合不要被个别极端值带偏太多。统计意义当观测误差服从正态分布时最小二乘估计等价于极大似然估计。这一点非常关键它让最小二乘法不只是一个几何技巧而是有严格统计理论支撑的估计方法。注意平方对离群点敏感是双刃剑。如果数据里混入了明显的错误点普通最小二乘会被严重拉偏这时候要考虑鲁棒回归或者先做异常值剔除。1.3 最小二乘法的适用边界不是所有拟合问题都适合直接上最小二乘。我总结了几条判断标准场景特征是否适合普通最小二乘说明误差主要在因变量方向适合经典假设自变量测量精确自变量和因变量都有明显误差需谨慎应考虑整体最小二乘或正交回归存在强离群点不适合改用Huber、RANSAC等鲁棒方法特征数远大于样本数不适合直接求解需引入正则化如岭回归关系明显非线性需变换可通过对数、多项式扩展转为线性这张表是我踩过坑之后整理的。早期做传感器标定时我一度忽略了自变量本身也有噪声结果拟合出来的系数系统性偏小后来换成整体最小二乘才把问题解决。2. 原理详解从几何直觉到统计基础2.1 几何视角投影与正交理解最小二乘法最漂亮的方式是几何。把观测值向量 y 想象成空间里的一个点把模型能表示的所有可能结果构成一个子空间列空间。最小二乘要做的就是在这个子空间里找一个离 y 最近的点。而空间中点到子空间的最短距离就是垂线——也就是投影。这个结论直接推出一个重要性质残差向量与拟合值向量正交。用公式写就是 X 转置乘以残差等于零这正是正规方程的来源。我第一次真正理解这个正交性的时候感觉整个最小二乘的推导都变得顺理成章了因为它不再是一堆代数变形而是一个清晰的几何事实。2.2 代数视角目标函数与求导设模型为 y Xβ ε其中 X 是设计矩阵β 是待求参数ε 是误差。残差平方和定义为S(β) (y - Xβ)ᵀ(y - Xβ)展开后得到S(β) yᵀy - 2βᵀXᵀy βᵀXᵀXβ对 β 求梯度并令其为零∂S/∂β -2Xᵀy 2XᵀXβ 0整理得到正规方程XᵀXβ Xᵀy如果 XᵀX 可逆解就是β (XᵀX)⁻¹Xᵀy这就是最小二乘法最核心的公式。整个推导过程只用了矩阵求导和二次函数极值两个工具非常干净。2.3 统计视角为什么它是最优的从统计角度看如果误差 ε 满足均值为零、方差为 σ²、且相互独立那么最小二乘估计 β 具有几个优良性质无偏性E(β) 真实参数估计不会系统性偏高或偏低。最小方差在所有线性无偏估计中最小二乘估计的方差最小这就是著名的高斯-马尔可夫定理。与极大似然等价当误差服从正态分布时最小化平方和等价于最大化似然函数。这三个性质解释了为什么最小二乘法能成为工程和统计领域的默认选择。它不是随便挑的一个损失函数而是在一组合理假设下被数学证明为最优的方案。2.4 与极大似然估计的推导联系很多人好奇最小二乘和正态分布到底怎么扯上关系。推导其实很直接。假设每个观测误差独立且服从 N(0, σ²)那么似然函数为L(β) ∏ (1/(√(2π)σ)) exp(-(yᵢ - xᵢβ)²/(2σ²))取对数后ln L -n/2 ln(2πσ²) - (1/(2σ²)) Σ(yᵢ - xᵢβ)²要让似然最大就要让最后那个求和项最小而它正好就是残差平方和。所以在正态误差假设下最小二乘就是极大似然。这也解释了为什么最小二乘对离群点敏感——正态分布的尾部很薄它天然不认为会出现极端偏差所以一旦出现就会被过度惩罚。3. 公式推导一步步手推正规方程3.1 一元线性回归的完整推导先从最简单的一元情况入手模型是 y a bx。残差平方和S(a, b) Σ(yᵢ - a - bxᵢ)²分别对 a 和 b 求偏导并令其为零∂S/∂a -2Σ(yᵢ - a - bxᵢ) 0∂S/∂b -2Σxᵢ(yᵢ - a - bxᵢ) 0由第一个式子得到a ȳ - b x̄把它代入第二个式子经过整理可以得到斜率b Σ(xᵢ - x̄)(yᵢ - ȳ) / Σ(xᵢ - x̄)²这个公式非常实用分子是协方差的未归一化形式分母是 x 的方差未归一化形式。我在做快速估算时经常直接用它不需要构建矩阵。3.2 多元线性回归的矩阵推导推广到多元情况模型是 y Xβ。前面已经推过核心结果是β (XᵀX)⁻¹Xᵀy这里要特别注意 X 的构造。如果模型带截距项X 的第一列必须全是 1。我见过不少初学者忘记加这一列结果拟合出来的线强行过原点误差大得离谱。设计矩阵 X 的结构如下样本截距列特征1特征2...特征p11x₁₁x₁₂...x₁ₚ21x₂₁x₂₂...x₂ₚ..................3.3 正规方程的数值陷阱理论上 β (XᵀX)⁻¹Xᵀy 完美但直接求逆在数值上非常危险。原因有两个条件数放大XᵀX 的条件数是 X 条件数的平方本来就不良的问题会变得更糟。显式求逆不稳定矩阵求逆的浮点误差在病态情况下会被急剧放大。实际工程中我从来不直接求逆而是用以下方法之一QR分解把 X 分解为正交矩阵 Q 和上三角矩阵 R然后解 Rβ Qᵀy。数值稳定性最好。SVD分解奇异值分解对病态问题最鲁棒还能顺便给出伪逆解。Cholesky分解当 XᵀX 正定时可用速度比QR快但稳定性略差。提示如果你用 numpy直接调用 np.linalg.lstsq 就行它内部用的是SVD稳定性和精度都有保障不要自己写求逆。3.4 从平方和到卡方分布的延伸在统计推断里残差平方和除以真实方差后服从卡方分布自由度是 n - p样本数减参数个数。这个结论是构造置信区间和做假设检验的基础。推导思路是残差是观测值在正交补空间上的投影而正态向量在正交变换下仍然是正态向量其平方和自然服从卡方分布。理解了这一层你就能明白为什么回归分析里总要看自由度以及为什么参数越多残差平方和的期望越小。4. 实操过程从零实现最小二乘拟合4.1 环境准备与数据构造我用 Python 演示依赖只有 numpy 和 matplotlib。先构造一组带噪声的数据import numpy as np import matplotlib.pyplot as plt np.random.seed(42) n 50 x np.linspace(0, 10, n) true_a, true_b 2.5, 1.8 noise np.random.normal(0, 1.5, n) y true_a true_b * x noise这里真实截距是2.5真实斜率是1.8噪声标准差1.5。构造数据时固定随机种子方便复现结果这是做实验的基本习惯。4.2 手写正规方程求解先按公式手写一遍理解每一步X np.column_stack([np.ones(n), x]) XtX X.T X Xty X.T y beta np.linalg.solve(XtX, Xty) print(截距:, beta[0], 斜率:, beta[1])注意我用的是 np.linalg.solve 而不是求逆这样数值上更稳。运行后你会看到截距接近2.5斜率接近1.8说明拟合成功。4.3 用lstsq做对比验证再用官方稳定实现跑一遍beta_lstsq, residuals, rank, sv np.linalg.lstsq(X, y, rcondNone) print(lstsq截距:, beta_lstsq[0], 斜率:, beta_lstsq[1])两次结果应该几乎一致。如果差异明显说明你的矩阵条件数可能有问题需要检查特征是否共线。4.4 拟合效果可视化与残差分析y_pred X beta plt.scatter(x, y, label观测点) plt.plot(x, y_pred, colorred, label拟合线) plt.legend() plt.show() residuals y - y_pred print(残差均值:, residuals.mean()) print(残差标准差:, residuals.std())残差均值应该接近零残差图不应该呈现明显的喇叭形或弯曲趋势。如果残差随 x 增大而扩散说明存在异方差普通最小二乘的标准误估计会失真。4.5 多项式拟合的扩展最小二乘不只做直线。把 X 扩展成多项式特征即可degree 3 X_poly np.column_stack([x**i for i in range(degree 1)]) beta_poly np.linalg.lstsq(X_poly, y, rcondNone)[0]但要注意高阶多项式会导致特征之间高度相关XᵀX 条件数急剧上升。我一般会先对 x 做标准化或者改用正交多项式基否则数值结果会很难看。5. 常见问题与排查技巧实录5.1 拟合结果完全不对怎么排查这是新手最常遇到的问题。我整理了一个排查顺序表现象可能原因排查方法拟合线过原点忘记加截距列检查X第一列是否全为1系数巨大且不稳定特征共线或未标准化计算条件数做标准化残差有明显规律模型形式不对画残差图考虑非线性项结果每次运行都变数据未固定随机种子设置seed预测值偏离极远外推超出训练范围限制预测区间5.2 多重共线性的识别与处理多重共线性是最小二乘最隐蔽的敌人。识别方法是计算方差膨胀因子VIF如果某个特征的VIF超过10说明它和其他特征高度相关。处理方法有删除冗余特征使用主成分回归引入岭回归在正规方程里加一个 λI 项岭回归的解是 β (XᵀX λI)⁻¹Xᵀy这个小小的 λI 能让矩阵从接近奇异变得良态代价是引入一点偏差。我在特征多、样本少的项目里几乎默认用岭回归。5.3 数值精度问题的实战经验有一次我处理一组量纲差异极大的数据一个特征范围是0到1另一个是0到一百万结果拟合系数完全乱套。原因是 XᵀX 的条件数被量纲差异撑到了天文数字。解决办法很简单标准化。把每个特征减去均值再除以标准差拟合完再把系数变换回原始尺度。这个操作我后来做成了标准流程几乎每个项目都会先跑一遍。5.4 加权最小二乘的使用时机当不同观测点的可信度不同时普通最小二乘一视同仁就不合理了。加权最小二乘给每个点一个权重 wᵢ目标函数变成 Σwᵢ(yᵢ - xᵢβ)²。解的形式是β (XᵀWX)⁻¹XᵀWy其中 W 是对角权重矩阵。权重通常取观测方差的倒数方差越大的点权重越小。这个技巧在处理传感器融合数据时特别有用因为不同传感器的精度本来就不一样。5.5 最小二乘与RANSAC的配合数据里混入少量严重离群点时最小二乘会被带偏。我的做法是先用RANSAC随机采样出一批内点再用最小二乘在这批内点上做精细拟合。RANSAC负责抗噪最小二乘负责精度两者配合效果很好。这个组合在计算机视觉的直线检测、平面拟合里是标准套路。6. 影响范围与延伸应用6.1 在机器学习中的位置线性回归就是最小二乘的直接应用。逻辑回归虽然用的是交叉熵损失但在某些近似下也和最小二乘有联系。神经网络训练用的梯度下降本质上是在做非线性最小二乘的迭代求解。理解了最小二乘你就理解了损失函数、优化、正则化这一整条线的起点。6.2 在信号处理与控制系统中的应用系统辨识里最小二乘用来估计传递函数参数。卡尔曼滤波的更新步骤在特定条件下也等价于递推最小二乘。我在做传感器标定时用最小二乘拟合温度补偿曲线把测量误差从百分之几降到了千分之几效果立竿见影。6.3 在计算机视觉中的延伸相机标定、单应矩阵估计、光流计算背后都有最小二乘的影子。比如拟合椭圆、直线检测都是先建立残差方程再用最小二乘求解。当问题变成非线性时就用高斯-牛顿或列文伯格-马夸尔特方法迭代求解而每一步迭代内部仍然是一个线性最小二乘问题。6.4 从最小二乘到正则化的演进普通最小二乘在特征多的时候会过拟合。岭回归加了L2正则Lasso加了L1正则弹性网把两者结合。这些方法的共同点是在残差平方和后面加一个惩罚项本质上是给最小二乘加了一个先验约束。理解了最小二乘这些正则化方法就都是自然的延伸而不是孤立的技巧。6.5 暴力枚举与解析解的对比思考有人会问既然参数空间不大为什么不直接暴力枚举所有可能的参数组合选平方和最小的那个理论上可行但维度一高就彻底不可行。假设有10个参数每个参数枚举100个值总组合是100的10次方天文数字。而解析解一步到位这就是数学推导的价值。暴力枚举可以作为验证手段用来检查解析解是否正确但绝不能作为主算法。7. 实操心得与避坑清单7.1 我踩过的三个典型坑第一个坑是忘记加截距列导致拟合线强行过原点误差大得离谱排查了半天才发现是设计矩阵的问题。第二个坑是特征未标准化量纲差异导致系数完全不可解释后来养成习惯拟合前先做标准化。第三个坑是直接用求逆解正规方程遇到病态矩阵时结果完全不可信换成lstsq之后问题消失。这三个坑现在都写进了我的检查清单。7.2 参数选择的经验法则样本数至少是参数个数的10倍否则过拟合风险很高。多项式阶数不要超过5再高就该考虑其他模型形式了。条件数超过1000就要警惕超过10000基本可以判定病态。残差标准差应该和观测噪声量级相当差太多说明模型有问题。7.3 验证拟合质量的实用方法除了看R²我更看重残差图。R²高不代表模型对残差图有规律就说明模型漏掉了某些结构。另外做交叉验证把数据分成训练集和测试集看测试集上的预测误差这比单看训练集R²可靠得多。我一般还会留一部分数据做外推测试看看模型在训练范围之外的稳定性。7.4 代码复用的封装建议把最小二乘拟合封装成一个函数输入X和y输出系数、残差、R²、条件数。这样每次用的时候直接调用不用重复写。我自己的工具函数里还加了自动标准化和自动加截距列的选项减少人为失误。封装的时候记得把中间量也返回出来方便排查问题。7.5 什么时候不该用最小二乘如果误差不满足独立同分布假设或者数据里离群点比例超过10%或者特征严重共线且无法通过正则化解决那就该考虑其他方法了。最小二乘是强大的默认选项但不是万能钥匙。判断标准很简单先跑一遍最小二乘看残差图如果残差有结构或者被离群点主导就换方法。最后分享一个我常用的快速验证技巧拿到任何一组数据先画散点图再跑最小二乘再画残差图这三步走完数据的基本结构和你模型的问题基本就暴露得差不多了。这个流程我用了很多年几乎没失手过。
返回列表