1. 从一次拟合失败的实战经历说起为什么我们需要最小二乘法我先讲个自己的经历。几年前做传感器标定采集到一组电压与温度的对应数据我拿高次多项式硬怼结果拟合曲线像条疯狗一样在采样点之间乱窜完全没法用。后来改成最小二乘拟合一条直线下去误差分布漂亮得像个正态分布。那是我第一次真切理解拟合不是穿过每个点而是在整体上最接近所有点。这个思想就是最小二乘法的灵魂。最小二乘法Least Squares Method是数学、统计、工程领域最基础也最被低估的工具之一。它要解决的核心问题可以一句话说清**给定一批散点数据如何找到一条曲线或一个函数使得它跟所有点的整体距离最小。**这里的距离不是简单的垂直距离而是误差的平方和。为什么用平方而不是绝对值为什么最小二乘解恰好是那个正规方程很多教程直接甩出公式就完事了但对真正想搞懂的人来说缺了中间的推导过程就永远只是会用不是懂。这篇博文的目标只有一个从零开始把最小二乘法的原理、推导、几何意义、概率解释、实战操作全部串起来让你看完之后不仅能独立推导出公式还能根据实际问题判断什么时候该用、怎么用、用了之后结果怎么解读。无论你是刚学线性代数的学生、搞数据分析的从业者、做机器学习建模的工程师还是需要处理实验数据的科研人员只要你在和数据打交道这篇内容就值得你花二十分钟彻底读一遍。2. 核心思想拆解误差最小化与平方损失函数2.1 拟合的起点一组方程一个未知的规律先从最经典的场景说起——一元线性回归。假设我们观测到一组数据[ (x_1, y_1), (x_2, y_2), \dots, (x_n, y_n) ]并且我们相信 (y) 和 (x) 之间大致是线性关系[ y ax b ]但因为观测误差、噪声、未建模因素的存在这 (n) 个点并不会老老实实全落在同一条直线上。所以我们面临的问题是怎么确定参数 (a)斜率和 (b)截距让这条直线最合理地代表这组数据最朴素的想法是让每个数据点到直线的距离尽量小。于是定义第 (i) 个点的误差残差为[ e_i y_i - (ax_i b) ]这个 (e_i) 可能是正的也可能是负的。如果把所有误差直接相加正负会互相抵消哪怕拟合效果一塌糊涂误差之和可能依然很小这显然不合理。所以需要一种能累积误差大小的方式。处理方式通常有三种候选误差绝对值之和(\sum |e_i|)误差平方之和(\sum e_i^2)误差四次方之和(\sum e_i^4)最小二乘法选择的是第二种——平方和即最小化[ S(a, b) \sum_{i1}^{n} (y_i - ax_i - b)^2 ]2.2 为什么偏偏是平方两个关键原因第一个原因是数学上的平方函数处处可导而且是凸函数。这意味着我们可以直接用微积分求极值——对 (a)、(b) 求偏导令其为零就能得到唯一的最优解。绝对值函数在零点不可导四次方函数也会让求解过程变得复杂得多。第二个原因是统计上的在误差服从正态分布的假设下最小化平方和等价于最大化似然函数这个我在第4节会详细推导。也就是说最小二乘法在统计意义上具有坚实的理论基础。平方不是拍脑袋定的而是从概率原理中自然涌现出来的。打个生活中的比方你要估计一群人平均身高如果目标是让猜测值和所有人的差距总和最小用平方误差得到的就是平均值用绝对值误差得到的是中位数。二者各有用途但平均值在数学性质上更好处理也是大多数人直觉中的代表值。2.3 从一元到多元问题升级为矩阵形式一元线性回归解决了一条直线的问题。但当变量增多时比如我们要拟合[ y \beta_0 \beta_1 x_1 \beta_2 x_2 \cdots \beta_p x_p ]就不能再用 (a, b) 两个字母了需要统一的矩阵语言。把所有观测数据写成矩阵形式[ \mathbf{y} \begin{bmatrix} y_1 \ y_2 \ \vdots \ y_n \end{bmatrix}, \quad \mathbf{X} \begin{bmatrix} 1 x_{11} x_{12} \cdots x_{1p} \ 1 x_{21} x_{22} \cdots x_{2p} \ \vdots \vdots \vdots \ddots \vdots \ 1 x_{n1} x_{n2} \cdots x_{np} \end{bmatrix}, \quad \boldsymbol{\beta} \begin{bmatrix} \beta_0 \ \beta_1 \ \vdots \ \beta_p \end{bmatrix} ]这里 (\mathbf{X}) 的第一列全是 1对应截距项 (\beta_0)。于是模型可以简洁地写成[ \mathbf{y} \mathbf{X}\boldsymbol{\beta} \boldsymbol{\varepsilon} ]其中 (\boldsymbol{\varepsilon}) 是误差向量。我们的目标变成找到一个 (\boldsymbol{\beta})使得残差平方和[ S(\boldsymbol{\beta}) |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2 ]最小。这个形式是所有最小二乘问题线性的统一表达后面的一切推导都从这里出发。3. 代数视角推导从偏导为零到正规方程3.1 标量形式的求导解法先做一遍最简单的一元线性回归推导把直觉建立起来。目标函数为[ S(a, b) \sum_{i1}^{n} (y_i - ax_i - b)^2 ]对 (b) 求偏导[ \frac{\partial S}{\partial b} -2 \sum_{i1}^{n} (y_i - ax_i - b) 0 ]化简得[ \sum_{i1}^{n} y_i - a\sum_{i1}^{n} x_i - nb 0 ]即[ b \bar{y} - a\bar{x} \tag{1} ]其中 (\bar{x})、(\bar{y}) 分别是 (x)、(y) 的均值。对 (a) 求偏导[ \frac{\partial S}{\partial a} -2 \sum_{i1}^{n} x_i (y_i - ax_i - b) 0 ]即[ \sum_{i1}^{n} x_i y_i - a\sum_{i1}^{n} x_i^2 - b\sum_{i1}^{n} x_i 0 \tag{2} ]把 (1) 代入 (2)[ \sum_{i1}^{n} x_i y_i - a\sum_{i1}^{n} x_i^2 - (\bar{y} - a\bar{x})\sum_{i1}^{n} x_i 0 ]利用 (\sum x_i n\bar{x})整理过程如下[ \sum x_i y_i - a\sum x_i^2 - n\bar{x}\bar{y} a n \bar{x}^2 0 ][ \sum x_i y_i - n\bar{x}\bar{y} a\left(\sum x_i^2 - n\bar{x}^2\right) ]所以[ a \frac{\sum x_i y_i - n\bar{x}\bar{y}}{\sum x_i^2 - n\bar{x}^2} ]进一步利用恒等式[ \sum (x_i - \bar{x})(y_i - \bar{y}) \sum x_i y_i - n\bar{x}\bar{y} ][ \sum (x_i - \bar{x})^2 \sum x_i^2 - n\bar{x}^2 ]于是得到教科书中常见的形式[ \boxed{a \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2}} ][ \boxed{b \bar{y} - a\bar{x}} ]这就是一元线性回归的最小二乘解。分子的含义是(x) 和 (y) 的协方差分母是(x) 的方差。3.2 矩阵形式的完整推导过程现在挑战更通用的多元情形。我们把残差平方和展开[ S(\boldsymbol{\beta}) (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) ]展开后得到三项[ S(\boldsymbol{\beta}) \mathbf{y}^T\mathbf{y} - \boldsymbol{\beta}^T\mathbf{X}^T\mathbf{y} - \mathbf{y}^T\mathbf{X}\boldsymbol{\beta} \boldsymbol{\beta}^T\mathbf{X}^T\mathbf{X}\boldsymbol{\beta} ]由于 (\boldsymbol{\beta}^T\mathbf{X}^T\mathbf{y}) 和 (\mathbf{y}^T\mathbf{X}\boldsymbol{\beta}) 都是标量1×1 矩阵它们互为转置且相等所以合并为一项[ S(\boldsymbol{\beta}) \mathbf{y}^T\mathbf{y} - 2\boldsymbol{\beta}^T\mathbf{X}^T\mathbf{y} \boldsymbol{\beta}^T\mathbf{X}^T\mathbf{X}\boldsymbol{\beta} ]接下来对 (\boldsymbol{\beta}) 求梯度。这里用到矩阵微积分中的几个常用公式(\frac{\partial (\boldsymbol{\beta}^T\mathbf{A})}{\partial \boldsymbol{\beta}} \mathbf{A})(\frac{\partial (\boldsymbol{\beta}^T\mathbf{B}\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} 2\mathbf{B}\boldsymbol{\beta})当 (\mathbf{B}) 对称注意 (\mathbf{X}^T\mathbf{X}) 是实对称矩阵所以[ \frac{\partial S}{\partial \boldsymbol{\beta}} -2\mathbf{X}^T\mathbf{y} 2\mathbf{X}^T\mathbf{X}\boldsymbol{\beta} ]令梯度等于零得到[ \mathbf{X}^T\mathbf{X}\boldsymbol{\beta} \mathbf{X}^T\mathbf{y} ]这就是著名的正规方程Normal Equation。若 (\mathbf{X}^T\mathbf{X}) 可逆则最小二乘估计量为[ \boxed{\hat{\boldsymbol{\beta}} (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}} ]整个推导过程干净利落只用到了矩阵转置、乘法和一个求导规则。3.3 驻点为什么一定是最小值凸性论证有人可能会问令梯度为零得到的点凭什么是最小值而不是最大值或鞍点关键在于 (S(\boldsymbol{\beta})) 的二阶导Hessian 矩阵[ \mathbf{H} \frac{\partial^2 S}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta}^T} 2\mathbf{X}^T\mathbf{X} ]对于任意非零向量 (\mathbf{v})有[ \mathbf{v}^T \mathbf{X}^T\mathbf{X} \mathbf{v} |\mathbf{X}\mathbf{v}|^2 \geq 0 ]这说明 (\mathbf{X}^T\mathbf{X}) 是半正定矩阵。如果 (\mathbf{X}) 列满秩则它是正定矩阵此时 (S(\boldsymbol{\beta})) 是严格凸函数有唯一的全局最小值点。所以正规方程的驻点一定是最小值。这个性质在实际中意义很大最小二乘问题永远不会遇到局部最优解陷阱。只要 (\mathbf{X}^T\mathbf{X}) 可逆解就是唯一的、全局最优的。这也是为什么最小二乘法在工程中如此受欢迎——稳定、可靠、不依赖初始化。4. 几何视角解读最小二乘解为什么是投影4.1 把拟合理解为向列空间投影代数推导虽然严谨但很多人会觉得抽象。我换个角度——从线性代数的几何意义来看。回忆 (\mathbf{X}) 是 (n \times (p1)) 矩阵它的每一列可以看作高维空间 (\mathbb{R}^n) 中的一个向量。由这些列向量张成的子空间记作 (\text{Col}(\mathbf{X}))它包含了所有 (\mathbf{X}\boldsymbol{\beta}) 可能的取值——也就是所有模型能够表达出来的预测值。真正观测到的 (\mathbf{y}) 一般不在这个子空间里因为存在噪声。我们做的事情是找一个子空间中的点 (\hat{\mathbf{y}} \mathbf{X}\hat{\boldsymbol{\beta}})让它与 (\mathbf{y}) 的距离最小。这恰恰就是投影的定义。高等代数中有个基本结论(\mathbf{y}) 到子空间 (\text{Col}(\mathbf{X})) 的投影向量 (\hat{\mathbf{y}})满足残差向量 (\mathbf{y} - \hat{\mathbf{y}}) 与子空间正交。正交意味着对 (\mathbf{X}) 的每一列都垂直[ \mathbf{X}^T (\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}) \mathbf{0} ]即[ \mathbf{X}^T \mathbf{y} \mathbf{X}^T\mathbf{X}\hat{\boldsymbol{\beta}} ]这不是别的正是正规方程投影矩阵 (\mathbf{P} \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T) 把 (\mathbf{y}) 变成 (\hat{\mathbf{y}})[ \hat{\mathbf{y}} \mathbf{P}\mathbf{y} ]4.2 投影视角下的解的唯一性条件如果 (\mathbf{X}) 列满秩那么 (\mathbf{X}^T\mathbf{X}) 可逆投影矩阵唯一确定最小二乘解也唯一。如果 (\mathbf{X}) 列不满秩比如存在多重共线性那么 (\mathbf{X}^T\mathbf{X}) 奇异正规方程有无穷多解但投影向量 (\hat{\mathbf{y}}) 依然是唯一的。通俗地说解不唯一但预测值唯一。这一点在实操中非常重要。很多人发现 (\mathbf{X}^T\mathbf{X}) 不可逆时一脸懵其实这不代表数据没有拟合的意义只是说明特征之间存在线性相关。此时可以用岭回归加正则项或伪逆来求解。4.3 一个简单的二维几何直觉想象三维空间中有一张纸这相当于 (\mathbf{X}) 列张成的二维子空间纸外有一个点相当于 (\mathbf{y})。你从那个点垂直向纸面打一束光影子落在纸上的位置就是投影 (\hat{\mathbf{y}})。从这个点到影子点的连线就是残差向量它跟纸面上的任何方向都垂直——这就是正交性。最小二乘法做的本质上就是打光打出来找影子。它不要求数据点在直线上点在纸面上而是承认噪声的存在接受一种最优的近似。5. 概率学基础正态分布假设与最大似然推导5.1 凭什么最小化平方和从概率角度重新发现很多人不知道的是最小二乘法的平方损失函数不是人为规定的而是在误差服从正态分布的假设下从最大似然估计中自然推导出来的。假设数据生成过程为[ y_i \mathbf{x}_i^T\boldsymbol{\beta} \varepsilon_i ]其中误差项 (\varepsilon_i \sim \mathcal{N}(0, \sigma^2))且彼此独立同分布这是高斯-马尔可夫定理的基础假设。那么给定 (\mathbf{x}_i) 时(y_i) 的条件分布为[ y_i | \mathbf{x}_i \sim \mathcal{N}(\mathbf{x}_i^T\boldsymbol{\beta}, \sigma^2) ]其概率密度函数为[ p(y_i | \mathbf{x}_i, \boldsymbol{\beta}) \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - \mathbf{x}_i^T\boldsymbol{\beta})^2}{2\sigma^2} \right) ]5.2 最大似然估计的完整推导给定全部 (n) 个独立样本似然函数为[ L(\boldsymbol{\beta}) \prod_{i1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - \mathbf{x}_i^T\boldsymbol{\beta})^2}{2\sigma^2} \right) ]取对数对数函数单调递增不影响最大化问题[ \ell(\boldsymbol{\beta}) -\frac{n}{2}\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{i1}^{n} (y_i - \mathbf{x}_i^T\boldsymbol{\beta})^2 ]要最大化 (\ell(\boldsymbol{\beta}))只需最小化最后一项[ \sum_{i1}^{n} (y_i - \mathbf{x}_i^T\boldsymbol{\beta})^2 ]这不就是最小二乘的损失函数吗换句话说正态噪声假设下最小二乘估计 最大似然估计。平方和不仅方便而且在统计意义上是最优的。5.3 最小一乘呢拉普拉斯噪声下的世界如果误差服从的不是正态分布而是拉普拉斯分布厚尾分布那么最大似然估计对应的损失函数就变成绝对值损失[ \sum_{i1}^{n} |y_i - \mathbf{x}_i^T\boldsymbol{\beta}| ]这就是最小一乘法最小绝对偏差。它的优点是抗离群值干扰能力强缺点是目标函数在零点不可导求解需要使用线性规划等方法计算难度更大。两相对比就能理解最小二乘法的平方选择既是数学上便于求解的结果也是统计上对正态噪声的最优回应。5.4 高斯-马尔可夫定理为什么最小二乘是最好的高斯-马尔可夫定理Gauss-Markov Theorem是理解最小二乘法地位的关键结论。它说在误差满足以下条件时——误差期望为零(E[\varepsilon_i] 0)误差方差恒定(\text{Var}(\varepsilon_i) \sigma^2)同方差性误差互不相关(\text{Cov}(\varepsilon_i, \varepsilon_j) 0)独立不相关——最小二乘估计 (\hat{\boldsymbol{\beta}}) 是所有线性无偏估计中方差最小的。也就是说在所有形如 (\tilde{\boldsymbol{\beta}} \mathbf{C}\mathbf{y}) 的无偏估计中最小二乘估计具有最小的方差。不用 (\sqrt{|e|})、不用四次方就是平方误差下的解最优。这个定理把最小二乘法从一个合理的方案提升到了一个在特定假设下不可超越的方案。当然如果数据违反同方差性出现异方差或者存在严重的离群值那么加权最小二乘或鲁棒回归会更合适——后面会说到。6. 手算实例从数据到参数完整走一遍6.1 问题设定假设我们要研究广告投入与销售额之间的关系。收集到5条数据广告投入 (x)万元销售额 (y)万元132537410512目标是拟合一元线性模型(y ax b)。6.2 分步计算先算各变量[ \sum x_i 15, \quad \sum y_i 37, \quad \bar{x} 3, \quad \bar{y} 7.4 ][ \sum x_i^2 1491625 55 ][ \sum x_i y_i 1\times3 2\times5 3\times7 4\times10 5\times12 310214060 134 ]接着[ \sum (x_i - \bar{x})(y_i - \bar{y}) \sum x_i y_i - n\bar{x}\bar{y} 134 - 5 \times 3 \times 7.4 134 - 111 23 ][ \sum (x_i - \bar{x})^2 \sum x_i^2 - n\bar{x}^2 55 - 5 \times 9 55 - 45 10 ]于是[ a \frac{23}{10} 2.3 ][ b 7.4 - 2.3 \times 3 7.4 - 6.9 0.5 ]最终拟合直线为[ \boxed{y 2.3x 0.5} ]直观检查当 (x1) 时预测值 (2.8)实际值 (3)当 (x5) 时预测值 (12)实际值 (12)。这组数据拟合效果相当好。残差平方和((3-2.8)^2 (5-5.1)^2 (7-7.4)^2 (10-9.7)^2 (12-12)^2 0.040.010.160.090 0.30)。若用绝对值方法最优解可能就不是这个值了这正好呼应前面概率部分的讨论。6.3 系数解读斜率与截距的业务含义斜率 (a 2.3)广告投入每增加1万元销售额平均增加2.3万元。截距 (b 0.5)当广告投入为0时基础销售额约为0.5万元。注意这里是外推区间0不在数据范围(x \in [1,5])内所以这个解读要谨慎——它只是数学上的截距业务上未必有实际含义。这就是最小二乘结果解读中常见的一个坑截距项的解读要结合数据范围不能机械套用。7. 从正规方程到数值计算矩阵求逆的局限与解法升级7.1 正规方程在实际计算中的三个问题正规方程 (\boldsymbol{\beta} (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}) 形式非常漂亮但直接按这个思路编程实现会遇到三个现实问题第一个是计算复杂度。求一个 (p \times p) 矩阵的逆计算量大约是 (O(p^3))。当特征数 (p) 很大时比如几万维这个运算量会爆炸。机器学习中经常出现特征数远大于样本量的情况(p \gg n)(\mathbf{X}^T\mathbf{X}) 的维度比数据还大正规方程直接被淘汰。第二个是数值稳定性。(\mathbf{X}^T\mathbf{X}) 的条件数是原矩阵 (\mathbf{X}) 条件数的平方。如果彼此相关的特征让 (\mathbf{X}) 条件数偏大(\mathbf{X}^T\mathbf{X}) 就可能出现严重的数值误差求逆结果不可靠。第三个是奇异不可逆。当特征之间存在完全共线性或者样本量小于特征数时(\mathbf{X}^T\mathbf{X}) 不可逆正规方程直接无解。7.2 稳定解法QR分解、SVD与伪逆针对以上问题工程上通常不直接计算矩阵逆而是用矩阵分解QR分解把 (\mathbf{X}) 分解为 (\mathbf{Q}\mathbf{R})(\mathbf{Q}) 是正交矩阵(\mathbf{R}) 是上三角矩阵然后解三角方程组 (\mathbf{R}\boldsymbol{\beta} \mathbf{Q}^T\mathbf{y})。计算量 (O(np^2))数值稳定性比正规方程好得多。SVD分解把 (\mathbf{X}) 分解为 (\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^T)最小二乘解为 (\boldsymbol{\beta} \mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{U}^T\mathbf{y})。SVD能处理奇异矩阵用伪逆而且对病态矩阵最鲁棒是今天所有主流数值计算库的标准底层方案。伪逆(\mathbf{X}^ \mathbf{V}\boldsymbol{\Sigma}^{-1}\mathbf{U}^T)也叫 Moore-Penrose 逆。即使 (\mathbf{X}^T\mathbf{X}) 不可逆伪逆依然能给出最小范数解。Python 的numpy.linalg.lstsq、Matlab 反斜杠运算符、R 语言的lm()底层用的都是这类分解方法而不是显式求逆。理解这一点很重要——写代码时永远不要手动计算 ((\mathbf{X}^T\mathbf{X})^{-1})。7.3 正则化当最小二乘无能为力时当特征数 (p n)高维数据或存在严重多重共线性时最小二乘解不再可靠。标准做法是加正则项岭回归L2最小化 (|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2 \lambda|\boldsymbol{\beta}|^2)解为 ((\mathbf{X}^T\mathbf{X} \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y})。(\lambda\mathbf{I}) 相当于给对角线加了个垫片保证矩阵可逆。LassoL1最小化 (|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2 \lambda|\boldsymbol{\beta}|_1)会得到稀疏解常用于特征筛选。正则化的实质是放弃无偏这一理想性质换取更低的方差和更好的泛化能力。这也印证了第5节内容最小二乘的最优性是有前提的数据一旦不合适算法也要跟着变。8. 模型评估拟合优度与假设检验8.1 R²决定系数的计算与解读拟合完模型之后紧接着的问题是这模型到底好不好定义几个量[ \text{SST总平方和} \sum (y_i - \bar{y})^2 ][ \text{SSR回归平方和} \sum (\hat{y}_i - \bar{y})^2 ][ \text{SSE误差平方和} \sum (y_i - \hat{y}_i)^2 ]三者满足关系[ \text{SST} \text{SSR} \text{SSE} ]决定系数定义为[ R^2 \frac{\text{SSR}}{\text{SST}} 1 - \frac{\text{SSE}}{\text{SST}} ](R^2) 表示模型解释了多少比例的数据波动。比如 (R^2 0.85)意思是 85% 的变异性被模型解释了剩余 15% 属于噪声或未建模因素。用前面例子计算(\hat{y}) 值分别是 (2.8, 5.1, 7.4, 9.7, 12)(\bar{y} 7.4)。SSE 0.30。SST ((3-7.4)^2(5-7.4)^2(7-7.4)^2(10-7.4)^2(12-7.4)^2 19.365.760.166.7621.16 53.2)。于是 (R^2 1 - 0.30/53.2 \approx 0.9944)。拟合效果极好。8.2 调整R²惩罚过多特征(R^2) 有一个特性只要往模型里加变量哪怕这个变量完全是噪音(R^2) 也不会下降通常还会略微上升。为了解决这个问题引入调整决定系数[ \bar{R}^2 1 - \frac{\text{SSE} / (n-p-1)}{\text{SST} / (n-1)} ]其中 (p) 是特征数量。调整 (R^2) 会惩罚多余特征适合比较不同变量数量的模型。8.3 残差诊断最小二乘有效性的最后检查拟合完后强烈建议把所有残差 (e_i y_i - \hat{y}_i) 画出来看看。正常的残差应该在零线附近随机波动没有明显的模式。需要注意三种异常残差随拟合值增大而增大喇叭状存在异方差误差方差不是常数违背高斯-马尔可夫假设。考虑对 (y) 做对数变换或者用加权最小二乘。残差出现曲线模式模型形式不对可能存在非线性关系考虑添加二次项、交互项或换非线性模型。残差出现明显的离群点某条数据严重偏离整体规律检查是不是记录错误或者考虑使用鲁棒回归。我有一个习惯每次做最小二乘拟合必画三张图——残差vs拟合值、残差QQ图、残差vs各特征。三张图看完模型的毛病基本都藏不住。这一步虽然不复杂但能省下后面很多返工的时间。9. 非线性与多变量扩展从最小二乘到更广阔的拟合世界9.1 多项式拟合变量变换仍是最小二乘最小二乘法不限于线性关系。如果数据的趋势是弯曲的可以用多项式[ y \beta_0 \beta_1 x \beta_2 x^2 \cdots \beta_k x^k ]虽然 (y) 和 (x) 的关系非线性但只要把 (x, x^2, \dots, x^k) 看作独立的特征问题就变回线性最小二乘。把原始的 (\mathbf{X}) 矩阵每一列改造一下就行[ \mathbf{X} \begin{bmatrix} 1 x_1 x_1^2 \cdots x_1^k \ 1 x_2 x_2^2 \cdots x_2^k \ \vdots \vdots \vdots \ddots \vdots \ 1 x_n x_n^2 \cdots x_n^k \end{bmatrix} ]接下来所有推导一律照旧。这就是线性在数学中的真实含义对参数线性而非对变量线性。这里有个坑要提醒多项式阶数 (k) 不能太高。高阶多项式会给出很好的训练集拟合(R^2) 接近1但在新数据上会剧烈振荡这就是过拟合。回到我开篇说的标定数据用 7 次多项式硬拟合后插值区间内的曲线像过山车正是这个原因。经验上工程问题优先尝试 (k1) 或 (k2)除非有明确的理论依据支持更高的阶数。9.2 加权最小二乘给不同数据不同的信任度实际数据中不同观测的可靠性往往不一样。比如物理实验中某些测量设备的精度比其他设备高金融数据中近期数据比远期数据更有参考价值。加权最小二乘引入权重 (w_i)[ S(\boldsymbol{\beta}) \sum_{i1}^{n} w_i (y_i - \mathbf{x}_i^T\boldsymbol{\beta})^2 ]矩阵形式写成[ S (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{W} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) ]其中 (\mathbf{W}) 是对角权重矩阵。求梯度并令其为零得到解[ \hat{\boldsymbol{\beta}} (\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}\mathbf{y} ]使用加权最小二乘的典型场景是误差方差与某个变量成正比异方差。此时设置 (w_i 1/\sigma_i^2)方差越大的观测权重越小能得到有效的估计。这本质上也是最大似然思想的应用——每个数据点的可信度不同贡献的权重就不同。9.3 迭代重加权鲁棒回归的基本思路当数据中存在不少离群点时普通最小二乘会被个别极端值拉偏。迭代重加权最小二乘IRLS的做法是反复执行加权最小二乘→根据残差更新权重→再拟合给残差大的点自动降低权重从而抵消离群值的影响。Huber损失就是一种经典的综合损失函数小残差用平方惩罚高斯假设大残差用线性惩罚拉普拉斯假设兼顾效率与鲁棒性。如果数据常常混入异常值这类方法比纯最小二乘稳妥得多。10. 基于 NumPy 的最小二乘实战演练10.1 用 lstsq 完成线性回归说了这么多理论下面动手实践。使用numpy.linalg.lstsq完成最小二乘求解。import numpy as np import matplotlib.pyplot as plt # 模拟数据真实关系 y 2.3x 0.5加一点高斯噪声 rng np.random.default_rng(42) x np.linspace(0, 10, 50) true_slope 2.3 true_intercept 0.5 y true_slope * x true_intercept rng.normal(0, 1.5, sizex.shape) # 构造设计矩阵 X两列第一列全1对应截距 X np.column_stack([np.ones_like(x), x]) # 最小二乘求解 beta, residuals, rank, s np.linalg.lstsq(X, y, rcondNone) print(f拟合斜率: {beta[1]:.4f}, 拟合截距: {beta[0]:.4f})输出效果拟合斜率接近 2.3截距接近 0.5噪声越大误差越大。lstsq底层使用 LAPACK 的gelsd例程基于SVD数值稳健性远超手动求逆。10.2 正规方程手动求解与对比作为验证可以用正规方程手动算一遍# 正规方程 beta_normal np.linalg.inv(X.T X) X.T y print(f正规方程斜率: {beta_normal[1]:.4f}, 截距: {beta_normal[0]:.4f}) # 与 lstsq 结果对比 print(f参数差异: {np.max(np.abs(beta - beta_normal)):.2e})两种方法的数值差异通常在 (10^{-12}) 量级在数据量小时但一旦特征数增多或矩阵病态差异就会拉开。10.3 多项式拟合与过拟合演示用同一个数据集尝试不同阶数的多项式拟合def poly_fit(x, y, degree): # 构造范德蒙德矩阵 X_poly np.vander(x, degree 1, increasingTrue) coef, _, _, _ np.linalg.lstsq(X_poly, y, rcondNone) return coef # 在训练集上拟合不同阶数 for degree in [1, 3, 9]: coef poly_fit(x, y, degree) y_pred np.polyval(coef[::-1], x) r2 1 - np.sum((y - y_pred) ** 2) / np.sum((y - np.mean(y)) ** 2) print(f阶数 {degree}: 训练集 R² {r2:.4f})试试就知道阶数越高训练集上 (R^2) 越接近 1但一旦在测试集上评估高阶模型的泛化能力反而大幅恶化。这是最小二乘实践中最重要的从理论到应用的提醒最小二乘给出的是数据拟合最优解不是预测最优解。两者之间的差距全靠你对模型复杂度的控制来弥合。10.4 实操中的几个细节第一数据标准化。当不同特征的量纲差异巨大比如一个特征是 0~1另一个是 0~100000模型会对量纲大的特征更敏感数值计算也会不稳定。建议先做标准化减均值除以标准差拟合后再换算回原始尺度。第二截距项别漏。很多人用np.linalg.lstsq时忘记加常数列导致拟合直线被迫过原点结果惨不忍睹。除非有明确的物理依据比如数据本身保证过零点否则永远加一列1。第三检查条件数。用np.linalg.cond(X)查看设计矩阵的条件数。条件数过大比如超过 (10^{12})说明矩阵严重病态结果并不可靠这时候应该做特征筛选或改用正则化方法。cond_number np.linalg.cond(X) print(f设计矩阵条件数: {cond_number:.2e})11. 通用求解模板遇到拟合问题直接用最后给一个可以直接套用的最小二乘求解模板覆盖从数据清洗到模型评估的全流程import numpy as np def least_squares_fit(X, y, add_interceptTrue, degree1): 通用最小二乘拟合函数 X: 原始特征矩阵 (n_samples, n_features) y: 目标向量 (n_samples,) add_intercept: 是否添加截距项 degree: 多项式扩展的阶数 X_use X.copy() # 添加截距列 if add_intercept: X_use np.column_stack([np.ones(X_use.shape[0]), X_use]) # 多项式特征扩展 if degree 1: features [X_use] for d in range(2, degree 1): features.append(X ** d) X_use np.column_stack(features) # 最小二乘求解 beta, residuals, rank, s np.linalg.lstsq(X_use, y, rcondNone) # 预测 y_pred X_use beta # 评估 ss_res np.sum((y - y_pred) ** 2) ss_tot np.sum((y - np.mean(y)) ** 2) r2 1 - ss_res / ss_tot # 调整R² n X_use.shape[0] p X_use.shape[1] - 1 adj_r2 1 - (ss_res / (n - p - 1)) / (ss_tot / (n - 1)) return { coefficients: beta, r_squared: r2, adjusted_r_squared: adj_r2, residuals: y - y_pred, y_pred: y_pred }这套流程适合绝大多数线性最小二乘问题。需要特别说明的是真正做研究或工业落地时我建议直接用scikit-learn的LinearRegression或statsmodels的OLS内部实现经过大量优化且带有完整的统计检验输出。自己写lstsq是为了理解原理生产环境用成熟库是为了少踩坑。12. 最后再交代几句实战经验做最小二乘这些年有几条经验一直刻在脑子里借这个机会分享给读到这里的你。第一拟合前先画散点图。不要急着算系数先肉眼观察数据长什么样是线性是弯曲有没有明显离群点这一步能避免80%的模型形式错误。数据可视化永远是最便宜的诊断工具。第二小心外推。最小二乘模型在训练数据范围内是可靠的一旦跨越出这个范围预测结果就可能完全失真。尤其是多项式拟合高阶项在外推时会迅速发散。我见过太多初学的人拿回归模型预测训练区间之外的数据然后对离谱的结果百思不得其解。第三不要迷信 (R^2)。(R^2 0.99) 不一定代表模型好可能只是数据太简单或者过拟合(R^2 0.3) 也不一定代表模型没用可能问题本身的信噪比就低。模型好不好终究要看残差是否随机、预测是否准确、业务上是否站得住脚。第四理解每一行代码背后的算法。现代数据分析工具越来越傻瓜化一行fit就出结果。但工具越黑盒越需要你理解底层原理。真遇到数据异常、结果诡异的情况能帮你定位问题的不是更高的库版本而是对最小二乘法本质的把握。最小二乘法看起来简单但它背后串起的线性代数、概率统计、数值计算方法是数据科学最重要的基础设施。把这套推导从头到尾自己动手推一遍你的数据处理能力会上升一个台阶。这值得你花掉的每一个小时。