简介这份资源是广东工业大学数值分析课程第四次上机实验的完整报告面向正在学习数值分析、需要完成Python上机实验的高校学生。内容围绕超定方程组的最小二乘解、最小二乘拟合以及石油产量预测三个典型问题展开涵盖算法原理推导、正规方程构建、结果分析与完整源代码读者可据此理解最小二乘法的数学表述与Python实现路径。资源包共1个docx文件约112KB以实验报告形式组织包含题目、算法介绍、结果及分析和源代码等模块结构清晰便于对照学习。目前已有213人学习下载适合需要参考实验思路、核对求解过程或补充数值计算实践的学生使用。1. 超定方程组为什么方程比未知数多反而有解四个方程、两个未知数按线性代数的常识这组方程应该无解——四条直线在平面上不可能交于一点。但数值分析给出的答案是它有一个最优解。这个反直觉的结论正是最小二乘法的起点。现实里的测量数据几乎不可能完美满足所有方程传感器有噪声、采样有误差、模型本身是近似。方程个数多于未知数个数恰恰是信息冗余的常态而不是异常。超定方程组的最小二乘解就是在这种冗余里找一组让所有方程残差平方和最小的参数。搞懂它等于搞懂了从数据拟合、线性回归到神经网络损失函数的一整条底层逻辑链。这份实验报告把超定方程组、非线性拟合、石油产量预测三个问题串在一起正好构成一条从解方程到做预测的完整路径。2. 正规方程ATAxATb 是怎么把无解变成有解的2.1 从无解到最近点的几何直觉先把问题翻译成几何语言。Axb 中A 是一个 m×n 的矩阵mn它的列向量张成 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 # 超定方程组 a4 个方程2 个未知数 A np.array([[2, 4], [3, -5], [1, 2], [2, 1]]) b np.array([[11], [3], [6], [7]]) # 第一步构造正规方程的系数矩阵 A^T A2x2 x np.dot(A.T, A) # 第二步构造正规方程的右端向量 A^T b2x1 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) 且 mnb 的形状是 (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, rcondNone)它内部走 SVD不显式构造 AᵀA条件数只放大一次。我一般会这么理解正规方程是理解最小二乘的入口lstsq 是生产环境的出口两者算的是同一个解只是数值路径不同。提示判断该用哪个看 A 的列是否近似线性相关。列几乎共线时AᵀA 接近奇异solve 可能报错或给出巨大数值此时必须换 lstsq。3. 非线性拟合把 1/yc₀c₁/x 掰回线性模型3.1 变量替换的本质换的是坐标系不是模型题目六给的拟合函数是 y x/(c₀xc₁)对参数 c₀、c₁ 而言这是非线性的正规方程那套直接套不上。实验报告的处理很巧妙取倒数令 Y1/y、X1/x函数变成 Yc₀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 # 基函数 φ01, φ11/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.0143x1.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₀]。实验报告给出的三组表达式值得逐项看拟合方式函数表达式二范数平方残差标准差绝对值均值直线0.9693x − 18659.53200.97630.8107抛物线−0.07235x² 290.2x − 2.908e056.76800.82270.7044立方曲线0.01821x³ − 109.3x² 2.185e05x − 1.457e085.74340.75780.6915三个指标同向递减立方曲线全面占优。但这里有个数值细节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, dtypefloat) 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, rcondNone)它走 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 默认走 SVDpolyfit 走 QR 或正规方程浮点误差量级不同。解决差异在 10⁻⁶ 以内属正常不必强行对齐需要严格复现时统一用同一种求解器并在文档里注明。6. 设计取舍为什么能化成线性是一条有代价的捷径回头看这三个问题它们共享同一个骨架把待求参数整理成线性组合构造设计矩阵解正规方程。超定方程组是设计矩阵就是 A 本身非线性拟合是变量替换后重新定义基函数多项式预测是基函数取 x 的各次幂。这套框架的威力在于统一代价在于每次化线性都偷偷换了优化目标。倒数变换最小化的是变换后空间的残差不是原始空间的残差多项式升阶降低的是区间内残差不是外推误差正规方程求解快但条件数平方。这些取舍没有绝对对错取决于你要的是数学上最优还是业务上可信。我一般会这么判断先问残差是在哪个空间里算的再问这个空间和真实误差是否同构最后才看数值大小。验证方法上我习惯做两件事一是用 lstsq 和 solve 各算一遍结果差超过 10⁻⁶ 就查条件数二是把拟合曲线和原始散点画在同一张图上肉眼过一遍残差分布有没有系统性弯曲、有没有被单点拉偏图上一眼就能看出来比盯数字快得多。从那以后我每次做最小二乘都强制走一遍中心化 → 看条件数 → 双求解器对照 → 画残差图这四步哪怕数据只有十个点。这套习惯帮我躲过了好几次数值最优但结论荒谬的坑。希望帮到你。本文还有配套的精品资源点击获取