最小二乘法

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

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

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

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

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

对每个点，采用“观测减预测”的残差约定： $$r_i=h_i-(b+a{t}_i).$$ 直接把残差相加会相互抵消。若取 $b=a=1$，残差为 $(0,0,-1,0)$；若取 $b=a=0.9$，残差为 $(0.1,0.2,-0.7,0.4)$，后者和为零，却仍没有精确经过各点。

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

先比较左图的正负抵消，再看右图的平方贡献。残差和为零，不等于每个残差为零。

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

从残差推导两个系数方程

先固定斜率，只改变截距 $b$。每个残差对 $b$ 的导数都是负一，因此 $$\frac{\partial S}{\partial b}=-2\sum _{i=1}^4{r}_i.$$ 再改变斜率 $a$，第 $i$ 个残差的变化率为 $-t_i$，所以 $$\frac{\partial S}{\partial a}=-2\sum _{i=1}^4{t}_i{r}_i.$$ 在最小点两者为零，得到 $\sum r_i=0$ 和 $\sum t_i{r}_i=0$。代入数据： $$4b+6a=9,6b+14a=18.$$ 第一式乘 $3/2$，得到 $6b+9a=13.5$；第二式减去它，得 $5a=4.5$，所以 $a=0.9$。代回第一式，$4b=9-5.4=3.6$，所以 $b=0.9$。

因此候选拟合线为 $\hat{h}=0.9+0.9t$。它使残差总和为零，且 $$\sum t_i{r}_i=0\cdot 0.1+1\cdot 0.2+2\cdot (-0.7)+3\cdot 0.4=0.$$ 这两个条件分别说明，继续平移或转动这条拟合线，都没有一阶降低平方和的方向。

先看四个实心点到直线的竖直距离，再看右侧两个空心点。空心点用于检验，没有参与确定斜率和截距。

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

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

导数为零一般不足以证明最小值，本例可以直接比较任何另一条直线。把参数改为 $b+u,a+v$，新预测比旧预测增加 $u+v{t}_i$，新残差为 $r_i-u-v{t}_i$。展开平方和： $$S(b+u,a+v)=\sum r_i^2-2u\sum r_i-2v\sum t_i{r}_i+\sum (u+v{t}_i{)}^2.$$ 在刚才求得的参数处，中间两项都为零，于是 $$S(0.9+u,0.9+v)=0.7+\sum _{i=1}^4(u+v{t}_i{)}^2\ge \mathrm{0.7.}$$ 所以没有任何参数能更好。要让等号成立，每个 $u+v{t}_i$ 都须为零；取 $t_1=0$ 得 $u=0$，再取 $t_2=1$ 得 $v=0$。因此最优参数也唯一。

含截距拟合还有一个方便性质。从 $\sum r_i=0$ 得 $$\overline{h}=b+a\overline{t},b=\overline{h}-a\overline{t}.$$ 所以拟合线经过样本平均点。本例 $\overline{t}=1.5,\overline{h}=2.25$。把 $b$ 消去后，另一条系数方程给出 $$a=\frac{\sum (t_i-\overline{t})(h_i-\overline{h})}{\sum (t_i-\overline{t}{)}^2}=\frac{4.5}{5}=\mathrm{0.9.}$$ 分母描述时间点的分散程度；如果所有时间都相同，分母为零，就不能用这些数据分别确定斜率和截距。

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

令 $c=b+1.5a$，把同一条直线改写成 $\hat{h}=c+a(t-1.5)$。参数 $c$ 是平均观测时刻的预测水位，而不是原来零时刻的截距。令 $z_i=t_i-1.5$，有 $\sum z_i=0$、$\sum z_i^2=5$。

用最优残差的正交条件展开，可以把整个目标写成 $$S(c,a)=0.7+4(c-2.25{)}^2+5(a-0.9{)}^2.$$ 具体地，改变参数 $c$ 与 $a$ 后，额外平方和是 $\sum _i[(c-2.25)+(a-0.9)z_i{]}^2$。交叉项含 $\sum z_i=0$，因此消失；余下两项系数分别为观测数4与 $\sum z_i^2=5$。

图中每条椭圆上的参数具有相同残差平方和。它们是优化目标的等高线，没有额外统计假设时不能直接称为置信区域。

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

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

把截距列与时间列排成设计矩阵： $$X=\left(\begin{array}{cc}1& 0\\ 1& 1\\ 1& 2\\ 1& 3\end{array}\right),\beta =\left(\begin{array}{c}b\\ a\end{array}\right),y=\left(\begin{array}{c}1\\ 2\\ 2\\ 4\end{array}\right).$$ 每行对应一次观测，每列对应一个模型项，所有预测为 $X\beta$。目标就是 $$\underset{\beta }{min}\Vert y-X\beta {\Vert }^2.$$ 前面的两条残差条件统一写成 $X^{𝖳}(y-X\hat{\beta })=0$，即正规方程 $$X^{𝖳}X\hat{\beta }=X^{𝖳}y.$$ 本例直接计算得到 $$X^{𝖳}X=\left(\begin{array}{cc}4& 6\\ 6& 14\end{array}\right),X^{𝖳}y=\left(\begin{array}{c}9\\ 18\end{array}\right),$$ 与上一节逐项求出的两个方程完全相同。

