四个观测点是 $(-1,0),(0,1),(1,1),(2,4)$。如果只允许用一条直线描述它们，方程已经多于未知数，而且四点并不共线。上一章的消元法不能让这些矛盾同时消失。

我们得换一个问题：允许直线与数据有差距，但怎样分配这些差距才算合适？本章采用残差平方和作为标准。这个标准会导出正规方程，也会让我们看见：描述最优解的方程，与适合在计算机上求它的算法，是两件需要分别检查的事。

把拟合写成矩阵

设直线为 $y=c_0+c_1t$。在四个位置代入，得到

$Ac$ 是四个预测值组成的向量。我们沿用上一章的减法顺序，记

因此，数据点在直线上方时，该点的残差为正。这里度量的是固定 $t_i$ 后的纵向差值，不是点到直线的最短几何距离。若横坐标也有不可忽略的测量误差，需要另立模型，不能悄悄把两种距离混在一起。

矩阵的列告诉我们可以调整哪些形状：常数列使全部预测一起上移或下移，第二列使预测按 $t_i$ 的大小改变。若用已知函数 $\phi_1,\ldots,\phi_n$ 拟合，只需取

“线性”说的是对待求系数 $c_j$ 线性。比如 $c_1\sin t+c_2\cos t$ 属于这里的范围； 的第二个未知量藏在正弦里，就不能直接这样处理。

下文令 $A\in\mathbb R^{m\times n}$、$m\ge n$，通常 $m>n$。所有长度和条件数都采用 2-范数；明确写别的范数时才例外。最小化 $\|b-Ac\|_2$ 与最小化它的平方有相同的解，因为平方函数在非负数上严格递增。

同一条最小二乘直线在四处的预测值，与实际观测逐项相减。第三个残差为负，表示预测偏高；平方和把这些偏差合在一起，得到 1.8。

最优残差为什么垂直于列空间

假设目前的系数是 $c$，我们试着把它改成 $c+h$。预测增加 $Ah$，残差变成 $r-Ah$。把长度平方展开：

这个式子同时包含了必要条件和充分条件。

如果 $c$ 已经最优，任取一个方向 $d$，令 $h=td$，则

若 $d^TA^Tr\ne0$，就能选一个与它同号、足够小的 $t$，使一次项压过二次项，右边小于零，与最优矛盾。所以对所有 $d$ 都必须有 $d^T{A}^{T}r=$，即 。

反过来，若 $A^Tr=0$，展开式中的交叉项消失，于是任意 $h$ 都满足

这已经证明 $c$ 是全局最小解。因此

最后一个系统叫正规方程。这里的“正规”与垂直关系相连：$A^Tr=0$ 表示残差与每一列正交，也就与这些列生成的整个列空间正交。注意图上的垂直发生在观测向量所在的 $\mathbb R^m$ 中；它不表示前面散点图中的纵向残差线段垂直于拟合直线。

正交投影的二维示意：预测向量 p 落在列空间里，残差 r 从预测端点指向观测端点，并与列空间垂直。图只表示向量关系，不对应正文四维观测的坐标。

满列秩决定系数能否唯一

刚才的推导没有使用满列秩。秩条件出现在唯一性上。

若 $A$ 满列秩，非零 $h$ 一定满足 $Ah\ne0$。在最优点处，$\Phi(c+h)=\Phi(c)+\|Ah\|_2^2$ 便严格变大，最优系数只能有一组。此时

所以 $A^TA$ 对称正定，正规方程可唯一求解。

若列相关，就有非零 $h$ 使 $Ah=0$。一旦 $c$ 最优，$c+h$ 也给出完全相同的预测。比如两列都是 $(1,2)^T$ 时，预测只依赖 ；数据无法区分两个系数各自的贡献。

最优预测向量仍然唯一。取列空间的一组标准正交基 $q_1,\ldots,q_k$，把 $b$ 在这些方向上的分量相加，

剩下的 $b-p$ 与每个基向量正交。对列空间任意 $y$，勾股分解给出

因此 $y=p$ 是唯一最近点；由于 $p$ 在列空间里，总能找到系数使 $Ac=p$。这既解释了解的存在，也说明“系数不唯一”不等于“拟合值不确定”。

把四点直线完整算一遍

对于开头的数据，

第一行来自所有残差之和为零，第二行来自残差与 $t$ 的内积为零。正规方程是

第一式减去第二式的两倍，可以求得 $c_1=6/5$，再代回得 $c_0=9/10$。拟合直线是

