企业官网建设流程全解析

简介：面向广东工业大学数值分析课程学习者，这份Python实验四上机报告涵盖三个核心专题：超定方程组的最小二乘解、最小二乘拟合以及石油产量预测。题目三要求求解两组具体超定方程组a、b的最小二乘解，报告中推导了正规方程ATAx=ATb，并以numpy.linalg.solve完成求解；题目六给定观察数据x为1、2、4、5，对应y为0.33、0.40、0.44、0.45，构造含基函数的线性组合模型，结合sklearn LinearRegression与matplotlib得到拟合函数并绘制图像；石油产量预测则演示了利用历史数据建立模型、通过最小二乘法优化参数的流程。资源为1个docx文档，压缩包大小仅112KB，内容包括实验题目、算法介绍、完整源代码、结果分析及结论，目录结构清晰，便于对照运行、复习提交或改写扩展。目前已有212人学习下载，适合正在完成同类数值分析实验或备考课程的同学参考。

1. 数值分析实验四的完整复现：超定方程组、最小二乘拟合到产量预测

数值分析这门课的上机实验，做到第四次基本就是和小二乘正面交锋了。这份实验报告覆盖了三件事：求超定方程组的最小二乘解、对一组观察数据做最小二乘拟合，以及用多项式预测石油产量并比较模型效果。整套代码用 Python 的 NumPy、scikit-learn 和 SciPy 写完，直接复制就能跑，结果也全部给出来了，省去大量推导时间。如果你正在做类似的上机作业，或者工作中遇到方程个数比未知数多、需要从观测数据里拟合曲线的问题，这份资源正好是现成的参考。读这篇文章，你会看到正规方程怎么从矩阵乘法里长出来，数据变换怎么做才不翻车，以及多项式预测里藏着的那个外推陷阱。

2. 超定方程组的最小二乘解：正规方程 AᵀAx = Aᵀb 的导出与验证

2.1 方程多了反而无解：矛盾方程组的本质

先看题目给的第一类问题。方程个数多于未知数个数，比如下面这个 4 个方程、2 个未知数的系统：

2x1 + 4x2 = 11 3x1 - 5x2 = 3 1x1 + 2x2 = 6 2x1 + 1x2 = 7

想找到一个 (x1, x2) 同时满足四个等式，绝大多数情况下不存在这样的解。原因很直观：两个未知数只能确定两条直线，多出来的方程如果恰好落在前两条直线的交点上，方程组才有精确解；可现实里的数据来自测量或观测，右边常数项带着误差，交点会被扰动，第三条、第四条等式自然就对不上了。这类方程组叫“矛盾方程组”，数学上称超定方程组。

但工程上不能因为没有精确解就放弃它。这里需要换个目标：不求 Ax = b 严格成立，而是求一个 x，让残差向量 r = Ax - b 的长度尽可能短。这就是最小二乘的含义，名字里的“二乘”指的是平方，也就是最小化残差平方和：

min ||Ax - b||² = min Σ(ai1·x1 + ai2·x2 + ... + ain·xn - bi)²

2.2 正规方程是怎么来的：从投影到可解条件

上面的最小化问题，核心结论是一个定理：x* 是超定方程组 Ax = b 的最小二乘解，当且仅当 x* 满足下面的正规方程：

AᵀAx = Aᵀb

理解它不用死记。A 是 m×n 矩阵，m > n，它的列向量张出一个 n 维子空间。Ax 的取值范围就是 A 的列空间，b 不一定落在这个空间里，所以我们要从列空间里找一个离 b 最近的向量，它就是 b 在列空间上的正交投影。既然是正交投影，残差 b - Ax 必须垂直于列空间的每个基向量，也就是 Aᵀ(b - Ax) = 0，整理一下就是 AᵀAx = Aᵀb。

从代数角度看，AᵀA 是一个 n×n 方阵，只要 A 的列线性无关，AᵀA 就可逆，于是 x 能被唯一解出来。这就是为什么很多教材和工程代码里都用正规方程来求最小二乘解——它把求解问题从一个矩形矩阵变成了方阵，直接套用高斯消元或 LU 分解就能处理。

