解方程组时，把第二行减去第一行的若干倍，是很熟悉的动作。问题出在这个“若干倍”上。如果倍数达到一万，原来不显眼的舍入误差也会一起被放大；而换一下方程的排列顺序，倍数可能就变成万分之一。

本章把消元过程展开来算。我们既要得到解，也要留下足够的信息，判断哪里发生了舍入、换行是否正确，以及计算结果代回原方程后还差多少。以下主要讨论 $A\in\mathbb R^{n\times n}$ 非奇异的系统 $Ax=b$；奇异系统需要另外检查相容性，不能照着非零主元公式一路除下去。

从一条方程中消去一个未知数

若第 $k$ 步当前的主元 $a_{kk}^{(k)}\ne0$，对下面第 $i$ 行取乘子

$m_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}},\qquad i>k.$

把第 $i$ 行减去 $m_{ik}$ 倍第 $k$ 行，就能把该行第 $k$ 列消为零。其他系数与右端也要同步更新：

$a_{ij}^{(k+1)}=a_{ij}^{(k)}-m_{ik}a_{kj}^{(k)},\qquad b_i^{(k+1)}=b_i^{(k)}-m_{ik}b_k^{(k)}.$

上标表示当前消元阶段，不能一直从原始矩阵里取乘子。行变换可以用相反操作还原，因此保持解集；只改矩阵、不改右端则会改变方程本身。

来看一个会在第二步换行的系统：

$A=\begin{bmatrix}4&1&1\\2&1&5\\1&3&1\end{bmatrix},\qquad b=\begin{bmatrix}5\\-1\\6\end{bmatrix}.$

第一列中绝对值最大的是 4，用它作主元。第二、三行的乘子分别为 $1/2$ 和 $1/4$，所以增广矩阵变成

$\left[\begin{array}{ccc|c}4&1&1&5\\0&1/2&9/2&-7/2\\0&11/4&3/4&19/4\end{array}\right].$

第二列只比较尚未处理的第二、三行：$11/4$ 比 $1/2$ 大，把它换到第二行。这里第一行已经处理完，不再参与候选。交换后用乘子

$m_{32}=\frac{1/2}{11/4}=\frac2{11}$

消去第三行，得到

$\left[\begin{array}{ccc|c}4&1&1&5\\0&11/4&3/4&19/4\\0&0&48/11&-48/11\end{array}\right].$

末行只有一个未知数，给出 $x_3=-1$；第二行给出 $(11/4)x_2-3/4=19/4$，故 ；第一行给出 ，故 。代回原始的三条方程，右端分别是 ，全部吻合。

五阶矩阵消元的几个阶段：浅蓝色表示已经消去的零，橙色标出主元位置。图只展示结构变化，不是正文三阶算例的数值矩阵。

这样的从末行往前算称为回代。对于上三角系统 $Ux=z$，一般公式是

$x_i=\frac{z_i-\sum_{j=i+1}^n u_{ij}x_j}{u_{ii}},\qquad i=n,n-1,\ldots,1.$

求第 $i$ 个分量时，右侧需要的 $x_{i+1},\ldots,x_n$ 已经算好。下三角系统 $Lz=c$ 则从首行开始：

$z_i=\frac{c_i-\sum_{j=1}^{i-1}l_{ij}z_j}{l_{ii}},\qquad i=1,2,\ldots,n.$

空和取零。两种代入法都要求对角元非零；三角矩阵的行列式就是对角元之积，所以这个条件也恰好保证三角系统有唯一解。

把消元留下的乘子保存成 LU

刚才每次消元都用了一些乘子。如果下一次只改变右端 $b$，这些乘子完全可以继续使用，因为它们只由矩阵 $A$ 决定。把它们保存下来，就得到矩阵分解的做法。

暂时不换行。把矩阵分成首行、首列和剩余部分，设首元为 $\alpha\ne0$，则直接相乘可以检查

$\begin{bmatrix}\alpha&v^T\\c&D\end{bmatrix}=\begin{bmatrix}1&0\\c/\alpha&I\end{bmatrix}\begin{bmatrix}\alpha&v^T\\0&D-cv^T/\alpha\end{bmatrix}.$

右下角 $D-cv^T/\alpha$ 正是消去第一列后留下的部分。左边那个下三角因子记录了怎样把减掉的行加回去，所以其中存的是乘子本身，不是乘子的相反数。对剩余部分继续分解，便得到

$A=LU,$

