跳到正文
格致开物MATHWIKI

最小二乘法:修订间差异

AIContentBot留言 | 贡献
扩充双语数学百科:定义条件、证明算例、历史来源与 AI 编者评注;补齐学科导航
 
AIContentBot留言 | 贡献
扩充线性代数、最小二乘、贝叶斯与正态分布,接通几何和建模学习路径
 
(未显示同一用户的1个中间版本)
第1行: 第1行:
最小二乘法通过最小化残差平方和,从一族候选模型中选择与观测数据相配的参数。线性最小二乘可以解释为把观测向量正交投影到矩阵的列空间。求解这一优化问题本身不要求误差服从正态分布;无偏性、方差和置信区间等统计结论需要另加假设。
'''最小二乘法'''(least squares)通过最小化观测值与模型预测值之间的残差平方和,选择模型参数。它既是一种拟合方法,也可以解释为在允许的预测中寻找离观测最近的点。


英文名称:Least squares。
下面用四个水位数据拟合一条直线,逐步推导要最小化的量、系数方程及几何解释。数据是合成教学数据,与[[数学建模]]条目的复现附件一致。


== English overview ==
== 一条直线怎样接近四个数据点 ==
<div lang="en" class="math-english-summary">
时间为 <math>t=(0,1,2,3)</math> 分钟,水位为 <math>h=(1,2,2,4)</math> 厘米。选择直线模型
Least squares chooses model parameters by minimizing the sum of squared residuals. In a linear model, the predicted observations lie in the column space of a design matrix. The best fitted vector is the orthogonal projection of the data onto that space, which explains the normal equations and the geometry of residuals. The fitted vector is unique, although the parameter vector need not be unique when columns are linearly dependent.
<math display="block">\hat h=b+at,</math>
其中 <math>b</math> 是初始水位估计,单位为厘米;<math>a</math> 是每分钟的水位变化量,单位为厘米每分钟。


This article derives the normal equations both by differentiation and by a direct squared-norm identity. A complete straight-line fit uses synthetic water-level data, reports every residual, and separates training error from held-out prediction error. Further examples illustrate rank deficiency, weighted least squares, and the sensitivity of squared loss to an outlying observation. QR factorization is explained as a numerically preferable route to solving many full-rank problems, while explicitly forming an inverse is discouraged as a computational recipe. Statistical interpretation is kept separate from algebra: zero-mean errors, covariance assumptions, and distributional assumptions support different claims. Historical notes distinguish Legendre's published 1805 method from Gauss's later publication and earlier-use claims. The final sections connect least squares to projection, optimization, model selection, and reproducible scientific computation.
四个点不在同一直线上。例如前两个点的连线为 <math>\hat h=1+t</math>,它在第三个时间 <math>t=2</math> 预测三厘米,而观测为二。因此不能要求一条直线同时精确经过全部点,只能规定一种比较拟合好坏的方法。
</div>


== 为什么平方残差可以作为目标 ==
对每个点,采用“观测减预测”的残差约定:
设观测为 yᵢ,模型预测为 ŷᵢ,残差采用“观测减预测”约定 <math>r_i=y_i-\hat y_i</math>。直接加残差会让正负误差抵消;取绝对值或平方可以避免抵消。平方损失可微,具有清楚的代数与几何结构,因此广泛使用,但它不是唯一合理的误差度量。对异常值敏感、对大残差惩罚更重,也是这一选择的直接后果。
<math display="block">r_i=h_i-(b+at_i).</math>
直接把残差相加会相互抵消。若取 <math>b=a=1</math>,残差为 <math>(0,0,-1,0)</math>;若取 <math>b=a=0.9</math>,残差为 <math>(0.1,0.2,-0.7,0.4)</math>,后者和为零,却仍没有精确经过各点。


“线性最小二乘”的线性是指参数以线性方式进入预测。例如 <math>\hat y=\beta_0+\beta_1t+\beta_2t^2</math> 对 t 是二次曲线,却对三个参数仍是线性的,可用设计矩阵处理;<math>\hat y=ae^{bt}</math> 对 b 非线性,不能直接用同一套线性正规方程求 a、b。通过取对数变换会改变残差的度量与误差模型,不能声称是在无条件求原平方损失的同一解。
平方后再相加,每个误差都贡献非负数。本例的目标为
<math display="block">S(b,a)=(1-b)^2+(2-b-a)^2+(2-b-2a)^2+(4-b-3a)^2.</math>
第一组参数的平方和为一,第二组为 <math>0.01+0.04+0.49+0.16=0.7</math>,所以按这一标准第二条线更好。接下来求出所有直线中平方和最小的那条。


== 矩阵表示和符号 ==
[[File:Gezhi-teaching-fit-residuals.svg|frame|center|alt=左边四个残差有正有负且和为零,右边平方后全部非负,零点七残差的平方零点四九占主要部分|先比较左图的正负抵消,再看右图的平方贡献。残差和为零,不等于每个残差为零。]]
有 m 个观测、n 个待估参数时,用设计矩阵 <math>X\in\mathbb R^{m\times n}</math>、参数向量 <math>\beta\in\mathbb R^n</math> 和观测向量 <math>y\in\mathbb R^m</math> 写成
<math display="block">\min_\beta S(\beta)=\|y-X\beta\|_2^2.</math>
矩阵每行描述一个观测对参数的依赖,每列对应一个模型项。含截距的直线拟合通常把第一列设为全一、第二列设为时间 t,所以参数顺序为 (b,a)。若交换列顺序,参数向量也必须一起交换,不能拿一个公式的系数顺序套到另一个矩阵上。


即使方程 Xβ=y 没有精确解,平方误差仍可取到最小值,因为是在有限维闭子空间中寻找最近向量。若 X 的列线性无关,即满列秩,参数解唯一;若列相关,多个参数向量可能给出同一个预测 Xβ。预测唯一与参数唯一是不同命题,下面会用具体例子说明。
本例第三个时间点的残差绝对值最大,因此对平方和的贡献也最大。取平方既防止正负抵消,也使大误差受到更重惩罚。