这里有一个容易忽略的前提：A 的列向量必须线性无关。如果某一列是其他列的线性组合，AᵀA 会奇异，正规方程求出来的解会是无穷多个，数值上也会出现奇怪的结果。后面避坑章节会再提到这一点。

2.3 两个算例的手算验证：矩阵组装与结果核对

题目的第一问给了两个超定方程组，我们在 Python 里把系数矩阵和右端向量组装出来。第一个：

A = [[2, 4], [3, -5], [1, 2], [2, 1]] b = [[11], [3], [6], [7]]

先手算一下正规方程的系数矩阵。AᵀA 是 2×2 矩阵，主对角线上的元素分别是第一列和第二列各自的平方和，副对角线是两列的内积：

AᵀA = [[2²+3²+1²+2², 2·4+3·(-5)+1·2+2·1], [2·4+3·(-5)+1·2+2·1, 4²+(-5)²+2²+1²]] = [[18, -3], [-3, 46]] Aᵀb = [[2·11+3·3+1·6+2·7], [4·11+(-5)·3+2·6+1·7]] = [[35], [46]]

解这个二元一次方程组，得到 x1 = 3.04029304，x2 = 1.24175824。把这两组数代回原方程的第一行看残差：

2×3.04029304 + 4×1.24175824 = 11.048，离 11 差了 0.048 3×3.04029304 - 5×1.24175824 = 2.905，离 3 差了 0.095

第二行残差相对大一些，这符合最小二乘的本质——它不保证每个方程都精确成立，而是让所有方程残差的平方和整体最小。手算这一步很有用，它能让你分辨“程序跑出来的数”和“可信的数”。程序输出不是拿到就直接信，通常我会先用这个规模小的算例手算一遍，确认矩阵组装没弄错，然后再加载更大的数据。

第二个方程组是 5 个方程、3 个未知数：

A1 = [[1, -2, 4], [1, -1, 1], [1, 0, 0], [1, 1, 1], [1, 2, 4]] b2 = [[0], [1], [2], [1], [0]]

A₁ᵀA₁ 是 3×3 矩阵，解出来 x1 = 1.65714286，x2 = 0，x3 = -0.42857143。x2 恰好等于 0，看着像巧合，其实是这个数据集本身的结构决定的——围绕 x = 0 对称的数据点在中间位置没有线性趋势，拟合出的二次曲线顶点就落在轴的中央。

3. 用 NumPy 求解正规方程：np.linalg.solve 完整拆解与验证

3.1 第一道题的代码流程：组装矩阵、一次性解出

报告里给的求解代码不长，核心步骤如下。先看第一个超定方程组：

import numpy as np # 超定方程组 a A = np.array([[2, 4], [3, -5], [1, 2], [2, 1]]) # 4x2 系数矩阵 b = np.array([[11], [3], [6], [7]]) # 4x1 右端向量 # 组装正规方程 A^T A x = A^T b ATA = np.dot(A.T, A) # 等价于 A.T @ A ATb = np.dot(A.T, b) # 等价于 A.T @ b # 求解 x = np.linalg.solve(ATA, ATb) print(x)

这里有两个写法上的习惯值得解释。第一，np.dot(A.T, A) 和 A.T @ A 完全等价，前者在旧代码里常见，后者更符合 Python 3.5 之后的可读性。第二，np.linalg.solve 接收两个参数，第一个是系数方阵，第二个是右端向量，要求第一个参数必须是方阵且可逆。AᵀA 恰好满足这个条件，所以能直接套用。

如果不用 np.linalg.solve，也可以手写高斯消元或者调 np.linalg.lstsq 一步到位。lstsq 在数值稳定性上其实比正规方程更好，因为正规方程会把矩阵条件数平方放大，后面避坑章会展开讲。但对于这种 2×2 或 3×3 的小矩阵，solve 完全够用，输出也直观。

3.2 第二道题的残差分布：五方程三未知数

第二个方程组的求解代码结构完全一致：