逐点检查比只看两个系数更有用：

残差平方和为

它并不是零，但

我们已找到所选直线模型的最优解。没有必要为了“消灭残差”，硬把数据改掉或把曲线拐过每一个点。

在下面的实验中，先自己调截距和斜率，让平方和下降。试着只消去残差之和，看看残差与 $t$ 的内积是否也已经为零；通常还没有。点击拟合后，查看两个正交条件是否同时满足，再上下拖动一个观测点，观察整条直线怎样回应。

残差平方和把大偏差放大后计入目标，所以一个明显偏离的数据点可能拉动全部系数。实验能显示这种影响，却不能判断那个点究竟是抄错了，还是真实而重要的观测。删除数据需要额外证据。

QR：换一组正交坐标

正规方程便于推导。不过，显式构造 $A^TA$ 会把原矩阵中很小的差异进一步压小。我们先建立一条不必形成这个乘积的计算路径，再讨论差异在哪一步丢失。

设 $A$ 满列秩。完整 QR 分解写成

其中 $Q$ 是 $m\times m$ 正交矩阵，$R$ 是可逆的 $n\times n$ 上三角矩阵，$Q_1$ 有 列， 有 列。去掉乘零的部分，就得到薄 QR：

$Q_1$ 的列是列空间的一组标准正交基。它通常不是方阵，因此不能写 $Q_1Q_1^T=I_m$。

正交方阵保持长度，因为

令 $Q^Tb=(d_1,d_2)^T$，前 $n$ 项组成 ，余下项组成 。把整个残差一起变换：

下半部分没有 $c$，调整任何系数都改变不了它。上半部分则可以通过回代令其为零：

这就解释了 QR 求解法。若只保存薄 QR，也能解出系数；但不能声称 $\|Q_1^T(b-Ac)\|_2=\|b-Ac\|_2$。在最优点，前者恰好为零，后者往往不为零。

完整正交变换把残差分成两段：上方 n 个分量可通过回代消去，下方 m−n 个分量不受系数 c 影响。两段长度的平方相加，仍是原残差的平方长度。

投影矩阵不是单位矩阵

由 $Rc=Q_1^Tb$ 得到

这个矩阵满足

所以投影一次之后再投影，位置不再变化。残差是 $(I-P)b$，而

如果观测只增加一个垂直于列空间的向量，系数不会改变，残差却会改变。若观测沿列空间移动，则可以完全通过调整系数吸收。这两个方向的差别，是判断参数敏感性时很有用的线索。

同一条直线的薄 QR

开头四个 $t_i$ 的平均值是 $1/2$。将第二列减去第一列的 $1/2$ 倍，得到

$s$ 的分量和为零，故它与常数列正交；长度是 $\sqrt5$。于是可选

检查第二列：$t=q_1+\sqrt5q_2$。再算

回代得到 $\sqrt5c_1=6/\sqrt5$，以及 ，仍是 。

把时间原点改到 $t=1/2$，只是换一种参数表达。以 $s=t-1/2$ 为变量时，直线写成 $1.5+1.2s$；预测向量与残差不变。中心化能减少这两列之间的相关性，却不会凭空增加观测信息。

用反射逐列制造零

刚才的 QR 特别容易手算，因为我们看出了两列的正交关系。一般矩阵需要一个可重复的构造。Householder 反射每次处理当前列的一段向量，同时保持长度。

取单位向量 $v$，定义

任意向量 $y$ 可拆成沿 $v$ 的分量 $(v^Ty)v$ 和垂直于 $v$ 的分量。乘以 $H$ 时，前者反号，后者不变。因此它是关于垂直于 $v$ 的超平面的反射。

代数检查也很短：

于是 $H^TH=I$。计算 $Hy$ 时无需造出整个矩阵，只要做一次内积和一次向量更新：

把当前列对准坐标轴

给定非零向量 $z$，希望找到 $H$，使 $Hz$ 只有第一项非零。保持长度决定了第一项只能是 $\pm\|z\|_2$。这两个符号在精确数学中都可以用，计算时却应该认真选择。

记 $\rho=\|z\|_2$，并约定

这里明确规定 $z_1=0$ 时取 $s=1$。如果 $z=0$，已经无需消元，直接跳过，不去除以 $\|w\|_2=0$。

为什么这个构造正确？由于 $\rho^2=z^Tz$，

因此

符号选择还保证

