超定方程组与最小二乘：从正规方程到石油产量预测的工程实践

简介：这份资源是广东工业大学数值分析课程第四次上机实验的完整报告，面向正在学习数值分析、需要完成Python上机实验的高校学生。内容围绕超定方程组的最小二乘解、最小二乘拟合以及石油产量预测三个典型问题展开，涵盖算法原理推导、正规方程构建、结果分析与完整源代码，读者可据此理解最小二乘法的数学表述与Python实现路径。资源包共1个docx文件，约112KB，以实验报告形式组织，包含题目、算法介绍、结果及分析和源代码等模块，结构清晰便于对照学习。目前已有213人学习下载，适合需要参考实验思路、核对求解过程或补充数值计算实践的学生使用。

1. 超定方程组：为什么方程比未知数多反而有解

四个方程、两个未知数，按线性代数的常识这组方程应该无解——四条直线在平面上不可能交于一点。但数值分析给出的答案是：它有一个"最优解"。这个反直觉的结论，正是最小二乘法的起点。

现实里的测量数据几乎不可能完美满足所有方程：传感器有噪声、采样有误差、模型本身是近似。方程个数多于未知数个数，恰恰是"信息冗余"的常态，而不是异常。超定方程组的最小二乘解，就是在这种冗余里找一组让所有方程残差平方和最小的参数。搞懂它，等于搞懂了从数据拟合、线性回归到神经网络损失函数的一整条底层逻辑链。这份实验报告把超定方程组、非线性拟合、石油产量预测三个问题串在一起，正好构成一条从"解方程"到"做预测"的完整路径。

2. 正规方程：ATAx=ATb 是怎么把无解变成有解的

2.1 从"无解"到"最近点"的几何直觉

先把问题翻译成几何语言。Ax=b 中，A 是一个 m×n 的矩阵（m>n），它的列向量张成 R^m 里的一个 n 维子空间。b 一般不在这个子空间里，所以方程无解。最小二乘做的事，是在这个子空间里找一个离 b 最近的点 Ax*，让 ‖Ax−b‖² 最小。

最近点的几何含义是：b−Ax* 必须垂直于 A 的列空间。垂直的代数表达就是 Aᵀ(b−Ax*)=0，展开即 AᵀAx*=Aᵀb。这就是正规方程的来源——它不是凭空规定的公式，而是"垂线最短"这一几何事实的直接翻译。

这个推导解释了一个关键点：最小二乘解一定存在，因为 AᵀA 是 n×n 的对称半正定矩阵，只要 A 的列满秩，AᵀA 就可逆，解唯一。实验里两个超定方程组都满足列满秩，所以 np.linalg.solve 能直接给出唯一解。

2.2 源码级拆解：三行 numpy 背后的两次矩阵乘法

实验报告里求超定方程组 a 的代码只有三行核心逻辑，但每一行都对应正规方程的一个组成部分：

import numpy as np # 超定方程组 a：4 个方程，2 个未知数 A = np.array([[2, 4], [3, -5], [1, 2], [2, 1]]) b = np.array([[11], [3], [6], [7]]) # 第一步：构造正规方程的系数矩阵 A^T A（2x2） x = np.dot(A.T, A) # 第二步：构造正规方程的右端向量 A^T b（2x1） b1 = np.dot(A.T, b) # 第三步：解 2x2 线性方程组，得到最小二乘解 print(np.linalg.solve(x, b1))

np.dot(A.T, A)把 4×2 的矩阵压成 2×2，这一步是"降维"的核心：原来 4 个方程互相矛盾，压缩后变成 2 个必然相容的方程。np.dot(A.T, b)把右端向量从 4 维投影到 2 维。最后np.linalg.solve解的是一个方阵系统，不存在无解问题。

参数上要注意：A 的形状必须是 (m, n) 且 m>n，b 的形状是 (m, 1)。如果 b 写成 (m,) 一维数组，np.dot(A.T, b)会得到 (n,) 而不是 (n,1)，solve仍能算但形状语义会变，后续拼接容易翻车。实验里 b 用[[11],[3],[6],[7]]的二维列向量写法，是更稳妥的选择。