# 超定方程组 b A1 = np.array([[1, -2, 4], [1, -1, 1], [1, 0, 0], [1, 1, 1], [1, 2, 4]]) # 5x3 系数矩阵 b2 = np.array([[0], [1], [2], [1], [0]]) # 5x1 右端向量 ATA1 = np.dot(A1.T, A1) ATb1 = np.dot(A1.T, b2) x1 = np.linalg.solve(ATA1, ATb1) print(x1)

运行结果：

[[ 1.65714286] [ 0. ] [-0.42857143]]

一个值得验证的细节是残差的符号分布。把解代进原方程，第一行 1×1.657 + (-2)×0 + 4×(-0.4286) = -0.057，第五行对称地 1×1.657 + 2×0 + 4×(-0.4286) = -0.057，中间的方程残差则分布在一侧。这种对称数据下，二次拟合的线性项系数自然退化为 0，是有数学依据的，不是代码 bug。

要更稳妥地验证解的正确性，可以直接算一下残差平方和，和正规方程目标函数的值相互印证：

residual = A1 @ x1 - b2 cost = (residual ** 2).sum() print("残差平方和:", cost)

正常的结果应该是一个小的正数，表示这个解在最小二乘意义下已经是全局最优。

3.3 solve 与 inv 的区别：为什么不要写 np.linalg.inv(AᵀA) @ Aᵀb

网上很多教程会把求解写成 np.linalg.inv(ATA) @ ATb。这条写法能跑，结果也正确，但在数值计算上属于坏习惯。先看两种写法：

# 写法一：solve（推荐） x1 = np.linalg.solve(ATA, ATb) # 写法二：inv + 矩阵乘法（不推荐） x2 = np.linalg.inv(ATA) @ ATb

两边的差别在于计算路径。np.linalg.solve 内部走的是 LU 分解，先分解再回代，整个过程只需要一次分解、两次三角回代，浮点误差被控制在非常小的范围。np.linalg.inv 则要先把逆矩阵完整求出来，再做一次矩阵乘法，这会多引入一轮舍入误差，而且求逆本身的计算量比 LU 分解更大。

矩阵求逆在数值分析里有一句经常被引用的话：不要显式求逆，除非你真的需要逆矩阵本身。实际项目里几乎没有人需要显式的 AᵀA 逆矩阵，大家需要的是“解”向量。正规方程的系数矩阵有时会接近奇异，矩阵越接近奇异，求逆引入的误差越大，而 LU 分解对同样的矩阵会更宽容一些。这个小细节在 2×2 矩阵上体现不出来，换到 50×50 的病态矩阵上，两条写法得到的结果可能在小数点后第三位就开始分叉。

提示：如果你面对的矩阵规模很大或条件数很高，直接调 np.linalg.lstsq(A, b) 或 scipy.linalg.lstsq 走 QR 或 SVD 路线更稳，正规方程目前更多被用在教学场景和中小规模问题上。

4. 最小二乘拟合的线性化：把 y = x/(c0·x + c1) 变成 1/y = c0 + c1·(1/x)

4.1 题目六的建模：广义多项式与基函数

题目六给了一组观察数据：

x = [1, 2, 4, 5] y = [0.33, 0.40, 0.44, 0.45]

要求找一个形如 y = x/(c0·x + c1) 的函数做最小二乘拟合。这个函数不是常规的线性组合形式，c0 和 c1 都藏在分母里，直接对原始变量做线性回归肯定不对。

处理办法是换元。把等式两边取倒数：

1/y = c0 + c1·(1/x)

令 Y = 1/y，X = 1/x，原问题就变成了标准的线性模型 Y = c0 + c1·X。这就是所谓的“可线性化的非线性拟合”——先做变量变换，把非线性问题变成线性问题，再套最小二乘解出参数，最后变换回来。这个思路在课程里叫广义多项式拟合，基函数取 φ0(x) = 1，φ1(x) = 1/x。

这里的关键坑在于：变换后的数据要全部取倒数。报告的代码很清楚地做了这件事：