不会在第一项里用两个很接近的同号数相减。比如 $z=(1,\varepsilon)^T$，若执意反射到正轴，构造中会出现 $1-\sqrt{1+\varepsilon^2}$；当 很小时，这个小差不容易准确计算。取负轴后变成相加。

当首分量为正时，把目标选在负半轴会让构造向量的首分量成为两个正量之和；若选正半轴，就会相减两个接近的数。两种目标在精确算术下都可以反射得到。

这不意味着可以忽略其他浮点问题。很大或很小的数据还需要稳定的范数计算与缩放；实际程序应采用成熟数值库。本章的手算与小规模实验用来弄清反射的作用和操作顺序。

后续反射只动剩下的行

第一步作用于整列。第二步只处理第 $2$ 行以下，第三步只处理第 $3$ 行以下。将一个小反射补成大矩阵，可以写成

先前列的活动部分已经全是零，小反射乘这些零仍得到零；前面的行则被单位矩阵保留。因此我们不会破坏已经得到的三角结构。

全部操作给出

左侧乘积是 $Q^T$，而 $Q=\widetilde H_1\widetilde H_2\cdots\widetilde H_n$。每个反射自己是对称的，不代表多个反射可以任意换顺序。

求最小二乘解时，把同样的反射按同样顺序作用于 $b$，就得到了 $Q^Tb$。不必显式储存完整 $Q$；保留反射向量，便可以重复处理不同右端。

一次完整的 Householder 求解

这次取

我们把两次反射、右端变化和最后的残差放在一起检查。

每次反射同时作用于矩阵与右端。右栏上方两行用于回代，最后一行留下无法消去的残差分量；其绝对值为 5/3。

运行下一个实验前，留意第二步的活动向量：它的第一项是零，但整个向量不是零。逐步执行两次反射，检查第一列的零是否保留，右端是否同步改变。也可以只变矩阵、故意漏变右端，比较回代结果在原问题中是否还满足正交条件。

对于 $m\times n$ 矩阵，一次反射作用于长度约 $m-k+1$ 的一个向量，需要一个内积与一次更新，约 $4(m-k+1)$ 次算术操作。第 $k$ 步处理约 列，合计主项为

右端跟着变换的成本约为 $4mn$，之后还有约 $n^2$ 的回代。这里估计的是操作次数，不是具体机器的运行时间。复用分解时，矩阵不变就不必再做那部分最昂贵的工作。

正规方程在哪里丢掉了信息

要说清“条件数平方”，先固定对象。对满列秩矩阵，记单位向量经 $A$ 作用后的最大、最小长度为

它们也是最大、最小奇异值。矩阵的 2-范数条件数为 $\kappa_2(A)=\sigma_{\max}/\sigma_{\min}$。

由于 $\|Az\|_2^2=z^TA^TAz$，对称正定矩阵 $A^TA$ 的最大、最小特征值分别是 。这里使用实对称矩阵的正交对角化：在单位向量上，二次型的极值就是最大、最小特征值。因此

这是上述条件下的精确等式，不是近似经验。而 $\|Q_1Rz\|_2=\|Rz\|_2$，说明 $A$ 与薄 QR 中的 有相同的这两个伸缩极值，故 。

一个能看清舍入位置的矩阵

取 $\varepsilon>0$，

精确解是 $(1,1)^T$，而

向量 $(1,1)^T$ 与 $(1,-1)^T$ 分别对应特征值 $2+\varepsilon^2$ 和 。后一方向表示两个系数一增一减：第一条观测只看到它们的和，难以察觉这种改变。

于是

在 binary64 中取 $\varepsilon=10^{-8}$，形成 Gram 矩阵时，$1+\varepsilon^2$ 会舍入到 $1$。计算所得矩阵变成两行相同的矩阵；原 $A_\varepsilon$ 中仍然存着的 ，已经在平方和加法中丢了。

精确 Gram 矩阵的对角线比非对角线大 10⁻¹⁶；用 binary64 形成它时，这份差别被舍去。原矩阵中的 10⁻⁸ 仍可表示，因此换一种计算路径具有实际意义。

Householder QR 直接处理 $A_\varepsilon$，可以避开这一次信息损失。规范实现具有良好的后向稳定性：通常可以把计算结果解释为附近矩阵、附近右端的最小二乘解。它并不保证附近问题与原问题有接近的系数，那仍取决于问题本身的敏感性。

算法稳定以后，数据仍然可能敏感

将刚才右端的第二项改成 $\varepsilon+\delta$。精确正规方程的右端增加 $(\varepsilon\delta,0)^T$。两行相加、相减可得