2.3 数值稳定性：为什么教科书推荐正规方程，工程上却要小心

正规方程在数学上完美，但在浮点数世界里有个隐患：条件数会被平方。A 的条件数是 κ(A)，AᵀA 的条件数是 κ(A)²。如果 A 本身条件数到了 10⁴，AᵀA 就到了 10⁸，双精度下有效数字会掉一大截。

实验数据规模小、矩阵良态，直接解正规方程没问题，结果 x₁=3.04029304、x₂=1.24175824 也验证了这一点。但常见的实现思路是：数据量大或矩阵病态时，改用 QR 分解或 SVD 求解，numpy 里对应np.linalg.lstsq(A, b, rcond=None)，它内部走 SVD，不显式构造 AᵀA，条件数只放大一次。我一般会这么理解：正规方程是"理解最小二乘的入口"，lstsq 是"生产环境的出口"，两者算的是同一个解，只是数值路径不同。

提示：判断该用哪个，看 A 的列是否近似线性相关。列几乎共线时，AᵀA 接近奇异，solve 可能报错或给出巨大数值，此时必须换 lstsq。

3. 非线性拟合：把 1/y=c₀+c₁/x 掰回线性模型

3.1 变量替换的本质：换的是坐标系，不是模型

题目六给的拟合函数是 y = x/(c₀x+c₁)，对参数 c₀、c₁ 而言这是非线性的，正规方程那套直接套不上。实验报告的处理很巧妙：取倒数，令 Y=1/y、X=1/x，函数变成 Y=c₀+c₁X，瞬间回到线性组合模型。

这一步的关键认知是：变量替换改变的是"数据所在的坐标系"，不是"模型对参数的非线性"。替换后，基函数是 φ₀(x)=1、φ₁(x)=1/x，拟合函数是这两个基函数的线性组合，参数 c₀、c₁ 以一次方出现，正规方程重新适用。

但这里埋着一个容易被忽略的代价：最小化的是 Σ(1/yᵢ − c₀ − c₁/xᵢ)²，而不是原始的 Σ(yᵢ − f(xᵢ))²。两者不等价，倒数变换会给小 y 值更大的权重。实验数据 y 在 0.33 到 0.45 之间，量级接近，这个偏差可以接受；如果 y 跨越几个数量级，变换后的拟合结果可能和原始残差意义下的最优解差很远。这是"能化成线性就化成线性"这条经验背后的取舍。

3.2 正规方程组的组装：求和表怎么变成矩阵

实验报告里列了一张求和表，把 Σ1、Σ(1/xᵢ)、Σ(1/xᵢ²)、Σ(1/yᵢ)、Σ(1/(xᵢyᵢ)) 都算出来，然后拼成 2×2 正规方程组：

import numpy as np # 原始观测数据 x = np.array([1.0, 2.0, 4.0, 5.0]) y = np.array([0.33, 0.40, 0.44, 0.45]) # 变量替换：X = 1/x, Y = 1/y X = 1.0 / x Y = 1.0 / y # 基函数 φ0=1, φ1=1/x，组装正规方程 A^T A c = A^T Y Phi = np.column_stack([np.ones_like(X), X]) # 设计矩阵，形状 (4, 2) ATA = Phi.T @ Phi # 2x2 系数矩阵 ATY = Phi.T @ Y # 2x1 右端向量 c = np.linalg.solve(ATA, ATY) # 解出 c0, c1 print("c0=%.4f, c1=%.4f" % (c[0], c[1]))

np.column_stack把两个基函数在四个采样点上的取值拼成设计矩阵，每一列是一个基函数，每一行是一个样本。Phi.T @ Phi就是求和表里那些 Σ 的矩阵化表达：对角元是各基函数的平方和，非对角元是基函数两两乘积和。这样写比手算求和表更不容易错，也更容易扩展到三个、四个基函数。