import numpy as np from sklearn.linear_model import LinearRegression X = np.array([[1, 0.33], [2, 0.40], [4, 0.44], [5, 0.45]]) # 原始数据，每行是 [x, y] X_n = 1 / X # 对每个元素取倒数 # X_n = [[1, 3.0303], [0.5, 2.5], [0.25, 2.2727], [0.2, 2.2222]] x_train = X_n[:, :1] # 取第一列作为变换后的自变量，即 1/x y_train = X_n[:, 1:2] # 取第二列作为变换后的因变量，即 1/y a = np.array(x_train).reshape(-1) # 转成一维，后面画散点图用 b = np.array(y_train).reshape(-1)

train 变量强调它们此时已经不再是原始的 x 和 y，变换后两个变量之间才具备线性关系。如果你在 fit 时错用了原始 y 而不是 1/y，拟合出来的模型就是冲着错误目标去的，后面再怎么调参数也是白费。

4.2 LinearRegression 拟合：参数怎么读

下面这段是拟合和参数提取：

lineModel = LinearRegression() lineModel.fit(x_train, y_train) y_predict = lineModel.predict(x_train) c1 = lineModel.coef_[0][0] # 斜率，也就是原表达式里的 c1 c0 = lineModel.intercept_[0] # 截距，也就是原表达式里的 c0 print("1/y=%.4f 1/x+%.4f" % (c1, c0))

参数含义要对应清楚：线性模型中 Y = c0 + c1·X，X 是 1/x，Y 是 1/y。sklearn 的 LinearRegression 输出里，coef_ 是斜率数组，intercept_ 是截距。这里 coef_[0][0] 读出来的是 X 对应的系数 c1，intercept_ 则是常数项 c0。名字容易搞混，我最初写这段代码时就把 coef_ 和 intercept_ 的顺序弄反过，最后拼公式时才发现拟合函数和图像完全不对。

运行结果得到：

c0 = 2.0143 c1 = 1.00931

也就是说：

1/y = 2.0143 + 1.00931·(1/x)

变换回原始形式：

y = x / (1.00931 + 2.0143·x)

这里再补一句对结果的直观验证：把 x = 1 代入上面的公式，1/y = 2.0143 + 1.00931 = 3.0236，取倒数得到 y = 0.3307，和原始数据第一个点 0.33 非常接近。四个点里误差最大的点也在 0.01 以内，说明线性化后的拟合效果是合理的。

4.3 作图验证：看拟合曲线而不是只读数

代码里最后有一段画图逻辑：

import matplotlib.pyplot as plt plt.scatter(a, b, c="blue") # 变换空间里的样本点 plt.plot(x_train, y_predict, c="red") # 拟合直线 plt.xlabel('1/x') plt.ylabel('1/y') plt.show()

注意这里画的是变换空间里的图，坐标轴是 1/x 和 1/y，不是原始数据。这样做的好处是能直观看到四个样本点在经过变换后确实落在一条直线上，拟合的线性假设是否合理一眼便知。如果想看原始空间里的拟合曲线，还得把预测值变换回去：

x_plot = np.linspace(0.5, 5.5, 100) y_plot = x_plot / (c1 + c0 * x_plot) plt.plot(x_plot, y_plot, c="red") plt.scatter([1, 2, 4, 5], [0.33, 0.40, 0.44, 0.45], c="blue") plt.show()

这一步很多初学者会漏掉——在变换空间里拟合完之后，直接拿拟合直线的 y 值当最终结果，导致画出来的曲线和原始数据完全对不上。数学上的变换是双向的，代入预测前记得要变回来。做这类拟合任务的经验是永远保留两个视图：一个看变换空间的线性关系是否成立，一个看原始空间的拟合曲线是否贴合数据。

5. 石油产量预测的避坑指南：多项式阶数、坐标平移与外推风险

5.1 三种多项式拟合的残差对比

实验的最后一部分是石油产量预测。数据是 1994 到 2003 十年间世界石油日产量，单位是百万桶每天：

年份: 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 产量: 67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777

要求分别用直线、抛物线、立方曲线去拟合这 10 个点，然后比较残差，并预测 2010 年的产量。报告里给出了残差统计表：