当 $\delta=\varepsilon$ 时，两个系数约为 $1.5$ 和 $0.5$。输入只改了一个很小的数，两个系数却各改了约 $0.5$。QR 即使把这个新问题算得十分准确，也不该“纠正”回原先的 $(1,1)$。

下面的实验把两件事分开显示：固定数据，比较正规方程与 Householder QR 的计算结果；再改变 $\delta$，比较新问题的精确参考解与旧问题的解。把 $\varepsilon$ 调小，找出 Gram 矩阵首次丢掉可见差异的位置，同时观察预测值是否像系数那样剧烈变化。

还有一条限制需要保留：最小二乘对矩阵扰动的敏感性不总由 $\kappa_2(A)$ 单独描述，最优残差也会参与。看一个只有一个系数的例子：

原问题的解为 $c=1$，$\kappa_2(A)=1$。若矩阵改为 $A_\eta=(1,\eta)^T$，则

$M$ 很大时，列方向的一点改变就能显著改变系数。原来垂直于列空间的那部分观测，经轻微旋转后进入了可拟合方向。所以不能把“QR 避免显式平方条件数”简化为“任何最小二乘的相对误差都不超过 $\kappa_2(A)u$”。

检查结果时，把问题分开

得到系数后，应回到原始数据计算预测和残差，再检查 $A^Tr$ 是否接近零。残差小不小，回答的是所选模型能否贴近观测；正交条件则用来核对求解是否满足最优性要求。数值上应结合列的尺度看这个量；要把它变成系数离最优解有多远的保证，还需要条件数等信息。这些检查不能互相替代。

即便正交条件满足，残差随 $t$ 呈现明显的弯曲规律，也可能说明直线模型漏掉了结构。增加一列已知函数，会扩大可拟合的空间，因此最优训练残差不会增加：旧系数加上一个零，就是新问题中的合法候选。但更小的训练残差本身不能保证新数据上的预测更好。

若矩阵列近相关，应检查数据扰动对系数的影响，而不仅仅是打印更多小数。若已经秩亏，前面的可逆 $R$ 回代不再适用；带列主元的 QR 或 SVD 可以进一步处理秩与解的选择，本章不把“除以一个极小对角元”当作替代方案。

改变观测的单位时也要留意目标。对 $t$ 做可逆的中心化、缩放，只要同时正确换回系数，所张成的模型空间可以不变。对观测 $b_i$ 取对数后再最小二乘，却是在最小化对数差的平方，已经换了误差度量，不能冒称仍解出了原始纵向误差问题。

留下几道需要动笔的问题

题 1：拿掉截距以后。 用 $y=ct$ 拟合 $(1,1),(2,2),(3,2)$。求 $c$ 与残差，检验 $t^Tr=0$。残差之和还是零吗？

题 2：两列完全相同。 对 $A=\begin{bmatrix}1&1\\2&2\end{bmatrix}$、$b=(1,0)^T$，给出全部最小二乘系数、唯一的预测和残差，并解释正规方程为什么仍然有效。

题 3：薄 QR 丢掉的那一部分。 已知

求最小二乘系数和残差。再将 $b$ 改成 $b+\delta(1,-1,0)^T$，说明系数与最小残差如何改变。

题 4：第一项为零的反射。 对 $z=(0,3,4)^T$，按本章规则求 $\alpha,w,H$，验证 $Hz$。如果把 $s(0)$ 错设成零，会发生什么？

题 5：保留反射，换一个右端。 沿用完整算例的 $3\times2$ 矩阵，把 $b$ 改为 $(1,1,1)^T$。复用两次反射求解，给出系数与最小残差长度。

题 6：拟合值与系数的变化不在一个尺度。 在 $A_\varepsilon$ 例子中，取 $\delta=\varepsilon$。求系数变化 $\Delta c$，并证明 $\|A_\varepsilon\Delta c\|_2\le|\delta|$。这与系数相差约 矛盾吗？

题 7：条件数为 1 也要看扰动对象。 对正文的 $A=(1,0)^T,b=(1,M)^T$，取 $M=1000$，把矩阵第二项改为 。求新系数。若只看 ，会漏掉什么？

题 8：从三角结果检查反射实现。 某程序报告完整 QR 的变换结果为

求系数与最小残差长度。若原输入的 $\|b\|_2=5$，这些结果能否同时来自一个正交 $Q$？指出检查依据。