代入数据后解得 c₀=2.0143、c₁=1.0093，拟合函数为 y = x/(2.0143x+1.0093)。实验报告里还用了 sklearn 的 LinearRegression 做对照，它内部同样解正规方程（默认走 SVD），系数一致，说明两条路径算的是同一个东西。

3.3 从两个基函数到多项式拟合：设计矩阵的统一视角

把视角拉高一点：直线拟合的设计矩阵是 [1, x]，抛物线是 [1, x, x²]，立方曲线是 [1, x, x², x³]。所谓"选拟合函数类型"，在实现层面就是"选设计矩阵的列"。实验四的石油产量预测正是这个思路的延伸——同一组数据，换不同的列组合，得到不同阶数的多项式，再用残差指标比较谁更合适。理解了设计矩阵，多项式拟合就不再是一堆 polyfit 调用，而是一个统一的线性最小二乘框架。

4. 石油产量预测：残差指标为什么比 R² 更能说明问题

4.1 三种拟合的系数从哪来：polyfit 内部做了什么

实验用np.polyfit(x, y, deg)分别拟合 1、2、3、4 阶多项式。polyfit 内部并不神秘：它构造范德蒙德矩阵 V（第 j 列是 x 的 j 次方），然后解 V c ≈ y 的最小二乘问题，返回的系数按高次到低次排列。所以np.polyfit(x, y, 1)返回 [k, b]，np.polyfit(x, y, 3)返回 [a₃, a₂, a₁, a₀]。

实验报告给出的三组表达式值得逐项看：

三个指标同向递减，立方曲线全面占优。但这里有个数值细节：x 用的是 1994 到 2003 的绝对年份，量级 10³，三次方后到 10⁹，系数被压到 10⁻² 甚至 10⁻⁵ 量级。这不是错误，而是"大数自变量"导致的系数尺度问题。常见的实现思路是先把年份中心化（减去 1994 或均值），系数会回到正常量级，数值稳定性也更好。实验代码里预测 2010 年时用2010-1994做偏移，其实已经隐含了中心化的意识，只是拟合阶段没做。

4.2 残差分析的三个口径：二范数、标准差、绝对均值

实验报告同时算了三个残差指标，它们回答的是不同问题：

• 二范数的平方 Σ(yᵢ−ŷᵢ)²：最小二乘直接优化的目标，对离群点最敏感，一个大残差能主导整个值。
• 残差的标准差：把残差当随机变量看离散程度，扣掉了均值偏移，反映"波动大小"。
• 残差绝对值的均值：L1 口径，对离群点不敏感，反映"平均偏差多大"。

三个口径一起看，能避免单一指标的误导。如果某个模型二范数很小但绝对均值不小，说明它被少数点"绑架"了，大部分点拟合得并不好。实验里立方曲线三项都最低，结论比较稳。

4.3 预测 2010：为什么数值最优不等于结论可信

三种模型对 2010 年的预测分别是 83.38、74.41、99.04。单看残差，立方曲线最好，预测值也最高。但实验报告的分析点出了关键：抛物线预测出下降趋势，与石油产量持续增长的现实不符；立方曲线增长快于直线，符合"增长加速"的直觉。

这里暴露的是外推的经典陷阱：多项式阶数越高，样本区间外的行为越不可控。三次多项式在 1994–2003 区间内拟合最好，但出了这个区间，三次项会主导，曲线可能急剧上翘或下坠。残差指标只衡量"区间内拟合得多好"，完全不衡量"区间外预测得多合理"。所以选模型时，残差是必要条件不是充分条件，还得叠加领域常识做约束。这也是为什么实验要求"讨论结果"而不是"报告最小值"——数值分析课在这里教的其实是建模判断力。

