我第一次真正用上最小二乘法是在处理一组传感器校准实验数据的时候。当时手里是一堆带噪声的采样点教科书里的公式能帮我“找一条最合理的直线”但工程上真正让人头疼的其实是后面这些事数据里混了离群点怎么办多项式阶数选多少才不会过拟合拟合出的结果到底可不可信。这篇文章想做的就是把最小二乘法从数学原理讲到代码实现再用一个带噪声的拟合案例从头到尾串一遍。适合刚接触数值计算、需要做回归拟合或者想看明白各种回归库背后在做什么的开发者。最小二乘法并不是什么高深技巧它更像是一套统一的“找最好参数”的思路。只要你手里有一组观测数据想用一个模型去描述它背后的规律最小二乘就是最常用的切入方式。下面我会从直觉出发把公式一步步拆开再给出可直接运行的代码最后聊几个我自己踩过的坑。1. 最小二乘法的核心思想从误差平方说起1.1 什么样的拟合才算“最好”先回到最简单的场景。给你一堆平面上的点大概分布在一条直线附近你要画一条直线穿过它们。问题来了什么叫“最合理”每个人都能画一条线但哪条线才是公认最好的最自然的想法是让这条直线尽可能靠近每一个点。也就是说对每一个数据点算出我预测的位置和实际位置之间的距离把这些距离加起来和最小的那条线就是最好的。这里的“距离”在统计里叫“残差”也就是观测值和模型预测值之差。最小二乘法的名字就是这么来的让残差的平方和最小。[ \min_{\beta_0,\beta_1} \sum_{i1}^{n} \left( y_i - (\beta_0 \beta_1 x_i) \right)^2 ]这个式子看着简单但它背后藏着两个关键决策一是为什么要用“平方”二是为什么求和而不是求平均。后一点比较好理解求和和求平均在参数优化上是等价的除以样本量也不会改变最优解的位置。真正需要琢磨的是“平方”。1.2 为什么选平方而不是绝对值你可能第一反应是既然是距离直接用残差的绝对值 ( |y_i - \hat{y}_i| ) 求和不就行了吗绝对值更好理解正负误差也不会抵消。确实把绝对值误差和最小化也是一种回归方法叫最小绝对偏差回归但它有几个工程上的麻烦。首先绝对值函数在零点不可导。最小二乘可以靠求导解出显式公式绝对值不行优化起来费劲。其次平方会让大误差被放大。如果一个点偏离模型特别远它的误差平方会以一个平方级的速度把损失“顶”上去这样优化算法会优先照顾那些偏差大的点结果对异常值特别敏感。说得更直白一点最小二乘等价于假设误差是高斯分布下的最大似然估计平方误差和最小化天然匹配“误差独立且服从正态分布”的场景。生活里其实也有这个直觉。你想用几个测量值估计一个物理量的真实值最简单的方法是取算术平均。为什么是平均而不是取中位数因为在平方误差的意义下均值恰好就是让所有点到它的平方距离之和最小的那个数。所以最小二乘并不只是一个花哨的数学技巧它和我们日常处理测量数据的方式是相通的。1.3 最小二乘不只是直线拟合还是一套基函数框架很多刚接触的人以为最小二乘只能拟合直线这是最大的误解。它真正的威力在于只要模型对参数是线性的就能用同一套方法求解。比如下面这些模型对参数 (\beta) 都是线性的[ y \beta_0 \beta_1 x \beta_2 x^2 ][ y \beta_0 \beta_1 \sin(x) \beta_2 \cos(x) ][ y \beta_0 \beta_1 x_1 \beta_2 x_2 \beta_3 x_1 x_2 ]这些模型里(x) 和 (x^2)、 (\sin(x)) 这些都是已知的“特征”模型对未知参数 (\beta) 仍然保持线性。所以你可以把多项式拟合、多元线性回归、三角基函数拟合全部塞进最小二乘的框架里。真正需要小心的是那些参数本身不是线性的模型比如 (y \beta_0 e^{\beta_1 x})它的指数里带了 (\beta_1)这才是真正的非线性问题得用非线性最小二乘或专门优化算法来处理。理解这一点非常关键因为它决定了你写代码时的思路先构造一个设计矩阵把特征都算好然后交给同一个求解器。2. 正规方程推导手把手写出最小二乘解2.1 一元线性回归的闭式解怎么来的先从一个特征的情况入手。假设模型是 (y \beta_0 \beta_1 x)我们有 (n) 组观测 ((x_i, y_i))。损失函数就是前面写的平方和[ L(\beta_0, \beta_1) \sum_{i1}^{n} (y_i - \beta_0 - \beta_1 x_i)^2 ]为了找到最小值我们让损失函数对 (\beta_0) 和 (\beta_1) 的偏导都等于 0。这一步是高等数学里很标准的操作极值点处的梯度为零。对 (\beta_0) 求偏导[ \frac{\partial L}{\partial \beta_0} -2 \sum_{i1}^{n} (y_i - \beta_0 - \beta_1 x_i) 0 ]对 (\beta_1) 求偏导[ \frac{\partial L}{\partial \beta_1} -2 \sum_{i1}^{n} x_i (y_i - \beta_0 - \beta_1 x_i) 0 ]整理之后可以得到我们都很熟悉的闭式解[ \beta_1 \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2} ][ \beta_0 \bar{y} - \beta_1 \bar{x} ]也就是说斜率等于 (x) 和 (y) 的协方差除以 (x) 的方差截距通过均值点反推。这个公式的好处是一目了然坏处是当特征变多时手推每个变量的偏导会越来越麻烦。所以实际工程里我们更习惯用矩阵形式。2.2 从矩阵形式到正规方程当特征有 (p) 个时把截距项并进参数向量设计矩阵 (X) 的第一列全部设为 1。模型写成[ y X\beta \varepsilon ]其中 (X) 是 (n \times (p1)) 矩阵(\beta) 是 (p1) 维列向量(y) 是 (n) 维列向量。目标是最小化[ L(\beta) | y - X\beta |^2 (y - X\beta)^T (y - X\beta) ]展开后得到[ L(\beta) y^T y - 2\beta^T X^T y \beta^T X^T X \beta ]对 (\beta) 求导并令其为零注意矩阵求导的规则就会有[ -2X^T y 2X^T X \beta 0 ]于是[ X^T X \beta X^T y ][ \beta (X^T X)^{-1} X^T y ]这个式子就是最小二乘的“正规方程”。只要 (X^T X) 可逆最小二乘解可以直接一步算出来不需要像梯度下降那样迭代。这也是为什么线性最小二乘在数据量不是特别大的时候速度快得很。补充一点上面的推导里我默认了矩阵求导的结果如果你手边没有现成的矩阵求导表可以直接把 (|y - X\beta|^2) 展开成求和形式然后一项一项对每个 (\beta_j) 求偏导得出的结果和矩阵形式完全一致。实际推导一遍之后你对正规方程的印象会深刻很多。2.3 正规方程在什么情况下会失效正规方程看着简洁实际应用时有几个致命问题。第一如果 (X^T X) 是奇异矩阵或者条件数非常大求解结果会很糟糕。这种情况通常发生在特征之间存在高度相关性时比如一个特征可以通过其他特征线性组合出来。这时 (X^T X) 的行列式接近零求逆会变得极不稳定微小的数据扰动就能让参数产生巨大波动。第二复杂度问题。((X^T X)^{-1}) 的计算复杂度大概是 (O(p^3))(p) 是特征数量。如果你有几百个特征还好几千几万个特征时求逆会非常慢。这也是为什么高维场景下更推荐用梯度下降或迭代方法而不是直接求解析解。第三数值稳定性。正规方程需要显式计算 (X^T X) 再求逆这中间会引入额外的数值误差。如果数据量不大且条件数正常直接用没问题但如果数据量较大或存在共线性我更推荐用奇异值分解或者 (np.linalg.lstsq) 去解最小二乘问题。(lstsq) 内部走的是数值更稳定的路线不显式构造 (X^T X)通常比先乘再求逆靠谱得多。3. 代码实现三种最常用的最小二乘写法3.1 numpy手写正规方程求解先看最“裸”的写法用 numpy 构造设计矩阵然后直接用正规方程求解。这里刻意用一个小数据集方便你观察每个参数的含义。import numpy as np x np.array([1.0, 2.0, 3.0, 4.0, 5.0]) y np.array([2.1, 3.9, 6.2, 8.1, 9.8]) # 构造设计矩阵第一列是1第二列是x X np.vstack([np.ones_like(x), x]).T # 正规方程解析解 theta np.linalg.inv(X.T X) X.T y print(截距和斜率, theta) # 结果约 [0.14, 1.96]用这组数据算出来的结果是截距约 0.14斜率约 1.96。也就是说拟合线大约是 (y 0.14 1.96x)。你可以用这个结果回去验算一下每个点的预测值和原始 (y) 对比会发现它们非常接近。不过实际工程里我不会直接用np.linalg.inv(X.T X)因为这种做法数值稳定性差。更推荐用np.linalg.lstsqtheta, residuals, rank, s np.linalg.lstsq(X, y, rcondNone) print(截距和斜率, theta)lstsq返回四个值第一个是最小二乘解第二个是残差平方和第三个是矩阵的秩第四个是奇异值。用它的好处是遇到近似奇异的矩阵也不会直接跑出无穷大而是会给出一个在数值上能接受的解。3.2 用多项式特征拟合带弯曲的数据如果数据分布不是直线而是明显弯曲的就可以在特征里加入 (x^2)、(x^3) 等多项式项。这里的关键不是模型复杂度而是“把非线性特征当成普通输入变量”。import numpy as np import matplotlib.pyplot as plt rng np.random.default_rng(42) x np.linspace(0, 2 * np.pi, 100) y_true 1.2 0.8 * np.sin(x) y y_true rng.normal(0, 0.15, x.size) # 构造多项式特征矩阵1, x, x^2, x^3 X np.column_stack([np.ones_like(x), x, x**2, x**3]) # 最小二乘求解 theta, _, _, _ np.linalg.lstsq(X, y, rcondNone) print(系数, theta) # 预测并可视化 x_smooth np.linspace(0, 2 * np.pi, 300) X_smooth np.column_stack([np.ones_like(x_smooth), x_smooth, x_smooth**2, x_smooth**3]) y_pred X_smooth theta plt.scatter(x, y, s12, alpha0.7, label带噪声数据) plt.plot(x_smooth, y_pred, colorred, label四次特征组合拟合) plt.legend() plt.show()这里用的是三次多项式所以 (X) 有四列常数项、一次项、二次项、三次项。最小二乘会自动估计出每个特征的权重。如果你的数据有周期性三次多项式通常拟合得不太好但这并不妨碍你理解流程只要构造了合适的特征矩阵最小二乘的求解逻辑没有任何变化。多项式拟合还有一个更偷懒的函数np.polyfit。它是专门做多项式最小二乘的接口返回从高次到低次的系数。coeffs np.polyfit(x, y, deg3) y_pred np.polyval(coeffs, x_smooth)np.polyfit底层也是最小二乘只是帮你把设计矩阵的构造过程封装好了。明白它的本质之后你就会知道它和上面手写设计矩阵的差别主要是接口层面的而不是数学层面的。3.3 sklearn一行完成线性回归如果你已经安装了 scikit-learn直接用它的线性回归模型会更省事尤其是在有多维特征时。import numpy as np from sklearn.linear_model import LinearRegression x np.array([1.0, 2.0, 3.0, 4.0, 5.0]).reshape(-1, 1) y np.array([2.1, 3.9, 6.2, 8.1, 9.8]) model LinearRegression() model.fit(x, y) print(斜率, model.coef_) print(截距, model.intercept_)注意这里x需要是二维数组因为 sklearn 的设计是所有特征都放在列上哪怕只有一个特征。fit之后coef_保存着每个特征的系数intercept_是截距。用这个模型还能继续调用predict对新数据进行预测。sklearn 的LinearRegression默认也会对设计矩阵做中心化处理然后使用类似分解的方式来求解所以大多数情况下不用手工处理均值。如果你拿到的特征量纲差异很大建议先做标准化因为这会直接影响后续的正则化或者解释性。4. 实战拟合一组带噪声的传感器数据并做诊断4.1 生成模拟数据和可视化散点只看玩具数据还是会觉得不过瘾所以这里模拟一个更接近实际的场景。假设某次传感器校准实验中真实输入输出关系是线性的但在采集过程中混入了随机噪声。我们用代码生成这样一组数据import numpy as np import matplotlib.pyplot as plt rng np.random.default_rng(42) x np.linspace(0, 10, 50) y_true 1.6 * x 0.5 y y_true rng.normal(0, 1.2, sizex.size) plt.scatter(x, y, alpha0.6, label采样数据) plt.plot(x, y_true, colorgreen, linestyle--, label真实关系) plt.xlabel(x) plt.ylabel(y) plt.legend() plt.show()选择随机噪声标准差为 1.2这会让数据点围绕真实直线上下波动但不会完全掩盖线性趋势。使用固定随机种子 42是为了让结果可复现。如果你换一个随机种子拟合结果会有细微差异但整体趋势应该差不多。实际的传感器数据往往比这更脏比如存在系统误差、重复采样的偏差、甚至部分工作状态异常导致的数据点突变。这里先从理想化的“高斯噪声”开始是为了后续能看到最小二乘的正常表现。4.2 求解回归系数并计算评估指标拟合本身很简单但评估不能少。我会同时计算三个指标均方误差、均方根误差和决定系数 (R^2)。X np.vstack([np.ones_like(x), x]).T theta, residuals, rank, s np.linalg.lstsq(X, y, rcondNone) beta0, beta1 theta y_pred X theta mse np.mean((y - y_pred) ** 2) rmse np.sqrt(mse) ss_res np.sum((y - y_pred) ** 2) ss_tot np.sum((y - np.mean(y)) ** 2) r2 1 - ss_res / ss_tot print(f拟合斜率为 {beta1:.3f}) print(f拟合截距为 {beta0:.3f}) print(fMSE 为 {mse:.4f}) print(fRMSE 为 {rmse:.4f}) print(fR^2 为 {r2:.4f})我在这组数据上跑出来的结果大概是斜率 (1.56)截距 (0.68)。真实值是斜率 (1.6)、截距 (0.5)可以看到即便加了噪声参数估计依然比较接近真实值。RMSE 大概在 (1.1) 左右接近甚至略小于生成数据时设定的噪声标准差 (1.2)说明拟合效果是合理的。至于 (R^2)它表示模型能解释数据方差的比例。(R^2) 越接近 1 表示拟合越好但要注意它并不是万能的。只要不断往模型里加特征(R^2) 一定会上升所以在后面多项式拟合里必须小心过拟合。4.3 通过残差图判断模型是否合理拟合结果出来后不要急着下结论还要看一眼残差图。残差是指 (y - \hat{y})也就是实际值和预测值的差。residuals y - y_pred plt.scatter(x, residuals, alpha0.7) plt.axhline(y0, colorred, linestyle--) plt.xlabel(x) plt.ylabel(residual) plt.show()如果模型是合理的残差应该随机地分布在零附近没有明显趋势。比如在某个区间内残差全是正的或者图像呈现明显的“喇叭口”形状那说明你的模型结构可能有问题比如漏掉了非线性项或者误差方差不恒定。在这个模拟数据里残差应该比较均匀地分布在零轴上下没有系统性偏离。这就验证了“线性模型 高斯噪声”的假设是成立的。其实这一步是很多刚接触回归的人容易忽略的。只看 (R^2) 不够残差图才是判断模型是否真正适配的重要工具。5. 避坑指南最小二乘使用中的高频问题5.1 系数矩阵病态导致结果不稳定实际工作中我遇到过一种很典型的情况特征之间相关性极高。比如同时把“温度”和“温度加上某个固定偏移”放进模型设计矩阵的列之间几乎线性相关此时 (X^T X) 的条件数会变得很大。特征相关时的典型表现是参数估计值很大而且正负交替比如某个特征系数是正的 1000另一个是负的 998两者之间互相抵消。模型在训练数据上可能照样拟合得不错但换一批数据结果就完全不一样。这种情况我也见过很多次所以在做线性回归之前建议先看特征的相关矩阵或者计算设计矩阵的条件数cond np.linalg.cond(X) print(cond)条件数非常大时就不要再硬用普通最小二乘了可以用岭回归也就是在 (X^T X) 的对角线上加一个小常数[ \beta (X^T X \lambda I)^{-1} X^T y ]这个 (\lambda) 会惩罚过大的系数让模型参数更稳定。实现起来也很简单lamb 0.1 p X.shape[1] theta_ridge np.linalg.inv(X.T X lamb * np.eye(p)) X.T y需要注意的是截距项通常不应该加惩罚否则结果会依赖于原点的位置。更严谨的做法是先对特征做中心化拟合后再还原截距。当然直接调用 sklearn 的Ridge会更省事。5.2 数据里混入了离群点最小二乘法对离群点非常敏感因为平方误差会让远离模型的数据点获得极大权重。想象三个点在同一条直线上第四个点在很远的位置拟合线会被那个远处的点强力拽过去导致模型对大多数点都不再适用。处理离群点的第一个办法是画散点图先看一眼但高维数据看不了所以更系统的方法是检查残差。如果某个点的残差是其他点的好几倍就要引起警惕。第二个办法是换用带鲁棒性的估计方法比如 RANSAC、Huber 回归、Theil-Sen 回归。sklearn 里都有现成实现from sklearn.linear_model import RANSACRegressor model RANSACRegressor() model.fit(x.reshape(-1, 1), y)RANSAC 的思路是反复随机抽取一部分数据拟合模型统计哪些点是“内点”最后用内点子集重新拟合。这种方法在有离群点的时候远比普通最小二乘可靠。代价是没有闭式解结果会带随机性但工程上这种取舍通常很划算。5.3 多项式阶数过高引发过拟合多项式拟合最容易犯的错就是把阶数拉得很高。比如用 15 次多项式去拟合 50 个带噪声的点模型会把每个噪声细节都学进去画出来的曲线剧烈震荡。用训练集评估时误差很小但一旦预测新数据效果往往很糟糕。我记得有次为了让训练集上的 (R^2) 更好看把多项式阶数从 3 调到了 12结果曲线在两个端点附近冲到了完全离谱的高度。那之后我给自己定了一条规矩永远不要只看训练集上的 (R^2)必须做交叉验证或者至少保留一部分验证数据。如果你必须使用高阶多项式可以考虑加正则化把损失函数改成[ L |y - X\beta|^2 \lambda |\beta|_2^2 ]这就是前面说的岭回归。它能在不显著降低拟合精度的前提下限制系数不要太大从而抑制过拟合。另一个做法是减少特征用更合理的基函数或者干脆换一种非参数模型。5.4 线性最小二乘真的不适用时怎么办不是所有问题都能套线性最小二乘。如果模型本身对参数非线性比如 (y a \exp(-b x) c)那正规方程就不能直接解。此时有两个选择一是对模型做变形比如取对数后变成线性形式但这会改变噪声的假设二是直接用非线性最小二乘优化算法scipy 提供了现成接口from scipy.optimize import curve_fit def model_func(x, a, b, c): return a * np.exp(-b * x) c popt, pcov curve_fit(model_func, x, y, p0[1, 0.5, 0.1])curve_fit内部会用莱文贝格-马夸特算法迭代求解好处是灵活坏处是需要给初始值。初始值给得太差可能收敛到局部最优所以尽量结合数据范围给一个合理猜测。还有一种情况是你要预测的目标不是连续数值而是类别标签。那就不再是回归问题应该用逻辑回归、支持向量机、决策树这些分类模型。最小二乘是回归家族的基石但它不是万能工具箱。判断问题类型是建模之前就要想清楚的事。最后分享一个我自己的习惯。拿到数据后我一般会先画散点图和残差图心里对数据形态有个大概判断再决定要不要用最小二乘、用什么样的特征组合。很多人一上来就敲代码调库最后参数出来了却不知道模型是不是真的合理。其实最小二乘最花时间的部分往往不是方程求解而是对数据和模型假设的检查。这个习惯帮我避开了很多看起来拟合得很好、实际上完全不能用的回归模型。