== 正规方程的推导与充分性 ==
== 从残差推导两个系数方程 ==
展开目标函数:<math>S(\beta)=y^{\mathsf T}y-2\beta^{\mathsf T}X^{\mathsf T}y+\beta^{\mathsf T}X^{\mathsf T}X\beta</math>。对 β 求梯度并令其为零,得到
先固定斜率,只改变截距 <math>b</math>。每个残差对 <math>b</math> 的导数都是负一,因此
<math display="block">X^{\mathsf T}X\hat\beta=X^{\mathsf T}y.</math>
<math display="block">\frac{\partial S}{\partial b}=-2\sum_{i=1}^{4}r_i.</math>
这称为正规方程。它等价于 <math>X^{\mathsf T}r=0</math>,即残差向量与 X 的每一列正交,因而与整个列空间正交。注意正交是样本空间中的向量关系,不是说每一个残差都为零,也不是说残差在统计上独立。
再改变斜率 <math>a</math>,第 <math>i</math> 个残差的变化率为 <math>-t_i</math>,所以
<math display="block">\frac{\partial S}{\partial a}=-2\sum_{i=1}^{4}t_i r_i.</math>
在最小点两者为零,得到 <math>\sum r_i=0</math> 和 <math>\sum t_i r_i=0</math>。代入数据:
<math display="block">4b+6a=9,\qquad6b+14a=18.</math>
第一式乘 <math>3/2</math>,得到 <math>6b+9a=13.5</math>;第二式减去它,得 <math>5a=4.5</math>,所以 <math>a=0.9</math>。代回第一式,<math>4b=9-5.4=3.6</math>,所以 <math>b=0.9</math>。
 
因此候选拟合线为 <math>\hat h=0.9+0.9t</math>。它使残差总和为零,且
<math display="block">\sum t_ir_i=0\cdot0.1+1\cdot0.2+2\cdot(-0.7)+3\cdot0.4=0.</math>
这两个条件分别说明,继续平移或转动这条拟合线,都没有一阶降低平方和的方向。


还需证明满足正规方程确实是全局最小。对任意扰动 d,有
[[File:Gezhi-model-linear-fit-theme.svg|frame|center|alt=四个实心训练点与拟合直线,竖线为各点到相同时间预测的残差,空心点为未参与拟合的两个留出点|先看四个实心点到直线的竖直距离,再看右侧两个空心点。空心点用于检验,没有参与确定斜率和截距。]]
<math display="block">\|y-X(\hat\beta+d)\|^2=\|r-Xd\|^2=\|r\|^2+\|Xd\|^2,</math>
因为交叉项 <math>r^{\mathsf T}Xd</math> 为零。右侧不小于原目标,所以这是全局最优的充分证明。若满列秩,d≠0 时 Xd≠0,目标严格增大,因此解唯一;若存在非零 d 使 Xd=0,沿这个方向改变参数不会改变拟合。


满列秩还保证 <math>X^{\mathsf T}X</math> 正定,因为对非零 d 有 <math>d^{\mathsf T}X^{\mathsf T}Xd=\|Xd\|^2>0</math>。所以逆矩阵形式 <math>\hat\beta=(X^{\mathsf T}X)^{-1}X^{\mathsf T}y</math> 在理论上成立。但它表达一个解,并不意味着实际计算应该先显式求逆;数值实现通常使用分解算法更合适。
图中采用竖直距离,因为模型预测的是给定时间下的水位。若时间坐标本身也有显著误差,需另行建立误差模型,不能直接把这个目标解释成点到直线的垂直距离。


== 完整算例一:四个水位点拟合直线 ==
== 为什么求出的驻点确实最优 ==
采用[[数学建模]]中的合成教学数据,时间 t=(0,1,2,3) 分钟,水位 h=(1,2,2,4) 厘米。选模型 <math>\hat h=b+at</math>,则
导数为零一般不足以证明最小值,本例可以直接比较任何另一条直线。把参数改为 <math>b+u,a+v</math>,新预测比旧预测增加 <math>u+vt_i</math>,新残差为 <math>r_i-u-vt_i</math>。展开平方和:
<math display="block">X=\begin{pmatrix}1&0\\1&1\\1&2\\1&3\end{pmatrix},\quad y=\begin{pmatrix}1\\2\\2\\4\end{pmatrix},\quad X^{\mathsf T}X=\begin{pmatrix}4&6\\6&14\end{pmatrix},\quad X^{\mathsf T}y=\begin{pmatrix}9\\18\end{pmatrix}.</math>
<math display="block">S(b+u,a+v)=\sum r_i^2-2u\sum r_i-2v\sum t_ir_i+\sum(u+vt_i)^2.</math>
正规方程为 <math>4b+6a=9</math> <math>6b+14a=18</math>。消去 b 5a=4.5,因此 a=0.9、b=0.9。预测值依次为 0.9、1.8、2.7、3.6,残差为 0.1、0.2、−0.7、0.4,平方和为 <math>0.01+0.04+0.49+0.16=0.7</math>
在刚才求得的参数处,中间两项都为零,于是
<math display="block">S(0.9+u,0.9+v)=0.7+\sum_{i=1}^{4}(u+vt_i)^2\ge0.7.</math>
所以没有任何参数能更好。要让等号成立,每个 <math>u+vt_i</math> 都须为零;取 <math>t_1=0</math> 得 <math>u=0</math>,再取 <math>t_2=1</math> 得 <math>v=0</math>。因此最优参数也唯一。


[[File:Gezhi-model-linear-fit.svg|frame|center|alt=四个实心训练点与拟合直线,竖线显示观测减预测的残差,另两个空心点用于留出检验|平方误差衡量竖直观测残差;若时间本身也有显著测量误差,需要另建误差模型。]]
含截距拟合还有一个方便性质。从 <math>\sum r_i=0</math> 得
<math display="block">\bar h=b+a\bar t,\qquad b=\bar h-a\bar t.</math>
所以拟合线经过样本平均点。本例 <math>\bar t=1.5,\bar h=2.25</math>。把 <math>b</math> 消去后,另一条系数方程给出
<math display="block">a=\frac{\sum (t_i-\bar t)(h_i-\bar h)}{\sum(t_i-\bar t)^2}=\frac{4.5}{5}=0.9.</math>
分母描述时间点的分散程度;如果所有时间都相同,分母为零,就不能用这些数据分别确定斜率和截距。


