我第一次真正用上最小二乘法，是在处理一组传感器校准实验数据的时候。当时手里是一堆带噪声的采样点，教科书里的公式能帮我“找一条最合理的直线”，但工程上真正让人头疼的，其实是后面这些事：数据里混了离群点怎么办，多项式阶数选多少才不会过拟合，拟合出的结果到底可不可信。这篇文章想做的，就是把最小二乘法从数学原理讲到代码实现，再用一个带噪声的拟合案例从头到尾串一遍。适合刚接触数值计算、需要做回归拟合，或者想看明白各种回归库背后在做什么的开发者。

最小二乘法并不是什么高深技巧，它更像是一套统一的“找最好参数”的思路。只要你手里有一组观测数据，想用一个模型去描述它背后的规律，最小二乘就是最常用的切入方式。下面我会从直觉出发，把公式一步步拆开，再给出可直接运行的代码，最后聊几个我自己踩过的坑。

1. 最小二乘法的核心思想：从误差平方说起

1.1 什么样的拟合才算“最好”

先回到最简单的场景。给你一堆平面上的点，大概分布在一条直线附近，你要画一条直线穿过它们。问题来了：什么叫“最合理”？每个人都能画一条线，但哪条线才是公认最好的？

最自然的想法是，让这条直线尽可能靠近每一个点。也就是说，对每一个数据点，算出我预测的位置和实际位置之间的距离，把这些距离加起来，和最小的那条线就是最好的。这里的“距离”在统计里叫“残差”，也就是观测值和模型预测值之差。最小二乘法的名字就是这么来的：让残差的平方和最小。

[ \min_{\beta_0,\beta_1} \sum_{i=1}^{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_{i=1}^{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_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i) = 0 ]

对 (\beta_1) 求偏导：

[ \frac{\partial L}{\partial \beta_1} = -2 \sum_{i=1}^{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 (p+1)) 矩阵，(\beta) 是 (p+1) 维列向量，(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.lstsq：

theta, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None) 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, rcond=None) 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, s=12, alpha=0.7, label="带噪声数据") plt.plot(x_smooth, y_pred, color="red", label="四次特征组合拟合") plt.legend() plt.show()

这里用的是三次多项式，所以 (X) 有四列：常数项、一次项、二次项、三次项。最小二乘会自动估计出每个特征的权重。如果你的数据有周期性，三次多项式通常拟合得不太好，但这并不妨碍你理解流程：只要构造了合适的特征矩阵，最小二乘的求解逻辑没有任何变化。

多项式拟合还有一个更偷懒的函数：np.polyfit。它是专门做多项式最小二乘的接口，返回从高次到低次的系数。

coeffs = np.polyfit(x, y, deg=3) 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, size=x.size) plt.scatter(x, y, alpha=0.6, label="采样数据") plt.plot(x, y_true, color="green", 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, rcond=None) 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(f"MSE 为 {mse:.4f}") print(f"RMSE 为 {rmse:.4f}") print(f"R^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, alpha=0.7) plt.axhline(y=0, color="red", 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内部会用莱文贝格-马夸特算法迭代求解，好处是灵活，坏处是需要给初始值。初始值给得太差可能收敛到局部最优，所以尽量结合数据范围给一个合理猜测。

还有一种情况是你要预测的目标不是连续数值，而是类别标签。那就不再是回归问题，应该用逻辑回归、支持向量机、决策树这些分类模型。最小二乘是回归家族的基石，但它不是万能工具箱。判断问题类型是建模之前就要想清楚的事。

最后分享一个我自己的习惯。拿到数据后我一般会先画散点图和残差图，心里对数据形态有个大概判断，再决定要不要用最小二乘、用什么样的特征组合。很多人一上来就敲代码调库，最后参数出来了却不知道模型是不是真的合理。其实最小二乘最花时间的部分往往不是方程求解，而是对数据和模型假设的检查。这个习惯帮我避开了很多看起来拟合得很好、实际上完全不能用的回归模型。