其中 $L$ 是对角元全为 1 的下三角矩阵，$U$ 是上三角矩阵。这个分块等式也说明，消元和 LU 并非两种互不相关的技巧：一个着眼于行操作，另一个保存这些操作的结果。

有换行时，约定用 $P$ 左乘来表示本次执行的行排列，写成

$PA=LU.$

对主例，$P$ 交换第二、三行，最终因子是

$L=\begin{bmatrix}1&0&0\\1/4&1&0\\1/2&2/11&1\end{bmatrix},\qquad U=\begin{bmatrix}4&1&1\\0&11/4&3/4\\0&0&48/11\end{bmatrix}.$

注意 $L$ 第一列现在的次序是 $1/4,1/2$，与第一次消元时写下的行次序相反。第二步交换方程时，它们已经经历过的第一步操作也跟着交换；如果只换当前工作矩阵、不换 $L$ 中已经完成的列，$PA=LU$ 就不成立了。检查乘积第二行是 $(1,3,1)$，第三行是 ，正好对应 。

第二步交换当前第二、三行时，第一列已经存下的乘子也要交换；这里仅展示已经完成的旧列，未计算的部分不参与交换。

实际保存分解时，可以同时维护一个排列向量。每步选择候选行后，交换工作矩阵的相应两行、排列向量的两个位置，以及 $L$ 的这两行中已经完成的列。尚未计算的列不用交换。上三角部分与严格下三角乘子也可以共用一块存储，但不能因此把原始 $A$ 丢掉：后面计算残差还需要它。

下面的实验按步展示消元，保留当前主元、乘子和排列。运行主例到第二列前，预测哪一行应该被选中；执行换行后，同时查看工作矩阵与 $L$ 的旧列。关闭主元策略可以作比较，但程序遇到零主元时必须停下，而不是继续生成无穷大。

分解存在的条件，和求解时的顺序

可逆不等于不换行也能消元。例如

$A=\begin{bmatrix}0&1\\1&0\end{bmatrix}$

可逆，但第一主元为零，不能直接计算 $m_{21}$。交换两行便成为单位矩阵。对非奇异 $A$，部分主元法在精确算术中总能找到可用候选：已经处理的部分形成非零对角的上三角块；如果剩余块的第一列全为零，剩余块便奇异，整个当前矩阵也奇异。这与可逆行变换保持原矩阵非奇异相矛盾。

若要求完全不换行，则对非奇异 $A$，其所有左上角 $k\times k$ 子矩阵都非奇异，是存在单位下三角 LU 分解的充要条件。必要性可以从

$\det A_k=u_{11}u_{22}\cdots u_{kk}$

看出来。充分性也可沿消元递推：第一个子式非零保证第一主元；消去第一列后，剩余块的各个左上角子式，乘以第一主元，等于原矩阵相应子式，所以仍非零，后续便能继续。这里不是拿行列式逐个计算主元的数值算法，而是在说明条件。

固定这个单位对角约定时，非奇异矩阵的 LU 若存在便唯一。假设 $L_1U_1=L_2U_2$，则

$L_2^{-1}L_1=U_2U_1^{-1}.$

左侧单位下三角，右侧上三角，一个矩阵同时具有这两种结构，只能是单位矩阵；因此两组因子相同。带不同排列的分解不必相同，不能把这里的唯一性扩大到任意 $P$。

得到 $PA=LU$ 后，原方程变为 $LUx=Pb$。必须按以下依赖求解：

$Lz=Pb,\qquad Ux=z.$

对主例，$Pb=(5,6,-1)^T$。前代得

$z_1=5,\qquad z_2=6-\frac54=\frac{19}4,$

$z_3=-1-\frac52-\frac2{11}\frac{19}4=-\frac{48}{11}.$

它正是增广消元产生的右端。回代于是重现 $(1,2,-1)^T$。这里的 $P$ 只换方程顺序，不换未知数顺序；若连列也交换，则还要另外追踪未知数的排列。

现在保持 $A$ 不变，改用 $b^{(2)}=(2,6,4)^T$。不重新消元，直接排列为 $(2,4,6)^T$，前代得

$z^{(2)}=\left(2,\frac72,\frac{48}{11}\right)^T.$

回代得到 $x^{(2)}=(0,1,1)^T$。这一组解用原始矩阵相乘确实给出 $(2,6,4)^T$。右端变了， 都保留下来。

主例的两个右端使用同一组 P、L、U。每条路径都先排列右端，再前代与回代；中间向量 z 与最终解 x 各有自己的位置。

运算花在哪里

