线性参数最小二乘全解析:从法方程建模到精度评估与Python实现

📅 发布时间:2026/9/7 1:29:58
线性参数最小二乘全解析:从法方程建模到精度评估与Python实现 简介这份PPT学习教案系统梳理了线性参数最小二乘法的核心内容围绕测量数据处理与曲线拟合的实际需求展开面向工程、物理、统计类专业学生及数据分析人员既可作课堂辅助教学也适合自主复习。课件从最小二乘原理出发依次讲解线性测量方程组的建立、正规方程组的构建与求解、不等权处理以及标准差和相关系数的估计公式。针对非线性参数问题进一步引入泰勒展开线性化与迭代收敛策略配合电容器电容量测定的完整算例逐步演示线性与非线性两种情形的计算流程与结果差异。资源为单个PPT演示课件共二十四页压缩包约306KB页面以公式推导和算例步骤为主结构紧凑、重点突出适合课堂讲授或自学复习。目前已有60人学习可作为测量数据处理、误差理论与最小二乘法应用的入门及进阶参考资料。 去年我为了给测量方向的本科生讲清楚线性参数的最小二乘法处理专门做了一套PPT学习教案。备课时翻阅了不少教材我发现一个很普遍的问题大多数教材把公式推导推完就结束了学生考试会算但换一个真实场景就不知道误差方程怎么列、法方程怎么用程序求、算完参数为什么还要评估精度。这篇文章我就把这套教案背后的完整思路整理出来聊聊线性参数最小二乘的建模过程、代码实现、精度校验以及在实际数据校核中会用到的递推最小二乘RLS和airPLS这类衍生方法。这套内容适合三类读者一是正在学测量平差、数据分析课的学生二是需要用最小二乘拟合参数却对原理有些模糊的工程师三是想设计同类教学材料的老师。我会从建模直觉出发再到手算验证、Python实现、精度评估最后延伸到动态数据的处理方法尽量用说人话的方式把链路打通。1. 为什么说最小二乘不是取平均的简单升级1.1 一个测量问题的两种处理思路假设我们要测一段距离用同一把尺子量了n次得到n个不完全相同的结果。如果只求这一个参数最合理的做法确实是取算术平均这是最小二乘在单参数、等精度情形下的一个特例。但工程中绝大多数问题不是这样。比如用全站仪测多个点的坐标或者用传感器采集一组数据去拟合直线yabx未知参数从一开始就不止一个。此时每个参数单独取平均的思路就行不通了因为$x_i$是共同参与决定多个观测值的参数之间相互耦合。我们需要一个统一的准则在所有可能的参数组合中找到一组参数使得所有观测残差的平方和最小。这就是最小二乘的出发点。这个准则用数学式写出来很简单设观测值为$L_i$对应求出的计算值为$\hat{L_i}$残差为$v_i \hat{L_i} - L_i$那么最小二乘准则就是$$\sum_{i1}^{n} v_i^2 \min$$1.2 为什么损失函数偏偏选平方这个问题我在课堂上几乎每年都会被问到为什么不用绝对值之和用绝对值看起来也能衡量偏差大小而且更直观。这里有几个层面的原因。首先平方函数处处可导方便求极值。我们最后要做的就是对参数求导并令导数为零绝对值函数在零点处不可导处理起来麻烦得多。其次平方会放大较大残差的权重这恰好符合测量的直觉——如果某个观测值偏离整体趋势很远它很可能是粗差或异常值应该让它在目标函数中承担更大的压力促使最优解主动避开它。再往深一层说当观测误差服从正态分布时最小二乘估计在统计上具有最优性质它是极大似然估计的结果。虽然这个结论要等学到数理统计时才会严格证明但对工程应用来说只要误差中没有明显的系统偏差和粗差最小二乘给出的参数估计就是最优线性无偏估计。这也是我在这套教案里第一个强调的点最小二乘不是凭空规定出来的数学技巧它和误差分布之间有深刻联系。2. 从观测方程到法方程线性参数求解的建模全流程2.1 线性参数问题长什么样所谓线性参数问题指的是观测值和参数之间是线性关系。常见的写法是把每个观测值表示成参数的线性组合$$L_i a_{i1}x_1 a_{i2}x_2 \cdots a_{it}x_t \Delta_i$$其中$x_1,\dots,x_t$是待求的t个参数$a_{i1},\dots,a_{it}$是已知系数由观测条件、几何关系决定$\Delta_i$是观测真误差。之所以单独强调线性是因为只有线性关系才能推导出解析解直接通过矩阵运算得到参数。非线性问题则需要先线性化比如泰勒展开取一次项再按线性最小二乘迭代求解。线性参数的最小二乘是整个测量数据处理的地基后面学到间接平差、条件平差都建立在这套框架上。2.2 误差方程与法方程的推导建模时因为观测值携带误差我们不用上式的等号形式而是写成误差方程$$V AX - L$$其中$V$是$n \times 1$残差向量$A$是$n \times t$已知系数矩阵$X$是$t \times 1$参数向量$L$是$n \times 1$观测值向量平差术语中通常用$L$表示观测向量注意此时$L$不是单个观测值而是所有观测值组成的列向量建模时看上下文即可区分。根据最小二乘准则要求$V^TV \min$。对$X$求导并令导数为零$$\frac{\partial (V^TV)}{\partial X} 2V^T \frac{\partial V}{\partial X} 2(AX-L)^T A 0$$整理后得到法方程$$A^TA X A^T L$$如果$A^TA$可逆参数解为$$X (A^TA)^{-1}A^TL$$这一步推导很简洁但就是这短短的几步把求最小转化成了解线性方程组把优化问题变成了矩阵运算问题。我在教学中特别强调法方程不是又一堆公式它本质上是对多维参数求偏导后得到的联立方程。理解了这个过程编程时心里就非常踏实。2.3 实战建模示例用5个点拟合直线为了让学生真正记住建模流程我习惯在PPT里放一个最经典但又不完全纸上谈兵的例子测了5个点$(x_i, y_i)$认为它们近似分布在一条直线上即$y a bx$求系数$a$和$b$。假设5个点的数据是点号xy11.03.222.05.133.06.844.09.355.010.9这里待求参数$X [a, b]^T$观测值向量$L [3.2, 5.1, 6.8, 9.3, 10.9]^T$系数矩阵$A$的每一行对应一个观测点根据$y a bx$写出$$A \begin{bmatrix} 1 1.0 \ 1 2.0 \ 1 3.0 \ 1 4.0 \ 1 5.0 \end{bmatrix}$$然后用公式$X (A^TA)^{-1}A^TL$求解。算出的$a \approx 1.04$$b \approx 1.99$拟合出来的直线是$y 1.04 1.99x$。这个例子最小的好处是手算能算、程序能验、画图能看。学生亲手计算一次再用代码跑一遍最后画个散点图看拟合效果整个流程闭环了后面再遇到复杂问题就不会觉得矩阵运算是一团迷雾。3. 手把手把法方程变成代码用Python告别手算3.1 核心代码从误差方程到参数求解我在这套教案里用的工具是Python NumPy因为它语法直观、生态丰富学生后续做数据分析也能接着用。先把上面这个例子用代码写出来import numpy as np # 观测数据 x np.array([1.0, 2.0, 3.0, 4.0, 5.0]) y np.array([3.2, 5.1, 6.8, 9.3, 10.9]) # 构造系数矩阵A每行对应一个观测方程: y a b*x A np.column_stack([np.ones_like(x), x]) # 方法一按法方程直接求解 ATA A.T A ATL A.T y result np.linalg.solve(ATA, ATL) a, b result print(fa {a:.4f}, b {b:.4f}) # 方法二直接用最小二乘函数 result2 np.linalg.lstsq(A, y, rcondNone) a2, b2 result2[0] print(flstsq: a {a2:.4f}, b {b2:.4f})值得说明的是第5个点5.0, 10.9相对整体趋势略有下偏但最终拟合结果仍然很好地反映了整体线性趋势这就是最小二乘综合权衡所有观测点之后的效果。3.2 为什么用lstsq而不是直接求逆很多初学者看到公式$X(A^TA)^{-1}A^TL$第一反应是用np.linalg.inv(ATA) ATL求逆。我在教案里特意列了三种办法让学生看到差别np.linalg.solve(ATA, ATL)解线性方程组性能好实际中最常用np.linalg.inv(ATA) ATL公式直观但求逆的数值稳定性通常不如直接求解当$A^TA$接近奇异时问题更明显np.linalg.lstsq(A, y, rcondNone)内部使用奇异值分解SVD对病态矩阵更稳健推荐优先使用。这里补充一个要点在实际工程中$A^TA$的条件数很大时法方程仍然可能产生数值不稳定。尤其当参数个数多、数据量级差异大时直接用lstsq可以绕开很多隐形坑。我在讲课时会对学生说公式写出来是为了理解原理代码里用lstsq是为了安稳睡觉。3.3 代码跑通后必须画图验证代码跑完之后我要求学生必须做一件事把观测点、拟合直线、残差画出来。这个习惯我在实际项目中一直保留因为数字不会主动告诉你拟合出了问题。import matplotlib.pyplot as plt residual y - (a b * x) print(残差:, residual) plt.figure(figsize(6, 4)) plt.scatter(x, y, label观测点) plt.plot(x, a b * x, colorred, label拟合直线) plt.xlabel(x) plt.ylabel(y) plt.legend() plt.title(最小二乘直线拟合) plt.grid(True) plt.show()如果残差呈现明显的喇叭状或波浪形说明模型选择有问题或者数据中存在未考虑的系统误差这时候即便参数算出来了也不能直接采用。画图不是形式主义是结果校验的第一道防线。4. 精度评估算出参数只是第一步4.1 单位权中误差衡量观测值整体精度的尺子最小二乘处理不能只给参数结果还要评估结果靠不靠谱。衡量观测值整体精度的指标叫单位权中误差通常记作$\sigma_0$。用后验残差计算$$\sigma_0 \sqrt{\frac{V^TV}{n - t}}$$其中$n$是观测值个数$t$是必要观测数也就是独立参数的个数$n-t$是多余观测数也叫自由度。为什么分母必须是$n-t$而不是$n$我的解释方式是这样的求$t$个未知参数时体系里已经有$t$个约束被消耗掉了这t个约束决定了无论数据怎样拟合残差都会被吃掉一部分自由度因此真正能反映随机误差波动的独立信息只有$n-t$个。如果少减了参数数量$\sigma_0$会被系统性低估让人误以为精度很高这是平差里常见的坑。对上文的5点直线拟合例子算出残差$V^TV$后单位权中误差约为0.13说明这组观测的离散程度约为0.13个单位整体拟合质量不错。4.2 参数中误差每个参数都有自己的精度圈单位权中误差描述的是整体水平实际中我们更关心某个参数自身的精度。参数向量的协因数阵为$$Q_{XX} (A^TA)^{-1}$$第$i$个参数的中误差为$$\sigma_{x_i} \sigma_0 \sqrt{Q_{ii}}$$其中$Q_{ii}$是协因数阵第$i$行第$i$列的对角元素。继续用上文的示例算出的参数中误差大概在$\sigma_a \approx 0.18$、$\sigma_b \approx 0.05$量级。这意味着截距$a$的不确定性相对较大而斜率$b$估计得很稳定。这个结论其实也有直观解释当$x$的取值范围远离0时截距是通过把直线外推到x0处得到的外推必然会放大不确定性。讲到这里学生通常会啊一声因为数值理解和直观感受终于对上了。这也是我觉得这套教案最值得保留的一个教学点精度评估不是冷冰冰的公式套用它背后有实实在在的几何直觉。5. 从静态到动态遗忘因子递推最小二乘与airPLS5.1 遗忘因子递推最小二乘RLS的数据校核应用场景常规最小二乘是批量处理一次拿到全部数据一次性解出参数。但在实际工程中数据往往是逐个到达的比如传感器每秒钟采集一个数据系统参数还可能缓慢漂移。这时候每来一个新数据都重新算一遍全量最小二乘计算量会越来越大也不利于实时跟踪参数变化。遗忘因子递推最小二乘RLS解决的就是这个问题。它通过递推公式在旧估计值的基础上用新观测值逐步修正参数$$K_k \frac{P_{k-1}h_k^T}{\lambda h_k P_{k-1} h_k^T}$$$$X_k X_{k-1} K_k(y_k - h_k X_{k-1})$$$$P_k \frac{P_{k-1} - K_k h_k P_{k-1}}{\lambda}$$其中$\lambda$是遗忘因子一般取0.95到1之间。$\lambda$越接近1旧数据的权重保留越多参数估计越平稳$\lambda$越小系统响应越快但对噪声也越敏感。实际的滑动权衡需要根据数据时变速度来定。我特别想说的是RLS在数据校核中的一个典型用法新观测值如果与当前模型预测值偏差过大残差$y_k - h_kX_{k-1}$会突然变大机制上说明这条数据很可能有问题。因此RLS在递推估计参数的同时天然附带了一个异常检测能力。实测中可以对残差设置阈值超限则报警人工复核该测点的原始数据这比事后批量处理更能及时发现问题。5.2 自适应迭代加权惩罚最小二乘airPLS在基线校正中的应用airPLS的全称是adaptive iteratively reweighted penalized least squares中文常叫自适应迭代加权惩罚最小二乘它在光谱数据基线校正场景中非常常用。它的目标函数是$$S \sum_{i1}^{n} w_i (y_i - z_i)^2 \lambda \sum_{i2}^{n-1} (\Delta^2 z_i)^2$$其中$z$是待求的基线信号第一项让拟合基线贴近原始数据第二项是平滑惩罚项约束基线不能剧烈起伏$\lambda$控制平滑程度。和普通最小二乘最大的区别在权重$w_i$上在迭代过程中对大于基线的数据点通常是信号峰赋予很小的权重甚至接近0对低于基线的点赋予较大权重这样拟合出来的$z$就是一条位于数据下侧的平滑基线而不是把峰也一起拟合掉。这个方法的巧妙之处在于非对称加权。普通最小二乘把正负残差同等对待而airPLS利用信号峰通常只往一个方向偏离基线这一特点通过迭代更新权重把峰的影响逐渐排除最终得到干净的基线。做拉曼光谱、红外光谱数据分析的人应该对这套思路很有共鸣。我在教案里把它作为最小二乘思想在现代数据处理中的延伸来介绍让学生看到同一个数学框架在不同场景下可以长出完全不同的应用形态。6. 把整套内容做成一份能讲得下去的PPT教案6.1 教学节奏与页面划分说是PPT学习教案那么在组织页面时就不能只是公式堆叠。我最终采用的节奏是这样的前几分钟用多次测量取平均够不够这个问题把学生带入场景然后花整块时间讲建模和误差方程代码演示放在推导之后、精度评估之前最后留10分钟做课堂小练习。课堂小练习我只改数据不换题型比如观测点增多到8个、x坐标变为非等间隔让学生独立完成建模、求解、画图、评精度这样能快速暴露他们是否真的理解了流程。教案中有一页我特意放了一个反面教材把误差方程的系数矩阵A列错导致法方程维度对不上程序直接报维度错误。我会现场演示这个报错再让学生去找问题出在哪。这种故意犯错再纠错的环节比单纯讲正确步骤更能加深印象因为错误本身就是最好的教材。6.2 以往教学中最容易卡壳的几个点根据我以往上课的反馈有四个点学生特别容易卡住。第一个是$L$到底是观测值向量还是某个观测值教材符号不统一时很混淆第二是不理解为什么精度评估要除以$n-t$容易死记第三是代码里矩阵维度搞不清A是$n \times t$结果乘出来维度不匹配第四个是RLS递推公式第一次看时觉得凭空多出来一个增益矩阵$K_k$理解不了来历。应对方法我都写进了教案。对第一个问题统一用观测向量为$L$、单个观测值为$L_i$并留足上下文说明对第二个问题用自由度坍缩的例子解释对第三个问题在代码注释里标出每个矩阵的形状对第四个问题我会先用标量形式的递推平均做类比再升级到矩阵形式。最后再分享一个我自己做教案时的小经验不要把所有推导步骤都放到同一页一页只放一个核心结论。设计时可以按照3分钟讲一页的节奏来控制内容量一环扣一环多留空行方便学生记笔记。这样一份PPT既能用来上课也能发出去给别人自学价值就翻倍了。本文还有配套的精品资源点击获取