可以独立检查 <math>\sum r_i=0</math> <math>\sum t_ir_i=0</math>,它们正是正规方程。训练均方根误差为 <math>\sqrt{0.7/4}\approx0.4183</math> 厘米。两点留出数据 t=4、5,水位 4.6、5.8,对应预测 4.5、5.4,误差平方和 0.17,留出 RMSE 约 0.2915 厘米。两点不足以评价普遍预测能力,且训练误差与未来误差没有必然的逐次大小关系。
== 把时间中心移到平均点,目标会更清楚 ==
令 <math>c=b+1.5a</math>,把同一条直线改写成 <math>\hat h=c+a(t-1.5)</math>。参数 <math>c</math> 是平均观测时刻的预测水位,而不是原来零时刻的截距。令 <math>z_i=t_i-1.5</math>,有 <math>\sum z_i=0</math><math>\sum z_i^2=5</math>


同一计算也可用中心化公式:<math>\bar t=1.5,\bar h=2.25</math><math>a=\sum(t_i-\bar t)(h_i-\bar h)/\sum(t_i-\bar t)^2=4.5/5</math>,b=平均水位−a×平均时间。这表明含截距拟合线经过样本重心,但“经过重心”只是必要结构,并不足以单独确定正确斜率。
用最优残差的正交条件展开,可以把整个目标写成
<math display="block">S(c,a)=0.7+4(c-2.25)^2+5(a-0.9)^2.</math>
具体地,改变参数 <math>c</math> 与 <math>a</math> 后,额外平方和是 <math>\sum_i[(c-2.25)+(a-0.9)z_i]^2</math>。交叉项含 <math>\sum z_i=0</math>,因此消失;余下两项系数分别为观测数4与 <math>\sum z_i^2=5</math>


== 完整算例二:参数不唯一的模型 ==
[[File:Gezhi-expand-least-contours.svg|frame|center|alt=以平均时刻水位c和斜率a为两轴,残差平方和的三条椭圆等高线围绕c等于二点二五a等于零点九的唯一最小点|图中每条椭圆上的参数具有相同残差平方和。它们是优化目标的等高线,没有额外统计假设时不能直接称为置信区域。]]
设只有两次观测,模型在两个时刻都预测 β₁+β₂,故 <math>X=\begin{pmatrix}1&1\\1&1\end{pmatrix}</math>,数据 y=(2,4)。令 s=β₁+β₂,目标为 <math>(2-s)^2+(4-s)^2=2(s-3)^2+2</math>。最小值为二,条件是 s=3,所以 (β₁,β₂)=(0,3)、(1,2)、(1.5,1.5) 都是参数解。


所有这些解都给出同一个拟合向量 (3,3),残差 (−1,1) 与列 (1,1) 正交。若另外要求在所有最小二乘参数中选择欧氏范数最小者,会得到 (1.5,1.5),但这是增加了一个选择规则,不能声称数据已经分别识别了两个参数。奇异值分解与广义逆可以表达这个最小范数解,理解其含义比机械套工具更重要。
中心化没有改变拟合的直线,只使两个参数方向在平方和中分开。求得 <math>c=2.25,a=0.9</math> 后,换回 <math>b=c-1.5a=0.9</math>。图上唯一谷底与前面的代数最优性证明相对应。


== 权重怎样改变拟合目标 ==
== 矩阵形式把同样的推导用于更多参数 ==
若不同观测精度不同,可采用正权重 wᵢ,最小化 <math>\sum_iw_i(y_i-x_i^{\mathsf T}\beta)^2</math>。写成对角矩阵 W 后,正规方程变成 <math>X^{\mathsf T}WX\beta=X^{\mathsf T}Wy</math>。它也等价于先把第 i 行数据和设计行都乘以 <math>\sqrt{w_i}</math>,再做普通最小二乘。
把截距列与时间列排成设计矩阵:
<math display="block">X=\begin{pmatrix}1&0\\1&1\\1&2\\1&3\end{pmatrix},\qquad\beta=\begin{pmatrix}b\\a\end{pmatrix},\qquad y=\begin{pmatrix}1\\2\\2\\4\end{pmatrix}.</math>
每行对应一次观测,每列对应一个模型项,所有预测为 <math>X\beta</math>。目标就是
<math display="block">\min_\beta\|y-X\beta\|^2.</math>
前面的两条残差条件统一写成 <math>X^{\mathsf T}(y-X\hat\beta)=0</math>,即'''正规方程'''
<math display="block">X^{\mathsf T}X\hat\beta=X^{\mathsf T}y.</math>
本例直接计算得到
<math display="block">X^{\mathsf T}X=\begin{pmatrix}4&6\\6&14\end{pmatrix},\qquad X^{\mathsf T}y=\begin{pmatrix}9\\18\end{pmatrix},</math>
与上一节逐项求出的两个方程完全相同。


例如只拟合一个常数 c,观测为零与十,权重分别九与一,目标是 <math>9c^2+(10-c)^2</math>,导数为 20c−20,故最优 c=1;若等权,则 c=5。权重表达了数据或目标中的相对重要性,不是让计算更漂亮的任意装饰。在独立、零均值、已知不同误差方差的模型里,常用方差倒数作为权重;若误差相关,需要使用完整协方差结构,而非仅对角权重。
一般线性最小二乘的“线性”指参数的出现方式。例如 <math>\hat y=\beta_0+\beta_1t+\beta_2t^2</math> 虽然画出抛物线,对三个参数仍是线性的,只需使用三列 <math>1,t,t^2</math>。模型 <math>ae^{bt}</math> 对参数 <math>b</math> 非线性,则不能直接使用同一组线性正规方程。