三个指标从不同角度描述的其实是同一件事：线性拟合误差最大，抛物线次之，立方曲线最小。从“哪个模型最好地代表这些数据”这个角度看，答案是明确的——立方曲线。但这里必须先泼一盆冷水：这是在样本区间内的拟合精度比较，不代表三次多项式就是预测 2010 年产量的可靠模型。拟合和预测是两回事，把这个概念混淆，是外推翻车的根本原因。

5.2 现象 → 原因 → 解决：四个必踩的坑

坑一：把年份直接当 x 拟合，系数变成天文数字

如果用原始年份做自变量：

x_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]) np.polyfit(x_year, y, 1)

得到的直线系数大约是这样的形式：y = 0.9693·year - 1865。截距的绝对值接近两千，看着吓人，但代回年份一算结果是对的。

问题在于，高次多项式用原始年份拟合时，系数会变得极其敏感。稍微调整一个数据点，高阶系数就可能在几个数量级之间震荡，这就是所谓的病态问题——polyfit 内部用的是范德蒙德矩阵，当自变量的数值范围很大时，矩阵条件数会急剧恶化。

解决方法是做坐标平移或归一化：让 x = 0 对应起始年份。教科书和报告代码里都采用了这个做法：

x = np.arange(0, 10, dtype='float') # 0 对应 1994 年

这样拟合出来的系数数量级就正常了。预测 2010 年时，对应的 x 是 2010 - 1994 = 16，而不是直接把 2010 塞进函数里。很多人在这一步直接翻车，把平移后的系数套到年份数上，结果完全不对。

坑二：polyfit 出来的多项式和手写函数对不上

报告代码里既定义了手写函数 quadratic_fun、cubic_fun，又调用了 np.polyfit，最后用 leastsq 做一轮精修。细看会发现 polyfit 输出的系数顺序和手写函数接收参数的顺序是反的——np.polyfit 返回从高次到低次的系数列表，而 quadratic_fun 自定义函数接收的顺序是 (k1, k2, b)，其中 k1 是二次项系数。如果直接拿 polyfit 的结果递给自定义函数，就会张冠李戴。

处理办法是写一个小工具函数把系数顺序转换一下，或者干脆统一用 np.poly1d 来生成多项式对象再求值：

p2 = np.poly1d(np.polyfit(x, y, 2)) print(p2)

这样既避免了顺序混淆，又能直接用 p2(16) 做预测，代码简洁得多。代码包里的写法其实已经用了 np.poly1d 做打印输出，但自定义函数和 polyfit 混在一起，初学者很容易在两种调用方式之间弄乱。

坑三：leastsq 的初值全设为 0，高次多项式可能收敛不到理想结果

代码里 param0、param1、param2 初始值全部是 [0, 0, ...] 或 [0, 0, 0, 0, 0]。leastsq 是局部优化算法，初始点选得不好，很容易落入局部极小值。对直线和抛物线来说，初值为 0 尚可通过迭代收敛，对四次多项式这种参数空间复杂的模型就不一定了。

常见做法是先用 np.polyfit 解一次作为初值，再交给 leastsq 微调，形成“解析解 + 数值精修”的组合：

init_guess = np.polyfit(x, y, 3) result = leastsq(dist, init_guess, args=(cubic_fun, x, y))

这样既快速又稳定，polyfit 给出的乘数已经接近全局最优，leastsq 只是在它附近再做一步细调。报告代码里把两条路径分开用了，初值却还保持着全零，这个细节不影响最终演示结果，但如果你自己换一组数据，翻车概率会明显上升。

坑四：预测 2010 年时，不同模型的预测值差距悬殊

三种模型给出的 2010 年预测值如下：

直线预测是稳步增长，抛物线的预测反而出现下降趋势，立方曲线则给出了一个明显偏高的数。结合经济和社会发展的大背景，报告结论认为立方曲线更符合实际。但冷静看，三次多项式在 x = 16 处的外推值高达 99，四次多项式的值会更离谱——多项式在样本区间之外往往是“放飞自我”的，次数越高越夸张。任何外推都要带着怀疑去看，样本数据只覆盖了 10 年，用它预测 16 年后的产量，本身置信度就有限。这个 99.04 不是置信区间的中心，更接近一个数学上的极限趋势外推，真要拿去指导决策，建议同时给出好几套模型的结果做区间参考。

