简介广东工业大学数值分析课程实验四的Python实现资源围绕超定方程组最小二乘求解与最小二乘拟合两大核心任务适合正在学习数值分析、科学计算或需要完成同类上机实验的本科生参考。实验报告完整记录了题目三两个超定方程组的最小二乘解推导过程以及题目六基于线性组合模型的拟合函数求解并拓展至石油产量预测应用场景。报告包含算法原理说明、正规方程组建立步骤、具体数值结果与图表分析同时附有基于NumPy、scikit-learn与Matplotlib的Python源代码方便读者对照复现。资源为单个docx格式文档压缩包共1个文件大小112KB内容结构清晰从问题描述、算法介绍到结果讨论层层递进。已有212人在CSDN学习或下载过这份实验报告对理解最小二乘方法在实际数据拟合与预测中的应用具有直接参考价值。1. 数值分析实验四一份把最小二乘从公式推到代码仓库的 Python 报告如果你正在抄数值分析作业或者正被超定方程组、最小二乘拟合、多项式预测三座大山压着这份某高校数值分析课程的实验四报告值得你下载下来逐行跑一遍。它不是一个教学 Demo而是一份“拿到题目、写出算法、用 numpy 和 sklearn 落地、再对结果做残差分析”的完整闭环覆盖了三类高频考点超定方程组的最小二乘解、可线性化的非线性拟合、以及多模型对比预测。整份报告只有三个 Python 脚本但把最小二乘的正规方程解法、线性回归模型、polyfit 与 leastsq 两种拟合路径全串起来了。适合正在修数值分析、需要交上机报告或者想看看别人怎么组织“算法 代码 结果分析”的开发者参考。2. 超定方程组求最小二乘解正规方程与 numpy 的配合2.1 先从理论说起为什么 Axb 无解时反而有“最优解”超定方程组就是方程个数大于未知数个数的线性方程组。比如题目给出的第一组2x1 4x2 113x1 - 5x2 3x1 2x2 62x1 x2 74 个方程、2 个未知数绝大多数情况下不存在一组 x 让四个等式同时成立。这在测量、拟合、反演问题里是常态——数据多了系统就“过定”了。数值分析里处理它的标准思路不是硬解而是找一个让所有方程残差平方和最小的解也就是最小二乘解。关键定理x* 是 Axb 的最小二乘解的充分必要条件是 x* 满足正规方程 A^T A x A^T b。推导不复杂目标函数是 ||Ax - b||²对 x 求导并令梯度为零就得到这个结论。所以代码实现的路径非常固定先算 A^T A 和 A^T b再解一个 n 元线性方程组。注意这里有个“小坑”藏在细节里A^T A 是方阵多数情况下是可逆的但条件数可能很差后面避坑章节我会专门讲。2.2 numpy 实现四行代码求解两个方程组报告里的源码用 numpy 把整个过程写得很干净。我把它拆开讲import numpy as np # 超定方程组 a A np.array([[2, 4], [3, -5], [1, 2], [2, 1]]) # 系数矩阵 A4x2 b np.array([[11], [3], [6], [7]]) # 右端向量 b4x1 # 构造正规方程 ATA np.dot(A.T, A) # A^T A得到 2x2 ATb np.dot(A.T, b) # A^T b得到 2x1 # 求解正规方程 print(np.linalg.solve(ATA, ATb))逻辑说明np.dot(A.T, A) 先做矩阵乘法得到 2x2 的系数矩阵np.dot(A.T, b) 得到 2x1 的右端项最后 np.linalg.solve 用高斯消元解这个 2 元方程组。这里必须注意 A.T 的维度——A 是 4x2转置后是 2x4dot 的结果才是 2x2如果顺序写反numpy 会直接抛维度不匹配的异常。第二个方程组也完全同一个套路只是 A 变成了 5x3 的矩阵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 右端向量 x1 np.dot(A1.T, A1) b3 np.dot(A1.T, b2) print(np.linalg.solve(x1, b3))参数说明A1 的行对应 5 个方程列对应 x1、x2、x3 三个未知数。这个结构很典型——列数等于未知数个数行数等于数据点数。运行结果分别是方程组 a 的最小二乘解x13.04029304x21.24175824方程组 b 的最小二乘解x11.65714286x20x3-0.428571432.3 结果验证别只信 solve 的输出拿到解之后我习惯立刻做一个残差验证这是很多人会跳过的步骤。把解代回原方程组看每个方程的残差x_opt np.linalg.solve(ATA, ATb) # 上一步解出的最优参数 residual np.dot(A, x_opt) - b # 原始方程组的残差向量 print(residual.reshape(-1)) # 观察每个方程的误差逻辑说明residual 是 4x1 的向量每个分量对应一个方程的残差。最小二乘解不保证某个残差为零但会让所有残差的平方和最小。验证的意义在于如果某个残差特别大比如比其他大一个数量级说明原始数据里可能存在异常点或者你的正规方程条件数很差解并不可靠。这是“数值解”和“解析解”之间微妙的地方——计算机给的答案你得自己判断该不该信。3. 最小二乘拟合实战数据变换让非线性回归变成线性回归3.1 题目本质形如 y c0 c1/x 的拟合题目六给了一组观察数据x 取 1、2、4、5对应 y 为 0.33、0.40、0.44、0.45。要求找一个形如 y c0 c1/x 的函数做最小二乘拟合。这里有个关键点y 关于参数 c0、c1 是线性的但关于自变量 x 不是线性的——它是个非线性函数。数值分析里管这类问题叫“可线性化的非线性拟合”处理手法非常经典做变量替换。令 X 1/xY 1/y原函数 y c0 c1/x 两边取倒数就变成 1/y 1/(c0 c1/x)。这里报告的做法是直接把数据点 (x, y) 变换成 (1/x, 1/y)然后对变换后的点做线性拟合。换句话说模型变成 Y a·X b其中 X 1/xY 1/y待求的 a、b 就是线性模型的斜率和截距。这个变换的思路是不改变参数的最小二乘意义但把问题从非线性回归降维成线性回归可以直接套用标准工具。3.2 sklearn 实现LinearRegression 的完整用法import numpy as np from sklearn.linear_model import LinearRegression import matplotlib.pyplot as plt # 原始观察数据[x, y] X np.array([[1, 0.33], [2, 0.40], [4, 0.44], [5, 0.45]]) X_n 1 / X # 对所有元素取倒数变成 [1/x, 1/y] 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) # 转成一维数组方便画图 lineModel LinearRegression() lineModel.fit(x_train, y_train) # 拟合线性模型 y_predict lineModel.predict(x_train) c1 lineModel.coef_[0][0] # 斜率 c0 lineModel.intercept_[0] # 截距 print(1/y%.4f 1/x%.4f % (c1, c0)) # 画图显示 plt.scatter(a, b, cblue) # 原始变换后的散点 plt.plot(x_train, y_predict, cred) # 拟合直线 plt.xlabel(1/x) plt.ylabel(1/y) plt.show()逻辑说明X_n 1 / X 利用 numpy 的广播机制对矩阵每个元素取倒数一步完成两个变量的变换。x_train 和 y_train 保持二维列向量形状这是 sklearn 的 LinearRegression 接口的要求——特征必须是二维数组哪怕只有一列。fit 之后coef_ 是斜率数组intercept_ 是截距两者就是变换后线性模型的参数。运行结果斜率 c11.00931截距 c02.0143。还原到原函数就是 1/y 1.00931/x 2.0143也就是说 y 1 / (1.00931/x 2.0143)。注意这里有个容易搞混的地方——sklearn 拟合的是变换后的空间不是原空间。如果想把参数带回原模型必须做逆变换不能直接拿 c0、c1 当原函数的系数用。3.3 为什么数据变换要小心损失函数的意义变了这里有一个值得反复琢磨的点原始最小二乘目标是最小化 Σ(y_i - f(x_i))²但变换后再拟合最小化的是 Σ(1/y_i - 拟合值)²。两个目标函数不一样得到的参数也不一样。严格来说报告里的做法是一种工程近似——它假设变换后的残差分布更接近高斯实际处理中这种“先变换再拟合”是常见手段但如果数据里 y 接近零1/y 会爆炸变换就翻车了。更好的做法是用 scipy 的 curve_fit 直接拟合原函数后面我会给出替代方案。4. 石油产量预测四种多项式拟合与残差对比4.1 任务描述与数据实验四的数据是 1994 到 2003 十年间全球石油日产量百万桶/天要求分别用直线、抛物线、立方曲线拟合这 10 个数据点写出具体的函数表达式做残差分析比较哪种拟合最好再外推预测 2010 年的产量。数据本身有轻微波动但整体上行从 67.052 涨到 76.777。题目问“哪一种拟合最好”——这不是看谁曲线贴得近而是看谁在“预测”这个任务上更可信。报告里用三个指标衡量二范数的平方即残差平方和、残差的标准差、残差绝对值的均值。这三个指标从不同侧面描述拟合误差二范数平方放大大误差标准差体现离散程度绝对均值给一个直观的平均偏差。4.2 polyfit 与手动定义函数两种写法的对比报告代码里同时出现了两种拟合实现np.polyfit 和 scipy.optimize.leastsq。前者是多项式拟合的快捷方式后者是通用最小二乘参数估计。两种写法各有适用场景import numpy as np import matplotlib.pyplot as plt from scipy.optimize import leastsq x np.arange(0, 10, 1, dtypefloat) # 年份偏移量0 代表 19949 代表 2003 y np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777], dtypefloat) # 写法一polyfit 直接拟合多项式 an np.polyfit(x, y, 1) # 1 次多项式直线 print(2010年预测%.4f % (np.polyval(an, 2010 - 1994))) p1 np.poly1d(an) print(p1) # 打印多项式表达式 # 写法二leastsq 通用拟合 param0 [0, 0] def linear_fun(s, x): k1, b s return k1 * x b def dist(a, fun, x, y): return fun(a, x) - y var leastsq(dist, param0, args(linear_fun, x, y)) print([linear] 二范数%.4f % ((y - linear_fun(var[0], x))**2).sum())逻辑说明polyfit 返回系数数组poly1d 把系数包装成可直接调用的多项式对象leastsq 则需要你手动定义函数形式传入初始参数猜测和残差函数。两者解出来的参数在数值上几乎一致区别在于 polyfit 只适合多项式leastsq 可以拟合任意自定义函数。param0 是初始值这里给 [0, 0] 在简单线性问题上没问题但如果拟合更复杂的模型初始值选不好 leastsq 可能收敛到局部最优。4.3 四种拟合的残差对比谁更好不能只看图报告的表 1、表 2 给出了明确的数值对比拟合方式二范数的平方残差标准差绝对值的均值直线9.53200.97630.8107抛物线6.76800.82270.7044立方曲线5.74340.75780.6915从数值上看立方曲线的三个指标都是最小的拟合精度最高。这符合直观预期——数据本身有波动多项式次数越高越能捕捉局部变化。但注意拟合精度高不等于预测可靠。表 3 的 2010 年预测值暴露了问题拟合方式预测 2010 年石油产量直线83.3823抛物线74.4106立方曲线99.0420抛物线预测值比 2003 年还低出现了下降趋势和实际经济发展趋势明显不符立方曲线给出 99.04增长幅度远大于前面十年的趋势。为什么立方曲线拟合误差最小、外推却最离谱这就是高次多项式外推的经典陷阱——多项式在数据区间外会按照最高次项主导的方向飞速增长次数越高外推越猛。报告结论选择了立方曲线理由是“结合经济社会发展趋势”但从纯数学角度看外推预测本身就带有很强的主观性。4.4 画图技巧坐标偏移的细节报告代码里有一个小技巧值得注意x 用的是 0 到 9 的偏移量但画图时用 x1994 把横轴还原成真实年份。plt.axis([1994, 2003, min(y)-0.5, max(y)0.5]) 设置坐标范围plt.grid(True) 打开网格方便读值。预测 2010 年时用 2010-199416 作为输入这是个容易出错的细节——如果你直接用年份 2010 代入 polyfit 得到的多项式结果会完全不对因为拟合时 x 的基准是偏移量而不是真实年份。5. 常见问题与避坑本实验最容易翻车的五个点5.1 plt.axis 范围设置导致图形空白现象代码设置了 plt.axis([1994, 2003, min(y)-0.5, max(y)0.5])但画出来是空白图。原因leastsq 循环里 plt.plot(x1994, ...) 用的是偏移量加年份而如果你在 plt.show() 之前先调用了 plt.axis坐标轴范围会被锁定而后面的 plot 数据超出范围时就只能看到局部。另一个更隐蔽的原因是axis 函数接收的参数顺序是 [xmin, xmax, ymin, ymax]如果写反了图形会直接错乱。解决先 plot 完所有曲线最后再设置 axis或者干脆不设置让 matplotlib 自动缩放。5.2 中文字体 SimHei 不生效现象plt.title 里的中文显示成方块负号显示成方块。原因matplotlib 默认字体不支持中文负号用的是 Unicode 的减号默认字体渲染不出来。代码里已经给出了标准解法mpl.rcParams[font.sans-serif] [SimHei] 指定黑体mpl.rcParams[axes.unicode_minus] False 让负号用 ASCII 的 -。解决如果换成 Mac 或 LinuxSimHei 可能不存在需要改成 WenQuanYi Zen Hei 或 Noto Sans CJK SC。建议用 plt.rcParams[font.family] sans-serif 加一个系统中文字体列表兜底。5.3 用年份直接拟合多项式预测值错得离谱现象用 np.polyfit(x, y, 1) 拟合时x 直接传真实年份得到的多项式系数大得吓人预测 2010 年输出一个完全不符合量级的值。原因真实年份是 1994 到 2003数值很大多项式拟合时高次项会产生严重的数值不稳定。np.polyfit 内部用的是以 x 的幂为基的范德蒙德矩阵x 太大条件数爆炸。报告中用 x np.arange(0, 10) 做偏移就是为了规避这个问题。解决所有拟合和预测都在偏移后的坐标下进行。如果非要显示真实年份只在画图时加偏移量。5.4 leastsq 初始参数给太差拟合结果不稳定现象换一组初始参数leastsq 拟合出来的结果差别很大甚至报错。原因leastsq 用的是 Levenberg-Marquardt 算法对初始值敏感。题目里的初始值给的是 [0, 0]对直线和抛物线问题足够了但如果拟合更复杂的函数比如带指数项或对数项的模型初始值离真实解太远会收敛到局部极小值。解决先用 polyfit 求一个低阶近似的参数作为初始值再传给 leastsq或者改用 scipy.optimize.curve_fit它对参数绑定的检查更友好报错信息也更明确。5.5 残差分析指标选错对比失去意义现象直接比较不同拟合方式的“误差平方和”得出高次多项式永远更好的结论。原因残差平方和会随着多项式次数升高单调下降——次数越高拟合能力越强9 次多项式甚至能精确穿过所有数据点但预测能力反而最差。这就是过拟合。报告里的三个指标只能衡量“拟合程度”不能衡量“预测能力”所以得到立方曲线最好这个结论本身就有前提限制。解决要做真正的模型对比需要用交叉验证或留一法把数据分成训练集和测试集看测试集上的预测误差。后面改造章节会给一个最小实现。6. 二次开发思路把验证型脚本改成可复用的拟合工具前四个实验脚本的价值在于演示但如果要在自己的数值分析作业或真实项目里复用代码还差点意思——所有数据都硬编码在数组里、模型函数写死、评估逻辑散落在主流程中。我一般会先做一次小重构把“拟合”和“评估”抽成通用函数。核心思路是封装一个拟合器类支持任意多项式次数并自动完成残差指标计算import numpy as np from scipy.optimize import leastsq class PolyFitter: def __init__(self, x, y, degree): self.x np.asarray(x, dtypefloat) self.y np.asarray(y, dtypefloat) self.degree degree self.params None def _model(self, coef, x): return np.polyval(coef, x) def _residual(self, coef): return self._model(coef, self.x) - self.y def fit(self, initNone): if init is None: init np.zeros(self.degree 1) # 初始参数全零 self.params leastsq(self._residual, init)[0] return self def predict(self, x): return np.polyval(self.params, x) def metrics(self): res self._residual(self.params) return { l2_sq: (res**2).sum(), std: res.std(), abs_mean: np.abs(res).mean(), } x 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]) fitter PolyFitter(x - 1994, y, degree3).fit() print(fitter.params) print(fitter.metrics()) print(2010 预测, fitter.predict(2010 - 1994))逻辑说明init接收原始数据fit 里用 leastsq 优化多项式系数predict 用 polyval 做预测metrics 一次性算出三个评估指标。这样改完换新数据、换多项式次数都只需要改一行代码。再进一步如果想把数据源换成文件而不是写死的数组就加一个 csv 读取函数如果想做模型对比就把不同 degree 的 PolyFitter 放进一个循环里分别算指标并输出对比表。最后一件事是验证对石油产量这类时间序列我建议至少做一次留一法交叉验证——每次拿 9 个点拟合预测剩下的 1 个点统计预测误差。你会发现结果和三指标排名很可能不一样。从那以后我每次拿到拟合任务都会强制走一遍“画散点图 → 试不同次数 → 交叉验证”的流程再下结论这个习惯帮我避开了不少“拟合效果好但预测全错”的坑。希望帮到你。本文还有配套的精品资源点击获取