对一般稠密 $n\times n$ 矩阵，第 $k$ 步计算 $n-k$ 个乘子，并更新 $(n-k)^2$ 个剩余元素，每个元素需要一次乘法和一次减法。因此分解运算数的主项来自

$\sum_{k=1}^{n-1}\left[(n-k)+2(n-k)^2\right]=\frac23n^3+O(n^2).$

主元比较与行交换增加的工作不改变这个三次主项。前代、回代处理的是逐行的内积，每个系统各为 $O(n^2)$，两次合起来的主项约为 $2n^2$。

于是同一矩阵配 $m$ 个右端，复用分解的成本主项约为 $\frac23n^3+2mn^2$。如果每个右端都重新分解，三次项就重复了 $m$ 遍。矩阵规模加倍时，分解主项约变成八倍。这里 表示系统规模，三次描述运算量的增长速度。

为求一个 $Ax=b$ 而显式形成 $A^{-1}$，通常没有必要。得到逆矩阵本身就需要额外求解与存储，最后乘 $b$ 又加入舍入。解三角系统已经能完成所需动作；后面写 $A^{-1}r$ 是数学关系，并不要求程序真的形成逆矩阵。

下面在已经分解好的主例上更换右端，逐行执行前代与回代。试着遗漏右端的排列，比较得到的解与原方程残差，再恢复正确步骤。你会看到，因子看起来都对，并不足以保证使用方式也对。

部分主元控制乘子，不能消除所有风险

部分主元法在每一步都比较当前列尚未处理的候选，选绝对值最大者。 它不是只在主元看起来很小时才启动的应急措施。相等时可约定选行号较小者，保证过程可复现。

在精确算术中，这个选择直接给出 $|m_{ik}|\le1$。它控制当前消元所用的倍数，与矩阵本身是否对数据敏感是两回事。

用一个刻意降低精度的实验把差别放大。设

$A=\begin{bmatrix}10^{-4}&1\\1&1\end{bmatrix},\qquad b=\begin{bmatrix}1\\2\end{bmatrix}.$

规定每次加减乘除后都舍入为三位十进制有效数字，以下算例不遇到中点平局。若不换行，乘子是一万。更新第二行时，$1-10000=-9999$ 舍入为 $-10000$，右端 $2-10000=-9998$ 也舍入为 $-10000$。回代得到 $$，第一行便给出 。

若先交换两行，乘子是 $10^{-4}$。第二行更新后的系数 $0.9999$ 和右端 $0.9998$ 在三位模型中都成为 1，回代得到 $(1,1)^T$。两条路径都会舍入，但代回原始数据的残差差别很大：