5.3 外推的数学直觉：为什么抛物线会拐弯

抛物线这个模型在 1994 到 2003 区间内拟合效果不错，但二次项的系数是负的，意味着抛物线有一个顶点，过了顶点之后趋势就会反转。2010 年对应的 x = 16 已经越过了顶点，预测值自然掉头向下。石油产量在持续增长几十年后突然开始下降，这个结论和现实逻辑不符，所以报告判定抛物线不适合做外推。

立方曲线没有这个问题，因为三次项系数是正的，x 足够大时增长率还在加大。但“增长越来越快”同样缺乏现实依据，资源类产品不可能永远指数式增长。所以在预测场景里，我的习惯是：先用数据做拟合比较，再回到行业背景做判断，模型输出只是众多参考之一。数值分析课程里关注的是算法本身——残差怎么算、比较什么指标——但拿到结论时要记得模型在样本外的行为只是数学行为，不是客观规律。

6. 把流程固化成习惯：平移坐标 + polyfit 定初值 + leastsq 精修 + 残差核对

这一章给一个可以直接迁移到其他拟合场景的完整流程。以石油产量数据为例，四步走完整套流程后，不会比直接抄代码多花时间，但能避免大多数隐藏问题。

第一步，坐标平移。无论拟合什么时间序列，都把起始时刻映射到 0，这样避免了范德蒙德矩阵条件数爆炸，也让系数有直观含义。

第二步，用 polyfit 快速解析求解一次。多项式拟合是一类特殊情况，polyfit 通过拟合线性方程组直接给出最优解，不需要迭代，用它来确定初值几乎不费成本。

import numpy as np from scipy.optimize import leastsq import matplotlib.pyplot as plt # 年份从 1994 到 2003，对应 x = 0 到 9 x = np.arange(0, 10, 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]) p_deg3 = np.polyfit(x, y, 3) # 三次数值解作为初值

第三步，把 polyfit 的结果作为 leastsq 的初始猜测，做精确到残差平方和的非线性精修。定义残差函数后调用 leastsq：

def cubic_fun(s, x): k1, k2, k3, b = s return k1 * x**3 + k2 * x**2 + k3 * x + b def dist(s, func, x, y): return func(s, x) - y init_guess = np.polyfit(x, y, 3) result = leastsq(dist, init_guess, args=(cubic_fun, x, y))

这样得到的参数在数值上更干净，也避开了 leastsq 收敛到局部极小的风险。注意此时 leastsq 的参数顺序和 polyfit 输出顺序一致，都是从高次到低次。

第四步，输出残差指标并绘图验证。用同一套函数算三个数：二范数的平方、残差标准差、绝对值的均值，和报告的表 2 对照。

y_fit = cubic_fun(result[0], x) ss_res = ((y - y_fit) ** 2).sum() # 二范数的平方 std_res = (y - y_fit).std() # 残差标准差 abs_res = (abs(y - y_fit)).mean() # 绝对值的均值 print("二范数平方: %.4f, 残差标准差: %.4f, 绝对值均值: %.4f" % (ss_res, std_res, abs_res))

这三个指标在样本内的排序关系通常一致：高次多项式残差更小。如果某个指标出现反常，说明拟合过程有 bug——比如初值没传对、多项式次数写错、或者数据里有异常点。这个“一致性检查”是我最常用的排查手段，比盯着曲线看还要快。

顺手还可以做一件很有价值的事：把原始数据画成散点图，再叠加三条拟合曲线。不是在变换空间，而是在原始坐标下，这样才能看出哪条曲线的趋势更贴近数据的真实走势。图一旦画出来，模型预测值之间的差距立刻变得直观，你会亲眼看到二次多项式为什么在末端掉头——它被迫用抛物线的形状去拟合数据，抛物线的拐点就落在数据末尾附近。

从那以后，我每次做拟合的第一件事就是把坐标平移写进脚本，初值一律从 polyfit 拿，leastsq 只做优化不做搜索，最后固定走一遍残差核对和双视图绘图。这套流程让我少掉了许多莫名其妙的拟合结果，希望也能帮到你。

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