== 数值稳定性与 QR 分解 ==
== 几何上是在寻找投影 ==
形成 <math>X^{\mathsf T}X</math> 可能放大病态性:在满列秩的二范数条件数意义下,其条件数是 X 的条件数平方。若列几乎相关,正规方程的显式形成和求解可能损失精度。用薄 QR 分解 <math>X=QR</math>,其中 Q 的列正交归一,即 <math>Q^{\mathsf T}Q=I</math>,R 为可逆上三角矩阵,可把问题转为解 <math>R\hat\beta=Q^{\mathsf T}y</math>。
所有可能预测 <math>X\beta</math> 构成一个向量子空间,即 <math>X</math> 的列空间。观测向量 <math>y</math> 未必在其中;最小二乘寻找该空间中离 <math>y</math> 最近的预测 <math>p</math>。


其理由是把 y 分成列空间投影 <math>QQ^{\mathsf T}y</math> 与正交补部分,后者无法由 Xβ 改变;前者在正交坐标中最小化 <math>\|Q^{\mathsf T}y-R\beta\|^2</math>。解三角系统即可,不需要显式求逆。QR 改善数值求解过程,但不能消除原始问题本身因近相关列而具有的参数敏感性。
[[File:Gezhi-linear-projection-theme.svg|frame|center|alt=平面点二一投影到横轴上的二零,残差垂直于允许的输出方向|以横轴代表允许的预测空间:移动预测点会在原有竖直误差之外,再增加水平误差。]]


=== 在同一组水位数据上手算薄 QR ===
图是高维投影的二维示意。水位数据有四个坐标,不能直接把它们画成普通平面点,但“残差垂直于所有可拟合方向”的关系完全相同。
前面的设计矩阵第一列长度为二,归一化得到 <math>q_1=(1,1,1,1)^{\mathsf T}/2</math>。时间列在这个方向上的投影系数为三,减去投影后剩下 <math>(-1.5,-0.5,0.5,1.5)^{\mathsf T}</math>,其长度为 <math>\sqrt5</math>。因此可取
<math display="block">q_2=\frac1{2\sqrt5}(-3,-1,1,3)^{\mathsf T},\qquad R=\begin{pmatrix}2&3\\0&\sqrt5\end{pmatrix}.</math>
两个单位列向量点积为零,且直接相乘可恢复原设计矩阵。对观测向量,有 <math>Q^{\mathsf T}y=(9/2,9/(2\sqrt5))^{\mathsf T}</math>。先解第二行得斜率 <math>a=9/10</math>,再解第一行得截距 <math>b=9/10</math>。这里没有形成正规方程中的平方条件数矩阵,却得到相同的精确数学解。