真解为 $(10000/9999,9998/9999)^T$。用第 2 章的无穷范数，$\|A\|_\infty=2$、$\mathrm{\parallel }{A}^{-1}{\mathrm{\parallel }}_{\mathrm{\infty }}=2\mathrm{/}(1-{10}^{-}$，所以 。这是一个条件不坏的系统，却能被不合适的计算顺序算坏。

每步运算舍入到三位十进制有效数字时，两条消元路径得到不同近似解。图中的残差另用原始数据计算，用来比较结果的实际偏差。

乘子小也不意味着后续元素不会长大。考虑

$W=\begin{bmatrix}1&0&0&1\\-1&1&0&1\\-1&-1&1&1\\-1&-1&-1&1\end{bmatrix}.$

按相等取上行的部分主元规则，不发生换行。各个消元乘子都是 $-1$，因此操作是把主元行加到下面。第一步后，下面各行末列从 1 变成 2；第二步后，剩余两行末列变成 4；第三步后，最后一个元素变成 8。最终 $U$ 的末列是 $(1,2,4,8)^T$。

定义增长因子为消元各阶段出现的系数绝对值最大值，除以原矩阵系数绝对值最大值，包含原始阶段但不包含右端。这个例子中 $\rho=8$。一般情况下，若上一阶段所有相关元素绝对值至多为 $M$，则更新满足

$|a_{ij}-m_{ik}a_{kj}|\le M+|m_{ik}|M\le2M.$

经过最多 $n-1$ 次这样的更新，得到粗上界 $\rho\le2^{n-1}$。部分主元避免任意大的乘子，却仍允许最坏情况下的指数增长。

这些中间量较大时，对它们的相对舍入就可能对应原始数据尺度上较大的绝对误差。增长因子说明风险如何产生，不能据此断言每个右端都一定算错；上面的整数例子在通常浮点运算中甚至可以完全精确。最终仍要检查实际残差与条件性。

把残差放回数据的尺度

按照第 2 章的约定，残差为 $r=b-A\widehat x$，所以

$A(x-\widehat x)=r,\qquad x-\widehat x=A^{-1}r.$

若 $b\ne0$，由范数不等式与 $\|b\|\le\|A\|\|x\|$，可得

$\frac{\|x-\widehat x\|}{\|x\|}\le\kappa(A)\frac{\|r\|}{\|b\|}.$

要使用同一向量范数及其相容矩阵范数。这个界可以很保守，但已经指出：小相对残差不必对应小相对解误差。若 $b=0$，就不能除以 $\|b\|$，应改看绝对误差或明确另外的尺度。

常用的整体检查量还包括

$\eta_\infty=\frac{\|r\|_\infty}{\|A\|_\infty\|\widehat x\|_\infty+\|b\|_\infty}.$

把 $A,b$ 同时乘任意非零常数，分子、分母一起缩放，这个值不变。报告裸残差“只有 $10^{-8}$”却不说明数据尺度，没有同样的解释力。若分母为零，在这里非奇异 $A$ 的范围内必有 $b=\widehat x=0$，已是精确解，可把指标记为零。

如果不同方程的量级差异很大，还应逐行看。令

$s_i=\sum_j|a_{ij}|\,|\widehat x_j|+|b_i|,\qquad \omega=\max_i\frac{|r_i|}{s_i}.$

$s_i=0$ 时该行所有参与量都为零，$r_i=0$，相应比值按零处理。$\omega$ 具有具体的后向误差含义：它是能使 $\widehat x$ 精确满足 的最小 ，其中要求逐元素满足 、。

为什么正好是这个比值？由扰动方程，$r=\Delta A\widehat x-\Delta b$，故 $|r_i|\le\varepsilon s_i$，任何可行扰动都必须让 。反过来，当 时设 ，选择

$\Delta a_{ij}=\theta_i|a_{ij}|\operatorname{sign}(\widehat x_j),\qquad \Delta b_i=-\theta_i|b_i|.$

取 $\operatorname{sign}(0)=0$，则该行的 $\Delta A\widehat x-\Delta b$ 正好为 $\theta_i s_i=r_i$，且所有相对变化都不超过 。 的行取零扰动即可。这给出了达到下界的构造。

第二条方程虽然整体很小，残差却与自身右端一样大。整体尺度指标与逐行尺度指标回答不同的问题，不能直接当作解误差。

例如 $A=\operatorname{diag}(1,10^{-8})$、$b=(1,10^{-8})^T$，程序返回 。真解是 ，第二分量完全没算对。整体指标 看起来很小，但第二行的 、，所以 ，立即暴露了这条方程相对自身尺度的偏差。

若把第二条方程乘 $10^8$，系统变成 $Ix=(1,1)^T$，解没有改变；整体指标变为 $1/2$，逐分量指标仍为 1。这解释了为什么有时需要行尺度调整或逐行检查。调整尺度不会凭空增加原始测量数据的准确位数。

计算残差应使用原始的 $A,b$，不能只拿算出的 $L,U$ 自己检验自己。若数据来自测量，还要区分“准确解了已存储的系统”和“准确描述了真实问题”；后者取决于数据误差及条件性。

用旧因子改进解，把残差算得更认真些

如果残差算得可靠，方程 $Ad=r$ 的精确解正好是缺少的修正 $d=x-\widehat x$。因此更新应为 $\widehat x+d$，不是减去 。已有 时，求修正只需

$Lw=Pr,\qquad Ud=w,\qquad \widehat x_{\text{新}}=\widehat x+d.$

不再分解 $A$，每轮主要成本是一次矩阵向量乘法与两次三角求解，仍为 $O(n^2)$。

例：低精度残差把问题藏起来了。 取

$A=\begin{bmatrix}4&1\\1&3\end{bmatrix},\qquad b=\begin{bmatrix}1\\2\end{bmatrix},\qquad x=\left(\frac1{11},\frac7{11}\right)^T.$

三位十进制有效数字下，因子中的乘子 $1/4$ 和第二主元 $2.75$ 恰好能表示。前代得 $z=(1,1.75)^T$，回代时 $1.75/2.75$ 舍入为 $0.636$，再算 得 。

对这个近似值，使用更高精度计算原始残差，可以看见

$r=\begin{bmatrix}1-(4\times0.091+0.636)\\2-(0.091+3\times0.636)\end{bmatrix}=\begin{bmatrix}0\\0.001\end{bmatrix}.$

仍用三位模型求修正：$w=(0,0.001)^T$，$d_2$ 舍入为 $0.000364$，$d_1=-d_2\mathrm{/}4=-0.0000910$。把更新本身保留在较高精度中，得到

$\widehat x_{\text{新}}=(0.090909,0.636364)^T.$

再次按这些十进制数核对，残差为 $(0,-0.000001)^T$，最大绝对解误差也从约 $3.64\times10^{-4}$ 降到约 $3.64\times10^{-7}$。修正方程并未精确求解，仍然带来了明显改善。

若第一轮残差也只按三位计算，则 $3\times0.636=1.908$ 先成为 $1.91$，加上 $0.091$ 后的 $2.001$ 又成为 $2.00$，第二个残差便被算成零。第一行的残差同样是零，于是程序误以为已经没有可改进之处。精细计算残差有用，是因为那里正在做接近量的相减。

同一个近似解代入第二条方程，逐步三位舍入会把残差 0.001 算成零。右栏前两处箭头表示舍入，不能读成精确等号。

也可以用一个理想化模型看清收敛条件。把一次近似求解的作用记作固定矩阵 $B$，暂时不计每轮的新舍入误差，令 $e=x-\widehat x$。由于 $r=Ae$，更新 $Br$ 后有

$e_{\text{新}}=e-BAe=(I-BA)e.$

若某个相容范数下 $q=\|I-BA\|<1$，就有 $\|e_k\|\le q^k\|e_0\|$。这个条件充分但不是必要；真实浮点求解器也并非严格固定的线性映射 ，这里用它分离“修正方向是否有效”与“新误差有多大”。

若计算出的残差是 $Ae+\rho_k$，更新另有误差 $\zeta_k$，同一模型给出

$\|e_{k+1}\|\le q\|e_k\|+\|B\|\,\|\rho_k\|+\|\zeta_k\|.$

即使 $q<1$，后两项持续存在，也只能降到相应的误差水平。因子太差、问题过于敏感、残差不可靠或更新小到无法表示，都可能使改进停滞。实际程序要记录原始残差、修正大小与迭代次数；达到轮数限制或不再改善时说明原因，不能用“又迭代了一次”代替精度证据。

下面保留低精度因子，分别用低精度和浏览器的普通浮点精度计算残差。预测哪一种会在初值处错误停止，再逐轮比较解、残差和真实误差。实验中的参考解已知，实际任务未必有这样的参照，因此不能省去条件性与尺度检查。

让每一步都能独立核对

1. 用部分主元法分解并求解

$A=\begin{bmatrix}0&2&1\\2&1&0\\1&0&2\end{bmatrix},\qquad b=(0,1,5)^T.$

写出 $P,L,U$ 与前代结果。再保持矩阵不变，求右端 $(3,3,3)^T$ 的解。

2. 在正文的三阶主例中，程序交换第二、三行后忘记交换 $L$ 已经存下的第一列，留下 $l_{21}=1/2,l_{31}=1/4$。只检查乘积第二行，说明它为何不满足 $PA=LU$。

3. 对 $A(t)=\begin{bmatrix}t&1\\1&1\end{bmatrix}$，哪些 $t$ 使矩阵可逆？其中哪些值允许不换行的单位下三角 LU 分解？写出因子，并解释 的情况。

4. 三位模型的小主元例子中，不换行得到 $(0,1)^T$。具体指出丢失信息的两次舍入。若把最终输出显示成十位小数，能否补回它们？

5. 把增长例子中严格下三角的 $-1$ 改为 $-t$，其中 $0\le t\le1$，其余结构不变并推广到 $n$ 阶。相等时仍选上行。证明最后一个主元为 $(1+t)^{n-1}$，并说明为什么 并不保证增长因子为 1。

6. 同一个 $300\times300$ 稠密矩阵要处理 20 个右端。只按本章主项估算，每次重新分解与复用一次分解分别需要多少次算术运算？为什么不能把这个比值直接当成实际加速倍数？

7. 对 $A=\operatorname{diag}(1,10^{-8})$、$b=(1,10^{-8})^T$，这次返回 。求无穷范数下相对解误差、相对右端残差、 与 ，比较各自解释的对象。

8. 在理想化改进模型中，取 $A=\operatorname{diag}(2,5)$、$b=(2,5)^T$，初值为零。若 $B=\operatorname{diag}(0.4,0.1)$，计算前两次更新并给出收缩因子。再把 的第一个对角元改为 1.2，说明会发生什么。