一般线性最小二乘的“线性”指参数的出现方式。例如 $\hat{y}={\beta }_0+{\beta }_{1}t+{\beta }_2{t}^2$ 虽然画出抛物线，对三个参数仍是线性的，只需使用三列 $1,t,t^2$。模型 $a{e}^{bt}$ 对参数 $b$ 非线性，则不能直接使用同一组线性正规方程。

几何上是在寻找投影

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

以横轴代表允许的预测空间：移动预测点会在原有竖直误差之外，再增加水平误差。

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

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

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

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

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

仍用水位数据。第一列 $(1,1,1,1)$ 的长度为二，归一化得 $q_1=(1,1,1,1)/2$。时间列 $t=(0,1,2,3)$ 沿它的投影系数是 $q_1^{𝖳}t=3$。减去投影后，剩余 $$t-3q_1=(-1.5,-0.5,0.5,1.5),$$ 长度为 $\sqrt{5}$，所以 $$q_2=\frac{(-3,-1,1,3)}{2\sqrt{5}},R=\left(\begin{array}{cc}2& 3\\ 0& \sqrt{5}\end{array}\right).$$ 于是第一列是 $2q_1$，第二列是 $3q_1+\sqrt{5}{q}_2$。观测在这两个单位方向上的坐标为 $$q_1^{𝖳}y=\frac{9}{2},q_2^{𝖳}y=\frac{-3-2+2+12}{2\sqrt{5}}=\frac{9}{2\sqrt{5}}.$$ 让预测具有相同的两个坐标，就得到三角方程 $$2b+3a=\frac{9}{2},\sqrt{5}a=\frac{9}{2\sqrt{5}}.$$ 先解第二式得 $a=0.9$，再解第一式得 $b=0.9$。无法投影到这两个方向上的观测部分，正是无法被模型消除的残差。

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

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

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

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

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

仍保留四个观测时刻，只把最后的水位从4增加到5厘米。这是一项人为扰动，用来研究估计对数据的反应。一般地，若只把第 $j$ 个水位增加 $\delta$，斜率公式中的分母不变。水位均值也增加了 $\delta /n$，所以分子的变化需要把这一项一起计入： $$\begin{array}{cc}\Delta \sum _i(t_i-\overline{t})(h_i-\overline{h})& =\sum _i(t_i-\overline{t})(\delta {\mathrm{𝟏}}_{\{i=j\}}-\frac{\delta }{n})\\ & =(t_j-\overline{t})\delta -\frac{\delta }{n}\sum _i(t_i-\overline{t})\\ & =(t_j-\overline{t})\delta .\end{array}$$ 这里 ${\mathrm{𝟏}}_{\{i=j\}}$ 在 $i=j$ 时为1，其余为0；最后一步用了中心化时刻之和为零。因此 $$\Delta a=\frac{t_j-\overline{t}}{\sum _i(t_i-\overline{t}{)}^2}\delta ,\Delta b=\frac{\delta }{n}-\overline{t}\Delta a.$$ 本例 $t_j=3,\overline{t}=1.5,\delta =1$，所以 $\Delta a=1.5/5=0.3$、$\Delta b=1/4-1.5(0.3)=-0.2$。新拟合线为 $\hat{h}=0.7+1.2t$。

虚线是原拟合，实线是扰动后的拟合，箭头标出唯一改动的读数。其他观测没有改变，所有拟合值却都可能改变。

同一点自己的拟合值增加了 $$\Delta {\hat{h}}_j=[\frac{1}{n}+\frac{(t_j-\overline{t}{)}^2}{\sum _i(t_i-\overline{t}{)}^2}]\delta .$$ 方括号称为该点的杠杆值。四个时刻的杠杆值为 $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$ 平方厘米，训练均方根误差为 $${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=\sqrt{0.7/4}\approx 0.4183\mathrm{c}\mathrm{m}.$$ 另外保留未参与拟合的两点：时间四、五分钟，水位 $4.6,5.8$ 厘米。模型预测 $4.5,5.4$ 厘米，残差为 $0.1,0.4$，留出平方和 $0.17$，均方根误差约为 $0.2915$ 厘米。两点只是说明检验流程，不足以判断一般预测性能。

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

历史与复现

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

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

参考来源与延伸阅读

• Boyd 与 Vandenberghe：Introduction to Applied Linear Algebra，最小二乘、QR、模型拟合。
• Driscoll 与 Braun：The QR factorization、Computing QR factorizations：薄 QR 与稳定实现。
• NIST：Weighted Least Squares Regression：权重估计与统计条件。
• MacTutor：Legendre，最小二乘发表与优先权背景。
• 先修：矩阵、向量空间、导数；相关：优化、统计推断、数学建模。