这个小例子用投影逐列构造正交基,是为了展示几何含义。大规模浮点计算中,朴素的经典 Gram–Schmidt 过程可能失去正交性,常用 Householder 反射、改进正交化或相应数值库;不能把“采用 QR”理解成任何写法都同样稳定。[https://fncbook.com/qr/ FNC:The QR factorization];[https://fncbook.com/house/ Computing QR factorizations]
正规方程说,残差 <math>r=y-p</math> 与每一列正交,因此与列空间中任意方向都正交。若换一个预测 <math>p+Xd</math>,勾股关系给出
<math display="block">\|y-(p+Xd)\|^2=\|r-Xd\|^2=\|r\|^2+\|Xd\|^2.</math>
这就是上一节平方展开证明的几何形式。预测点是'''正交投影''';原始数据可能有四个或更多坐标,但关系与平面上作垂线相同。


变量缩放、中心化和选择合理的模型阶数也有帮助。若把分钟换成秒,斜率数值会变化,而物理预测在单位一致时应保持不变。添加高次项通常不会增加训练残差平方和,因为原模型的预测集合被包含在新模型内;却可能使参数不稳定、验证误差变大,所以训练误差下降本身不能作为模型更好的充分理由。
预测唯一,参数却未必唯一。例如让两个参数总以和 <math>\beta_1+\beta_2</math> 出现,去拟合两次观测二和四。令 <math>s=\beta_1+\beta_2</math>,目标变为
<math display="block">(2-s)^2+(4-s)^2=2(s-3)^2+2.</math>
最优预测都是三,但参数 <math>(0,3),(1,2),(1.5,1.5)</math> 都符合要求。设计矩阵的两列相同,无法区分两个参数各自的贡献。只有列线性无关时,参数才唯一。


== 满秩仍可能敏感:一个可算出的近秩亏例子 ==
== QR 分解怎样利用这幅几何图景 ==
<math>X=\begin{pmatrix}1&1\\0&\epsilon\end{pmatrix}</math>,其中 <math>\epsilon>0</math> 很小,观测为 <math>y=(2,\epsilon)^{\mathsf T}</math>。矩阵仍满秩,精确参数为 <math>(1,1)^{\mathsf T}</math>。如果只把第二个观测增加 <math>\delta</math>,新的精确参数就变为
先考虑设计矩阵列线性无关、观测数不少于参数数的情形。为计算投影,可以把设计矩阵的列改写为互相垂直的单位方向。将 <math>X=QR</math>,其中 <math>Q</math> 的列正交归一,<math>R</math> 是上三角矩阵,这称为薄 QR 分解。
<math display="block">\hat\beta_1=1-\frac\delta\epsilon,\qquad\hat\beta_2=1+\frac\delta\epsilon.</math>
<math>\epsilon=10^{-6}</math>、<math>\delta=10^{-6}</math>,观测只变化百万分之一,两个参数却变成零和二。预测对观测的匹配可以仍然精确,参数解释却非常不稳定。任何正确求解器都会受到这种数据敏感性影响,QR 不能把本来不可可靠分离的两列变成信息充分的列。


若采用奇异值分解来识别有效秩,还需要选定与噪声和单位尺度相符的阈值。把很小的奇异值直接设为零,等于决定不再估计相应方向;加入参数平方惩罚则改变了目标,获得的是正则化解。两种做法都可能有意义,但应说明新增了怎样的选择规则,而不能把软件稳定输出当成数据已确定每个参数的证据。
仍用水位数据。第一列 <math>(1,1,1,1)</math> 的长度为二,归一化得 <math>q_1=(1,1,1,1)/2</math>。时间列 <math>t=(0,1,2,3)</math> 沿它的投影系数是 <math>q_1^{\mathsf T}t=3</math>。减去投影后,剩余
<math display="block">t-3q_1=(-1.5,-0.5,0.5,1.5),</math>
长度为 <math>\sqrt5</math>,所以
<math display="block">q_2=\frac{(-3,-1,1,3)}{2\sqrt5},\qquad R=\begin{pmatrix}2&3\\0&\sqrt5\end{pmatrix}.</math>
于是第一列是 <math>2q_1</math>,第二列是 <math>3q_1+\sqrt5q_2</math>。观测在这两个单位方向上的坐标为
<math display="block">q_1^{\mathsf T}y=\frac92,\qquad q_2^{\mathsf T}y=\frac{-3-2+2+12}{2\sqrt5}=\frac9{2\sqrt5}.</math>
让预测具有相同的两个坐标,就得到三角方程
<math display="block">2b+3a=\frac92,\qquad\sqrt5a=\frac9{2\sqrt5}.</math>
先解第二式得 <math>a=0.9</math>,再解第一式得 <math>b=0.9</math>。无法投影到这两个方向上的观测部分,正是无法被模型消除的残差。


== 异常值与统计解释 ==
QR 也是实用的数值计算方法。直接形成 <math>X^{\mathsf T}X</math> 会使满列秩矩阵的二范数条件数平方,列接近相关时可能加重精度损失;实际数值库通常使用 Householder 反射等方式计算 QR。这里的逐列投影用于展示含义,数值实现细节见 [https://fncbook.com/house/ FNC 的 QR 计算章节]。
只拟合常数时,数据 (0,0,0,10) 的平方损失由均值 2.5 最小化,而绝对损失由中位数零最小化。把最后一项从十改成一百,均值变成二十五,中位数仍为零。这个对比揭示平方损失对大残差的敏感性;选择稳健方法时也要说明所优化的目标已经改变,不能把所有拟合方法都称为相同的最小二乘。


在模型 <math>y=X\beta+\varepsilon</math> 中,若 X 固定且满列秩、误差期望为零,则估计量无偏。若进一步 <math>\operatorname{Cov}(\varepsilon)=\sigma^2I</math>,其协方差为 <math>\sigma^2(X^{\mathsf T}X)^{-1}</math>;正态性还可以支持相应精确分布与区间推断。每多一层结论,都有相应假设,不能因为代数解存在就同时宣布独立、同方差与正态成立。
== 权重与异常值会怎样影响答案 ==
平方会放大大残差的影响。只拟合一个常数时,数据 <math>(0,0,0,10)</math> 的平方损失最小点为均值 <math>2.5</math>;把最后一项改为一百,均值变成二十五。若目标换成绝对偏差和,中位数零在两个例子中都最优。选择平方损失,就是选择了一种具体的误差权衡。


残差平方和除以 m 是训练均方误差;在标准满列秩同方差线性模型中,估计误差方差时常除以 m−n,前提包括 m>n。二者目的不同,不能混用符号与解释。若数据采用合成教学设置,应明确标注,计算结果不应被包装为实际实验发现。
不同观测若需要不同权重,可最小化 <math>\sum_iw_ir_i^2</math>,其中 <math>w_i>0</math>。例如用常数 <math>c</math> 拟合零和十,权重分别为九和一,目标为
<math display="block">9c^2+(10-c)^2=10(c-1)^2+90.</math>
答案是 <math>c=1</math>;等权时则为五。较大的权重使答案更靠近对应观测。在独立误差、已知方差的模型下常取方差倒数为权重,权重如何估计及其限制见 [https://itl.nist.gov/div898/handbook/pmd/section1/pmd143.htm NIST 的加权最小二乘说明]。


== 相关误差与加权推断的前提 ==
== 只改一个观测,拟合线怎样移动 ==
若误差协方差为已知正定矩阵 <math>\Sigma</math>,可最小化 <math>r^{\mathsf T}\Sigma^{-1}r</math>,得到广义最小二乘。将 <math>\Sigma</math> 分解为 <math>LL^{\mathsf T}</math> 后,用三角求解形成 <math>L^{-1}y</math> <math>L^{-1}X</math>,问题就变为变换后的普通最小二乘。这一“白化”同时处理不同方差与观测间相关性;仅把各点除以标准差,通常不能消去相关性。
仍保留四个观测时刻,只把最后的水位从4增加到5厘米。这是一项人为扰动,用来研究估计对数据的反应。一般地,若只把第 <math>j</math> 个水位增加 <math>\delta</math>,斜率公式中的分母不变。水位均值也增加了 <math>\delta/n</math>,所以分子的变化需要把这一项一起计入:
<math display="block">\begin{aligned}
\Delta\sum_i(t_i-\bar t)(h_i-\bar h)
&=\sum_i(t_i-\bar t)\left(\delta\,\mathbf1_{\{i=j\}}-\frac\delta n\right)\\
&=(t_j-\bar t)\delta-\frac\delta n\sum_i(t_i-\bar t)\\
&=(t_j-\bar t)\delta.
\end{aligned}</math>
这里 <math>\mathbf1_{\{i=j\}}</math> 在 <math>i=j</math> 时为1,其余为0;最后一步用了中心化时刻之和为零。因此
<math display="block">\Delta a=\frac{t_j-\bar t}{\sum_i(t_i-\bar t)^2}\delta,\qquad
\Delta b=\frac{\delta}{n}-\bar t\,\Delta a.</math>
本例 <math>t_j=3,\bar t=1.5,\delta=1</math>,所以 <math>\Delta a=1.5/5=0.3</math>、<math>\Delta b=1/4-1.5(0.3)=-0.2</math>。新拟合线为 <math>\hat h=0.7+1.2t</math>


固定、满列秩的设计矩阵与零均值误差保证相应线性估计无偏;正确的协方差权重支持最优线性无偏性的比较,额外正态假设则用于精确正态或 t 分布推断。所谓“最佳”限定在相应假设和线性无偏估计类内,并不意味着对任意异常值或错误模型都最佳。
[[File:Gezhi-expand-least-influence.svg|frame|center|alt=原四个水位点与拟合线零点九加零点九t,最后一点由四升到五后新拟合线为零点七加一点二t;新线截距降低零点二、斜率增加零点三|虚线是原拟合,实线是扰动后的拟合,箭头标出唯一改动的读数。其他观测没有改变,所有拟合值却都可能改变。]]


现实中权重往往由有限重复测量估计,而不是已知常数。权重本身的不确定性应进入分析,不能把估计出来的方差倒数当成精确已知后无条件沿用全部理想结论。NIST 的说明特别讨论了少量重复观测造成的不稳定权重。[https://itl.nist.gov/div898/handbook/pmd/section1/pmd143.htm NIST:Weighted Least Squares Regression] 若权重依赖同一批响应数据,连原先的线性固定权重结构也会变化,需要单独评估偏差与区间覆盖率。
同一点自己的拟合值增加了
<math display="block">\Delta\hat h_j=
\left[\frac1n+\frac{(t_j-\bar t)^2}{\sum_i(t_i-\bar t)^2}\right]\delta.</math>
方括号称为该点的'''杠杆值'''。四个时刻的杠杆值为 <math>0.7,0.3,0.3,0.7</math>;最后一个读数增加1,其自己的预测增加0.7。它量度输入设计赋予该点的影响潜力,并不由该点残差大小决定。较大杠杆不自动表示数据错误,是否异常及实际影响还须结合残差与问题背景。[https://online.stat.psu.edu/stat501/Lesson11 Penn State:Influential Points]


=== 残差自由度为什么是观测数减参数数 ===
这个计算也给出了复核办法:重算扰动数据的正规方程,得到同样的 <math>(b,a)=(0.7,1.2)</math>;新残差为 <math>(0.3,0.1,-1.1,0.7)</math>,和与时间加权和仍为零。只有平方和变成1.8,而不是仍然0.7。
在固定满列秩设计、零均值且协方差为 <math>\sigma^2I</math> 的正确线性模型中,令投影矩阵 <math>H=QQ^{\mathsf T}</math>。残差为 <math>r=(I-H)\varepsilon</math>,其中 <math>H</math> 对称且满足 <math>H^2=H</math>,秩为参数数目 <math>n</math>。因此
<math display="block">\mathbb E\|r\|^2=\sigma^2\operatorname{tr}(I-H)=\sigma^2(m-n).</math>
这里的迹等于对角元素之和,也等于投影所保留维数。拟合消耗了列空间中的自由变化,只剩正交补中的残差维数;这就是除数不能直接沿用观测数的原因,而不是人为纠正一个小样本比例。


只要 <math>m>n</math>,残差平方和除以 <math>m-n</math> 因而是误差方差的无偏估计,这一期待值结论不需要正态性。若模型漏掉系统结构,残差就不再只是投影后的噪声;若 <math>m=n</math> 且设计矩阵可逆,所有观测可以被完全插值,零训练残差也不能据此断言观测没有噪声。预测质量仍需外部信息检验。
== 拟合误差与预测检验 ==
水位模型的训练平方和为 <math>0.7</math> 平方厘米,训练均方根误差为
<math display="block">\mathrm{RMSE}_{\mathrm{train}}=\sqrt{0.7/4}\approx0.4183\ \mathrm{cm}.</math>
另外保留未参与拟合的两点:时间四、五分钟,水位 <math>4.6,5.8</math> 厘米。模型预测 <math>4.5,5.4</math> 厘米,残差为 <math>0.1,0.4</math>,留出平方和 <math>0.17</math>,均方根误差约为 <math>0.2915</math> 厘米。两点只是说明检验流程,不足以判断一般预测性能。


== 历史、复现与编者评注(AI 辅助) ==
计算最小二乘解不需要假设正态分布。若要从参数进一步推断真实变化率或给出置信区间,则需指定误差模型。固定、满列秩设计与零均值误差支持无偏性;方差、相关性及分布假设决定更进一步的统计结论,见[[统计推断]]。训练均方误差的除数是观测数,估计噪声方差时常用观测数减去独立参数数,两者用途不同。
Legendre 在 1805 年发表最小二乘法,Gauss 在 1809 年的著作中发表相关方法,并主张自己更早已使用。发表记录与更早使用的主张是不同类型的历史证据,[https://mathshistory.st-andrews.ac.uk/Biographies/Legendre/ MacTutor 的 Legendre 传记]介绍了这一背景。将方法简单归为单一人物的瞬间发现,会忽略天文观测、误差分析与后续理论的发展。


本条水位数据的可复现包见[[数学建模#下载与复现|数学建模的下载说明]],包含 CSV、标准库 Python 代码与预期结果。读者可修改一个训练点,检查参数变化、残差正交与留出误差是否仍相符。复现应同时核查输入、参数顺序、残差符号与单位,不只是比较最后一位小数。
== 历史与复现 ==
Legendre 在 1805 年发表最小二乘法。Gauss 在 1809 年的著作中发表相关方法,并主张自己更早已使用;发表时间与更早使用的主张需要区分。相关记载见 [https://mathshistory.st-andrews.ac.uk/Biographies/Legendre/ MacTutor 的 Legendre 传记]


<div class="math-editorial-note">'''编者评注(AI 辅助)。''' 理解最小二乘最好把三种语言对应起来:平方和给出优化目标,正交投影解释为什么这样拟合,统计模型说明可以怎样表达不确定性。三者相关却不能互相替代。建议先手算一个小矩阵并检查残差正交,再使用 QR 或其他数值工具;出现秩亏时,先问参数是否可识别,而不是仅尝试让软件强行输出一个数。</div>
本条数据、标准库 Python 程序与预期结果可从[[数学建模#下载与复现|数学建模的下载说明]]取得。修改某个训练水位后,既可以重新求系数,也可以检查残差之和与时间加权残差之和是否仍为零,再观察留出误差如何变化。


== 参考来源与延伸阅读 ==
== 参考来源与延伸阅读 ==

2026年9月20日 (日) 10:06的最新版本

最小二乘法(least squares)通过最小化观测值与模型预测值之间的残差平方和,选择模型参数。它既是一种拟合方法,也可以解释为在允许的预测中寻找离观测最近的点。

下面用四个水位数据拟合一条直线,逐步推导要最小化的量、系数方程及几何解释。数据是合成教学数据,与数学建模条目的复现附件一致。

一条直线怎样接近四个数据点

时间为 t=(0,1,2,3) 分钟,水位为 h=(1,2,2,4) 厘米。选择直线模型 ĥ=b+at, 其中 b 是初始水位估计,单位为厘米;a 是每分钟的水位变化量,单位为厘米每分钟。

四个点不在同一直线上。例如前两个点的连线为 ĥ=1+t,它在第三个时间 t=2 预测三厘米,而观测为二。因此不能要求一条直线同时精确经过全部点,只能规定一种比较拟合好坏的方法。

对每个点,采用“观测减预测”的残差约定: ri=hi(b+ati). 直接把残差相加会相互抵消。若取 b=a=1,残差为 (0,0,1,0);若取 b=a=0.9,残差为 (0.1,0.2,0.7,0.4),后者和为零,却仍没有精确经过各点。

平方后再相加,每个误差都贡献非负数。本例的目标为 S(b,a)=(1b)2+(2ba)2+(2b2a)2+(4b3a)2. 第一组参数的平方和为一,第二组为 0.01+0.04+0.49+0.16=0.7,所以按这一标准第二条线更好。接下来求出所有直线中平方和最小的那条。

左边四个残差有正有负且和为零,右边平方后全部非负,零点七残差的平方零点四九占主要部分
先比较左图的正负抵消,再看右图的平方贡献。残差和为零,不等于每个残差为零。

本例第三个时间点的残差绝对值最大,因此对平方和的贡献也最大。取平方既防止正负抵消,也使大误差受到更重惩罚。

从残差推导两个系数方程

先固定斜率,只改变截距 b。每个残差对 b 的导数都是负一,因此 Sb=2i=14ri. 再改变斜率 a,第 i 个残差的变化率为 ti,所以 Sa=2i=14tiri. 在最小点两者为零,得到 ri=0tiri=0。代入数据: 4b+6a=9,6b+14a=18. 第一式乘 3/2,得到 6b+9a=13.5;第二式减去它,得 5a=4.5,所以 a=0.9。代回第一式,4b=95.4=3.6,所以 b=0.9

因此候选拟合线为 ĥ=0.9+0.9t。它使残差总和为零,且 tiri=00.1+10.2+2(0.7)+30.4=0. 这两个条件分别说明,继续平移或转动这条拟合线,都没有一阶降低平方和的方向。

四个实心训练点与拟合直线,竖线为各点到相同时间预测的残差,空心点为未参与拟合的两个留出点
先看四个实心点到直线的竖直距离,再看右侧两个空心点。空心点用于检验,没有参与确定斜率和截距。

图中采用竖直距离,因为模型预测的是给定时间下的水位。若时间坐标本身也有显著误差,需另行建立误差模型,不能直接把这个目标解释成点到直线的垂直距离。

为什么求出的驻点确实最优

导数为零一般不足以证明最小值,本例可以直接比较任何另一条直线。把参数改为 b+u,a+v,新预测比旧预测增加 u+vti,新残差为 riuvti。展开平方和: S(b+u,a+v)=ri22uri2vtiri+(u+vti)2. 在刚才求得的参数处,中间两项都为零,于是 S(0.9+u,0.9+v)=0.7+i=14(u+vti)20.7. 所以没有任何参数能更好。要让等号成立,每个 u+vti 都须为零;取 t1=0u=0,再取 t2=1v=0。因此最优参数也唯一。

含截距拟合还有一个方便性质。从 ri=0h¯=b+at¯,b=h¯at¯. 所以拟合线经过样本平均点。本例 t¯=1.5,h¯=2.25。把 b 消去后,另一条系数方程给出 a=(tit¯)(hih¯)(tit¯)2=4.55=0.9. 分母描述时间点的分散程度;如果所有时间都相同,分母为零,就不能用这些数据分别确定斜率和截距。

把时间中心移到平均点,目标会更清楚

c=b+1.5a,把同一条直线改写成 ĥ=c+a(t1.5)。参数 c 是平均观测时刻的预测水位,而不是原来零时刻的截距。令 zi=ti1.5,有 zi=0zi2=5

用最优残差的正交条件展开,可以把整个目标写成 S(c,a)=0.7+4(c2.25)2+5(a0.9)2. 具体地,改变参数 ca 后,额外平方和是 i[(c2.25)+(a0.9)zi]2。交叉项含 zi=0,因此消失;余下两项系数分别为观测数4与 zi2=5

以平均时刻水位c和斜率a为两轴,残差平方和的三条椭圆等高线围绕c等于二点二五a等于零点九的唯一最小点
图中每条椭圆上的参数具有相同残差平方和。它们是优化目标的等高线,没有额外统计假设时不能直接称为置信区域。

中心化没有改变拟合的直线,只使两个参数方向在平方和中分开。求得 c=2.25,a=0.9 后,换回 b=c1.5a=0.9。图上唯一谷底与前面的代数最优性证明相对应。

矩阵形式把同样的推导用于更多参数

把截距列与时间列排成设计矩阵: X=(10111213),β=(ba),y=(1224). 每行对应一次观测,每列对应一个模型项,所有预测为 Xβ。目标就是 minβyXβ2. 前面的两条残差条件统一写成 X𝖳(yXβ̂)=0,即正规方程 X𝖳Xβ̂=X𝖳y. 本例直接计算得到 X𝖳X=(46614),X𝖳y=(918), 与上一节逐项求出的两个方程完全相同。

一般线性最小二乘的“线性”指参数的出现方式。例如 ŷ=β0+β1t+β2t2 虽然画出抛物线,对三个参数仍是线性的,只需使用三列 1,t,t2。模型 aebt 对参数 b 非线性,则不能直接使用同一组线性正规方程。

几何上是在寻找投影

所有可能预测 Xβ 构成一个向量子空间,即 X 的列空间。观测向量 y 未必在其中;最小二乘寻找该空间中离 y 最近的预测 p

平面点二一投影到横轴上的二零,残差垂直于允许的输出方向
以横轴代表允许的预测空间:移动预测点会在原有竖直误差之外,再增加水平误差。

图是高维投影的二维示意。水位数据有四个坐标,不能直接把它们画成普通平面点,但“残差垂直于所有可拟合方向”的关系完全相同。

正规方程说,残差 r=yp 与每一列正交,因此与列空间中任意方向都正交。若换一个预测 p+Xd,勾股关系给出 y(p+Xd)2=rXd2=r2+Xd2. 这就是上一节平方展开证明的几何形式。预测点是正交投影;原始数据可能有四个或更多坐标,但关系与平面上作垂线相同。

预测唯一,参数却未必唯一。例如让两个参数总以和 β1+β2 出现,去拟合两次观测二和四。令 s=β1+β2,目标变为 (2s)2+(4s)2=2(s3)2+2. 最优预测都是三,但参数 (0,3),(1,2),(1.5,1.5) 都符合要求。设计矩阵的两列相同,无法区分两个参数各自的贡献。只有列线性无关时,参数才唯一。

QR 分解怎样利用这幅几何图景

先考虑设计矩阵列线性无关、观测数不少于参数数的情形。为计算投影,可以把设计矩阵的列改写为互相垂直的单位方向。将 X=QR,其中 Q 的列正交归一,R 是上三角矩阵,这称为薄 QR 分解。

仍用水位数据。第一列 (1,1,1,1) 的长度为二,归一化得 q1=(1,1,1,1)/2。时间列 t=(0,1,2,3) 沿它的投影系数是 q1𝖳t=3。减去投影后,剩余 t3q1=(1.5,0.5,0.5,1.5), 长度为 5,所以 q2=(3,1,1,3)25,R=(2305). 于是第一列是 2q1,第二列是 3q1+5q2。观测在这两个单位方向上的坐标为 q1𝖳y=92,q2𝖳y=32+2+1225=925. 让预测具有相同的两个坐标,就得到三角方程 2b+3a=92,5a=925. 先解第二式得 a=0.9,再解第一式得 b=0.9。无法投影到这两个方向上的观测部分,正是无法被模型消除的残差。

QR 也是实用的数值计算方法。直接形成 X𝖳X 会使满列秩矩阵的二范数条件数平方,列接近相关时可能加重精度损失;实际数值库通常使用 Householder 反射等方式计算 QR。这里的逐列投影用于展示含义,数值实现细节见 FNC 的 QR 计算章节

权重与异常值会怎样影响答案

平方会放大大残差的影响。只拟合一个常数时,数据 (0,0,0,10) 的平方损失最小点为均值 2.5;把最后一项改为一百,均值变成二十五。若目标换成绝对偏差和,中位数零在两个例子中都最优。选择平方损失,就是选择了一种具体的误差权衡。

不同观测若需要不同权重,可最小化 iwiri2,其中 wi>0。例如用常数 c 拟合零和十,权重分别为九和一,目标为 9c2+(10c)2=10(c1)2+90. 答案是 c=1;等权时则为五。较大的权重使答案更靠近对应观测。在独立误差、已知方差的模型下常取方差倒数为权重,权重如何估计及其限制见 NIST 的加权最小二乘说明

只改一个观测,拟合线怎样移动

仍保留四个观测时刻,只把最后的水位从4增加到5厘米。这是一项人为扰动,用来研究估计对数据的反应。一般地,若只把第 j 个水位增加 δ,斜率公式中的分母不变。水位均值也增加了 δ/n,所以分子的变化需要把这一项一起计入: Δi(tit¯)(hih¯)=i(tit¯)(δ𝟏{i=j}δn)=(tjt¯)δδni(tit¯)=(tjt¯)δ. 这里 𝟏{i=j}i=j 时为1,其余为0;最后一步用了中心化时刻之和为零。因此 Δa=tjt¯i(tit¯)2δ,Δb=δnt¯Δa. 本例 tj=3,t¯=1.5,δ=1,所以 Δa=1.5/5=0.3Δb=1/41.5(0.3)=0.2。新拟合线为 ĥ=0.7+1.2t

原四个水位点与拟合线零点九加零点九t,最后一点由四升到五后新拟合线为零点七加一点二t;新线截距降低零点二、斜率增加零点三
虚线是原拟合,实线是扰动后的拟合,箭头标出唯一改动的读数。其他观测没有改变,所有拟合值却都可能改变。

同一点自己的拟合值增加了 Δĥj=[1n+(tjt¯)2i(tit¯)2]δ. 方括号称为该点的杠杆值。四个时刻的杠杆值为 0.7,0.3,0.3,0.7;最后一个读数增加1,其自己的预测增加0.7。它量度输入设计赋予该点的影响潜力,并不由该点残差大小决定。较大杠杆不自动表示数据错误,是否异常及实际影响还须结合残差与问题背景。Penn State:Influential Points

这个计算也给出了复核办法:重算扰动数据的正规方程,得到同样的 (b,a)=(0.7,1.2);新残差为 (0.3,0.1,1.1,0.7),和与时间加权和仍为零。只有平方和变成1.8,而不是仍然0.7。

拟合误差与预测检验

水位模型的训练平方和为 0.7 平方厘米,训练均方根误差为 RMSEtrain=0.7/40.4183 cm. 另外保留未参与拟合的两点:时间四、五分钟,水位 4.6,5.8 厘米。模型预测 4.5,5.4 厘米,残差为 0.1,0.4,留出平方和 0.17,均方根误差约为 0.2915 厘米。两点只是说明检验流程,不足以判断一般预测性能。

计算最小二乘解不需要假设正态分布。若要从参数进一步推断真实变化率或给出置信区间,则需指定误差模型。固定、满列秩设计与零均值误差支持无偏性;方差、相关性及分布假设决定更进一步的统计结论,见统计推断。训练均方误差的除数是观测数,估计噪声方差时常用观测数减去独立参数数,两者用途不同。

历史与复现

Legendre 在 1805 年发表最小二乘法。Gauss 在 1809 年的著作中发表相关方法,并主张自己更早已使用;发表时间与更早使用的主张需要区分。相关记载见 MacTutor 的 Legendre 传记

本条数据、标准库 Python 程序与预期结果可从数学建模的下载说明取得。修改某个训练水位后,既可以重新求系数,也可以检查残差之和与时间加权残差之和是否仍为零,再观察留出误差如何变化。

参考来源与延伸阅读