import numpy as np # 年份中心化，避免大数导致的系数尺度问题 year = np.arange(1994, 2004, dtype=float) y = np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777]) t = year - 1994 # 中心化后 t 在 0..9 for deg in (1, 2, 3): coef = np.polyfit(t, y, deg) # 拟合，系数按高次到低次 pred = np.polyval(coef, 2010 - 1994) # 外推到 2010 年 resid = y - np.polyval(coef, t) # 区间内残差 print("deg=%d, 2010预测=%.4f, 残差平方和=%.4f" % (deg, pred, (resid ** 2).sum()))

np.polyval按系数从高到低求值，和 polyfit 的输出顺序严格对应，这是最容易写反的地方。中心化后 t 的量级是 10，三次方只有 10³，系数回到可读范围，外推时的数值误差也小得多。

5. 避坑与排查：最小二乘实验里最容易翻车的五处

5.1 现象：solve 报 "Singular matrix"，但矩阵看着没问题

原因：AᵀA 接近奇异，通常是 A 的列近似线性相关，或者数据量刚好等于未知数个数导致方阵退化。解决：先算np.linalg.cond(A)，条件数超过 10¹⁰ 就改用np.linalg.lstsq(A, b, rcond=None)，它走 SVD，能给出最小范数解而不是直接报错。

5.2 现象：拟合系数和手算求和表对不上

原因：变量替换时把 1/x 和 1/y 的对应关系搞反，或者设计矩阵列顺序和系数顺序不一致。解决：固定"第 i 列对应第 i 个基函数"的约定，组装完打印Phi.shape和Phi[:2]抽查前两行，确认每列确实是预期的基函数取值。

5.3 现象：polyfit 警告 "RankWarning: Polyfit may be poorly conditioned"

原因：用绝对年份做自变量，范德蒙德矩阵列之间量级差好几个数量级，条件数爆炸。解决：先中心化再拟合，t = year - year.mean()，拟合完预测时用同样的偏移量代入，系数立刻回到正常量级。

5.4 现象：残差平方和很小，但预测值离谱

原因：高阶多项式在样本区间外发散，区间内拟合好不代表外推可信。解决：外推前画图看曲线在目标区间附近的走向，同时用低阶模型做交叉验证；如果高阶和低阶预测差距巨大，说明数据不足以支撑高阶外推，应降阶或加约束。

5.5 现象：sklearn 和 numpy 算出的系数有微小差异

原因：两者数值路径不同，LinearRegression 默认走 SVD，polyfit 走 QR 或正规方程，浮点误差量级不同。解决：差异在 10⁻⁶ 以内属正常，不必强行对齐；需要严格复现时统一用同一种求解器，并在文档里注明。

6. 设计取舍：为什么"能化成线性"是一条有代价的捷径

回头看这三个问题，它们共享同一个骨架：把待求参数整理成线性组合，构造设计矩阵，解正规方程。超定方程组是"设计矩阵就是 A 本身"，非线性拟合是"变量替换后重新定义基函数"，多项式预测是"基函数取 x 的各次幂"。这套框架的威力在于统一，代价在于每次"化线性"都偷偷换了优化目标。

倒数变换最小化的是变换后空间的残差，不是原始空间的残差；多项式升阶降低的是区间内残差，不是外推误差；正规方程求解快，但条件数平方。这些取舍没有绝对对错，取决于你要的是"数学上最优"还是"业务上可信"。我一般会这么判断：先问残差是在哪个空间里算的，再问这个空间和真实误差是否同构，最后才看数值大小。

验证方法上，我习惯做两件事：一是用 lstsq 和 solve 各算一遍，结果差超过 10⁻⁶ 就查条件数；二是把拟合曲线和原始散点画在同一张图上，肉眼过一遍残差分布，有没有系统性弯曲、有没有被单点拉偏，图上一眼就能看出来，比盯数字快得多。

从那以后我每次做最小二乘，都强制走一遍"中心化 → 看条件数 → 双求解器对照 → 画残差图"这四步，哪怕数据只有十个点。这套习惯帮我躲过了好几次"数值最优但结论荒谬"的坑。希望帮到你。

本文还有配套的精品资源，点击获取