数值分析 04：线性方程组的直接解法

上一篇的方程组 Newton 法在每一步都要解一个线性方程组。本篇集中研究这个线性子问题：给定 $A\in\mathbb R^{n\times n}$ 和 $\mathbf b\in\mathbb R^n$，计算

$$A\mathbf x=\mathbf b.$$

在精确算术中，消元法经过有限次运算便可得到解，因此称为直接法。浮点计算却带来三个必须分开的疑问：矩阵是否可逆，数据扰动会使精确解变化多少，算法产生的舍入误差又有多大。可逆性、问题的条件性和算法的稳定性不是同一件事。

除特别说明外，本篇假设 $A$ 为实方阵。讨论相对误差时还假设相应的分母不为零；所有运算量均按加、减、乘、除各计一次浮点运算，只保留关于 $n$ 的主要项。

---

1. 向量范数与矩阵范数

1.1 用范数统一衡量误差

定义（向量范数）

映射 $\|\cdot\|:\mathbb R^n\to[0,\infty)$ 若满足下列条件，就称为向量范数：

1. $\|\mathbf x\|\ge0$，且 $\|\mathbf x\|=0$ 当且仅当 $\mathbf x=\mathbf0$；
2. 对任意实数 $c$，有 $\|c\mathbf x\|=|c|\,\|\mathbf x\|$；
3. $\|\mathbf x+\mathbf y\|\le\|\mathbf x\|+\|\mathbf y\|$。

最常用的三个向量范数是

$$\|\mathbf x\|_1=\sum_{i=1}^n|x_i|, \qquad \|\mathbf x\|_2=\left(\sum_{i=1}^n|x_i|^2\right)^{1/2}, \qquad \|\mathbf x\|_\infty=\max_i|x_i|.$$

有限维空间中的所有范数都给出相同的收敛概念，但数值大小和误差界中的常数会随范数改变。写条件数或相对误差时，应说明所用范数，不能把不同范数下的量直接代入同一个界。

1.2 从向量范数诱导矩阵范数

定义（诱导矩阵范数）

给定一个向量范数，与之对应的诱导矩阵范数为

$$\|A\|=\max_{\mathbf x\ne\mathbf0} \frac{\|A\mathbf x\|}{\|\mathbf x\|} =\max_{\|\mathbf x\|=1}\|A\mathbf x\|.$$

它立即满足相容性与次乘性：

$$\|A\mathbf x\|\le\|A\|\,\|\mathbf x\|, \qquad \|AB\|\le\|A\|\,\|B\|.$$

常用诱导范数可以写成

$$\|A\|_1=\max_j\sum_{i=1}^n|a_{ij}|, \qquad \|A\|_\infty=\max_i\sum_{j=1}^n|a_{ij}|,$$

以及

$$\|A\|_2=\sqrt{\lambda_{\max}(A^{\mathsf T}A)} =\sigma_{\max}(A).$$

因此 $1$ 范数是最大列绝对值和，$\infty$ 范数是最大行绝对值和，$2$ 范数是最大奇异值。

另一个常用量是 Frobenius 范数

$$\|A\|_{\mathrm F} =\left(\sum_{i,j}|a_{ij}|^2\right)^{1/2}.$$

它满足次乘性，也满足 $\|A\mathbf x\|_2\le\|A\|_{\mathrm F}\|\mathbf x\|_2$，但它不是由向量 $2$ 范数诱导的矩阵范数。例如 $\|I\|_{\mathrm F}=\sqrt n$，而任意诱导范数都有 $\|I\|=1$。

---

2. 三角方程组与 Gauss 消元

2.1 前代与回代

若 $L$ 为下三角矩阵且对角元非零，方程 $L\mathbf y=\mathbf b$ 可依次求解：

$$y_i=\frac{1}{l_{ii}} \left(b_i-\sum_{j=1}^{i-1}l_{ij}y_j\right), \qquad i=1,2,\ldots,n.$$

这称为 前代。若 $U$ 为对角元非零的上三角矩阵，方程 $U\mathbf x=\mathbf y$ 用 回代：

$$x_i=\frac{1}{u_{ii}} \left(y_i-\sum_{j=i+1}^{n}u_{ij}x_j\right), \qquad i=n,n-1,\ldots,1.$$

每次三角求解大约需要 $n^2$ 次运算，存储矩阵本身需要 $O(n^2)$ 个数。若三角矩阵有额外的带状或稀疏结构，实际成本应按非零元计算。

2.2 消去一个未知量

Gauss 消元的目标很具体：从左到右逐列工作，把主对角线下方的系数变成零，最终得到上三角方程组。为了让系数和右端始终接受同样的行变换，把方程组写成增广矩阵

$$[A\mid\mathbf b] = \left[ \begin{array}{cccc|c} a_{11}&a_{12}&\cdots&a_{1n}&b_1\\ a_{21}&a_{22}&\cdots&a_{2n}&b_2\\ \vdots&\vdots&&\vdots&\vdots\\ a_{n1}&a_{n2}&\cdots&a_{nn}&b_n \end{array} \right].$$

一次消元只做一件事。若要用第 $k$ 行消掉第 $i$ 行的第 $k$ 个系数，其中 $i>k$ 且主元 $a_{kk}^{(k)}\ne0$，先取乘子

$$m_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}.$$

再对整行执行

$$R_i\leftarrow R_i-m_{ik}R_k.$$

这里 $R_k$ 是提供消元基准的主元行，$R_i$ 是正在修改的目标行。乘子这样选择，是因为目标位置会恰好变为

$$a_{ik}^{(k)}-m_{ik}a_{kk}^{(k)} =a_{ik}^{(k)} -\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}a_{kk}^{(k)} =0.$$

行中的其余系数和右端项也必须减去同一倍数，不能只把目标系数直接写成零。这个行变换是可逆的：从新行加回 $m_{ik}R_k$ 就恢复原行，所以它不改变方程组的解集。

处理完第 $k$ 列主元下方的所有行后，第 $k$ 个未知量便不再出现在后面的方程中。随后转向右下角尚未处理的子矩阵，重复“选主元—消去下方元素”，直到矩阵成为上三角形。

2.3 一个完整的三阶消元过程

考虑方程组

$$\begin{aligned} 2x+y-z&=8,\\ -3x-y+2z&=-11,\\ -2x+y+2z&=-3. \end{aligned}$$

其初始增广矩阵为

$$\left[ \begin{array}{rrr|r} 2&1&-1&8\\ -3&-1&2&-11\\ -2&1&2&-3 \end{array} \right].$$

第一列中绝对值最大的系数是 $-3$，所以先交换 $R_1$ 与 $R_2$。此时用主元 $-3$ 消去下面的 $2$ 与 $-2$：

$$m_{21}=\frac{2}{-3}=-\frac23, \qquad m_{31}=\frac{-2}{-3}=\frac23.$$

例如第一项行变换为

$$R_2\leftarrow R_2-\left(-\frac23\right)R_1 =R_2+\frac23R_1,$$

目标系数确实成为 $2+\frac23(-3)=0$。完成第一列后得到

$$\left[ \begin{array}{rrr|r} -3&-1&2&-11\\ 0&\frac13&\frac13&\frac23\\ 0&\frac53&\frac23&\frac{13}{3} \end{array} \right].$$

接下来只看右下角 $2\times2$ 子矩阵。第二列的两个候选主元是 $1/3$ 和 $5/3$，所以交换 $R_2$ 与 $R_3$，再用

$$m_{32}=\frac{1/3}{5/3}=\frac15$$

执行 $R_3\leftarrow R_3-\frac15R_2$。上三角增广矩阵为

$$\left[ \begin{array}{rrr|r} -3&-1&2&-11\\ 0&\frac53&\frac23&\frac{13}{3}\\ 0&0&\frac15&-\frac15 \end{array} \right].$$

最后从最下面的方程向上回代：

$$z=-1,\qquad y=\frac{13/3-(2/3)z}{5/3}=3, \qquad x=\frac{-11+y-2z}{-3}=2.$$

下方实验把刚才的过程拆成八个动作。每次只改变一行或求出一个未知量，并同时说明“正在做什么”“为什么这一项会变成零”。建议先逐步前进，再用“上一步”回看矩阵变化。

2.4 运算量与原位存储

对稠密 $n\times n$ 矩阵，消元或等价的 LU 分解主要需要

$$\frac{2}{3}n^3+O(n^2)$$

次运算；一次回代需要 $n^2+O(n)$ 次运算。因此求解一个一般稠密系统的主要成本来自分解，而不是三角求解。

实现中通常用上三角部分存放 $U$，用严格下三角部分存放乘子，即 $L$ 的非对角元素。右端项可以与系数矩阵一起更新，也可以在分解完成后单独前代。后一种组织方式更适合多个右端项。

2.5 零主元不等于矩阵奇异

不换行的消元可能遇到零主元，即使原矩阵可逆。例如

$$A=\begin{pmatrix}0&1\\1&1\end{pmatrix}, \qquad \det A=-1.$$

第一步主元为零，公式无法继续；交换两行后便可正常消元。因此“消元失败”有时只说明当前行次序不合适。

若当前列从第 $k$ 行向下的候选元素全为零，则换行也找不到非零主元。对方阵的标准消元过程而言，这表示矩阵秩不足，系统不存在唯一解。

---

3. 列主元消元与置换矩阵

3.1 在当前列中选择主元

列主元策略也称部分选主元。在第 $k$ 步，从尚未处理的行中选择

$$p\in\operatorname*{arg\,max}_{k\le i\le n}|a_{ik}^{(k)}|,$$

交换第 $k$ 行与第 $p$ 行，再执行消元。若最大值为零，则矩阵奇异；否则每个本步乘子都满足

$$|m_{ik}|\le1.$$

这一策略避免除以同列中不必要的小数，并控制消元乘子的大小。它的比较成本为 $O(n^2)$，不改变分解的 $O(n^3)$ 主量级。

行交换可以用置换矩阵表示。把所有交换的累计结果记为 $P$，保留所有乘子后得到

$$\boxed{PA=LU},$$

其中 $L$ 是单位下三角矩阵，$U$ 是上三角矩阵。求解 $A\mathbf x=\mathbf b$ 时相应地计算

$$L\mathbf y=P\mathbf b, \qquad U\mathbf x=\mathbf y.$$

不能在右端漏掉同样的行置换。在原位分解中，如果后续又交换行，先前已经存入 $L$ 的乘子也要按正确的置换关系处理。

3.2 选主元解决什么问题

部分选主元不能把奇异矩阵变成可逆矩阵，也不能改善原问题的条件数。它控制的是计算过程，主要作用是避免因小主元产生巨大乘子和严重的中间误差放大。

为了单独观察这个作用，考虑

$$A_\varepsilon= \begin{pmatrix}\varepsilon&1\\1&1\end{pmatrix}, \qquad \mathbf b_\varepsilon= \begin{pmatrix}1+\varepsilon\\2\end{pmatrix}.$$

它的精确解为 $(1,1)^{\mathsf T}$。若 $\varepsilon=10^{-8}$，不换行时首个乘子为 $1/\varepsilon=10^8$；列主元策略则先交换两行，首个乘子变为 $\varepsilon=10^{-8}$。两种算法处理的是同一个系统，区别只在消元前是否换行。

下方实验把两种结果并排放置。输入读入以及每次乘、除、加减后都舍入到所选十进制有效位；这是教学模型，并不等同于浏览器实际的二进制浮点格式。最终残差和解误差使用原始数据与普通双精度重新计算。

“小首主元”预设直接展示巨大乘子造成的信息损失；“零首主元”预设说明不换行失败不等于矩阵奇异；“增长因子”预设则说明部分选主元也不能排除所有中间元素增长。三个预设回答的是不同问题，不应混成一个结论。

选主元之后仍需报告奇异或接近奇异情形。把某个固定绝对阈值机械地用于所有矩阵并不可靠；实际判断通常结合矩阵尺度、工作精度和条件估计。

---

4. LU 分解与多个右端项

4.1 不换行的 LU 分解

定义（LU 分解）

若

$$A=LU,$$

其中 $L$ 为单位下三角矩阵，$U$ 为上三角矩阵，就称之为 $A$ 的 LU 分解。Gauss 消元中的乘子正好构成 $L$ 的严格下三角部分。

定理（无主元 LU 分解的存在与唯一性）

设 $A\in\mathbb R^{n\times n}$。存在单位下三角 $L$ 与对角元非零的上三角 $U$ 使 $A=LU$，当且仅当 $A$ 的所有顺序主子式都非零：

$$\det A_{1:k,\,1:k}\ne0, \qquad k=1,2,\ldots,n.$$

在这些条件下，分解是唯一的。只知道 $A$ 可逆并不足以保证不换行的 LU 分解存在；上一节的零首主元矩阵就是反例。

实际求解一般矩阵时通常采用带主元的 $PA=LU$。若矩阵本身具有能保证无零主元的结构，才可以有依据地省去选主元。

4.2 分解一次，求解多次

若有多个右端项

$$AX=B, \qquad B=(\mathbf b_1,\ldots,\mathbf b_s),$$

可以先用约 $2n^3/3$ 次运算求出 $PA=LU$，再对每一列依次解

$$L\mathbf y_j=P\mathbf b_j, \qquad U\mathbf x_j=\mathbf y_j.$$

每增加一个右端项只需两次三角求解，即 $O(n^2)$ 成本。若 $A$ 不变而右端项不断变化，应保存分解结果，而不是对每个右端项重新消元。

LU 分解也可用来计算行列式：

$$\det A=\det(P)\prod_{i=1}^n u_{ii}.$$

其中 $\det(P)=(-1)^s$，$s$ 是行交换次数。虽然也能通过解 $n$ 个右端项形成 $A^{-1}$，但求解 $A\mathbf x=\mathbf b$ 时通常没有必要显式计算逆矩阵；直接使用分解更省计算和存储，也避免额外舍入误差。

---

5. 利用矩阵结构

5.1 对称正定矩阵与 Cholesky 分解

定义（对称正定矩阵）

若 $A=A^{\mathsf T}$，并且对任意 $\mathbf x\ne\mathbf0$ 都有

$$\mathbf x^{\mathsf T}A\mathbf x>0,$$

则称 $A$ 为对称正定矩阵。

定理（Cholesky 分解）

实矩阵 $A$ 对称正定，当且仅当存在对角元为正的下三角矩阵 $L$，使

$$\boxed{A=LL^{\mathsf T}}.$$

要求 $L$ 的对角元为正时，分解唯一。元素递推公式为

$$l_{jj}=\left(a_{jj}-\sum_{k=1}^{j-1}l_{jk}^2\right)^{1/2},$$

以及对 $i>j$，

$$l_{ij}=\frac{a_{ij}-\sum_{k=1}^{j-1}l_{ik}l_{jk}}{l_{jj}}.$$

精确算术下，每个根号内的量都为正。计算中若出现非正量，可能说明输入不满足正定性，也可能说明矩阵接近正定边界并受到舍入影响；不能继续对负数开实平方根。

Cholesky 分解约需 $n^3/3$ 次运算，约为一般 LU 分解的一半，并且只需存储一个三角因子。求解时先前代 $L\mathbf y=\mathbf b$，再回代 $L^{\mathsf T}\mathbf x=\mathbf y$。

正定性不能用“所有元素都为正”代替。例如

$$\begin{pmatrix}1&2\\2&1\end{pmatrix}$$

的元素全为正，却有特征值 $3$ 和 $-1$，因此不是正定矩阵。对称性同样是 Cholesky 定理的必要假设。

5.2 三对角方程组与追赶法

考虑三对角矩阵

$$A= \begin{pmatrix} d_1&u_1&&\\ \ell_2&d_2&\ddots&\\ &\ddots&\ddots&u_{n-1}\\ &&\ell_n&d_n \end{pmatrix}.$$

追赶法就是专门针对三对角结构的无主元消元。令

$$\widetilde d_1=d_1, \qquad \widetilde b_1=b_1.$$

对 $i=2,\ldots,n$ 计算

$$m_i=\frac{\ell_i}{\widetilde d_{i-1}}, \qquad \widetilde d_i=d_i-m_i u_{i-1}, \qquad \widetilde b_i=b_i-m_i\widetilde b_{i-1}.$$

随后回代

$$x_n=\frac{\widetilde b_n}{\widetilde d_n}, \qquad x_i=\frac{\widetilde b_i-u_i x_{i+1}}{\widetilde d_i}, \quad i=n-1,\ldots,1.$$

算法只需要 $O(n)$ 次运算和 $O(n)$ 存储，但前提是所有修改后的主元 $\widetilde d_i$ 都非零。严格行对角占优的三对角矩阵，以及对称正定三对角矩阵，都能给出不发生零主元的充分保证。一般三对角矩阵仍可能需要主元交换；换行还会改变简单的带状结构和存储方式。

结构化算法的优势来自利用已知结构。把三对角矩阵先扩展成稠密矩阵再调用一般消元，会把线性成本变回立方成本，并浪费大量存储。

---

6. 条件数与数据扰动

6.1 线性系统的条件数

给定一个诱导矩阵范数，对可逆矩阵定义

$$\boxed{\kappa(A)=\|A\|\,\|A^{-1}\|}.$$

若 $A$ 奇异，约定 $\kappa(A)=\infty$。常数倍缩放不改变条件数，且

$$\kappa(A)\ge1, \qquad \kappa(cA)=\kappa(A)\quad(c\ne0).$$

在 $2$ 范数下，

$$\kappa_2(A)=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}.$$

最小奇异值接近零意味着矩阵接近奇异，条件数很大。条件数依赖于范数，也是给定扰动模型下的最坏方向放大上界；一次具体扰动不一定达到这个上界。

6.2 只扰动右端项

设

$$A\mathbf x=\mathbf b, \qquad A(\mathbf x+\Delta\mathbf x)=\mathbf b+\Delta\mathbf b.$$

由 $A\Delta\mathbf x=\Delta\mathbf b$ 可得

$$\boxed{ \frac{\|\Delta\mathbf x\|}{\|\mathbf x\|} \le\kappa(A) \frac{\|\Delta\mathbf b\|}{\|\mathbf b\|} }.$$

这个上界在最坏扰动方向上可以达到。若 $A$ 的条件数很大，即使算法完全准确，输入中很小的相对扰动仍可能导致明显的解变化；这属于问题本身的敏感性。

6.3 系数矩阵和右端项同时扰动

设

$$(A+\Delta A)(\mathbf x+\Delta\mathbf x) =\mathbf b+\Delta\mathbf b.$$

记

$$\varepsilon_A=\frac{\|\Delta A\|}{\|A\|}, \qquad \varepsilon_b=\frac{\|\Delta\mathbf b\|}{\|\mathbf b\|}.$$

若 $\kappa(A)\varepsilon_A<1$，则 $A+\Delta A$ 仍可逆，并有有限扰动界

$$\frac{\|\Delta\mathbf x\|}{\|\mathbf x\|} \le \frac{\kappa(A)}{1-\kappa(A)\varepsilon_A} (\varepsilon_A+\varepsilon_b).$$

分母条件不可省略。当 $\kappa(A)\varepsilon_A$ 接近 $1$ 时，线性化的一阶近似已经不足以可靠描述解的变化，扰动后的矩阵也可能接近奇异。

6.4 实验：近乎平行的两条直线

考虑

$$A_\varepsilon= \begin{pmatrix}1&1\\1&1+\varepsilon\end{pmatrix}, \qquad \mathbf b_\varepsilon= \begin{pmatrix}2\\2+\varepsilon\end{pmatrix}.$$

只要 $\varepsilon\ne0$，原系统的精确解始终是 $(1,1)^{\mathsf T}$。现在只给第二个右端分量增加 $\delta=t\varepsilon$。扰动后精确解为

$$\mathbf x_\delta= \begin{pmatrix}1-t\\1+t\end{pmatrix}.$$

因此在 $2$ 范数下，解的相对变化恰为 $|t|$。固定 $t$ 并减小 $\varepsilon$ 时，输入扰动 $|t\varepsilon|$ 越来越小，解的变化却保持不变。

图中两条原始直线随着 $\varepsilon$ 减小而趋于平行。扰动后的第二条直线与第一条直线交于新解；即使从整幅图上几乎看不出右端项的移动，交点也可能明显改变。

---

7. 残差、前向误差与后向误差

7.1 残差恒等式

设 $\widehat{\mathbf x}$ 是一个计算结果，定义残差和前向误差

$$\mathbf r=\mathbf b-A\widehat{\mathbf x}, \qquad \mathbf e=\mathbf x-\widehat{\mathbf x}.$$

由于 $A\mathbf x=\mathbf b$，精确算术下有

$$A\mathbf e=\mathbf r, \qquad \mathbf e=A^{-1}\mathbf r.$$

对同一诱导范数，可以推出双边关系

$$\frac{1}{\kappa(A)} \frac{\|\mathbf r\|}{\|\mathbf b\|} \le \frac{\|\mathbf e\|}{\|\mathbf x\|} \le \kappa(A) \frac{\|\mathbf r\|}{\|\mathbf b\|}.$$

所以小的相对残差只有在条件数不过大时，才能保证小的相对前向误差。实验中的扰动解若拿回原系统检查，其残差正好是 $-\Delta\mathbf b$；原始残差很小，而相对于原始解的误差可能被病态性放大。

7.2 把残差解释为后向误差

同一个计算结果精确满足

$$A\widehat{\mathbf x}=\mathbf b-\mathbf r.$$

因此，只允许右端项改变时，$\|\mathbf r\|/\|\mathbf b\|$ 是一个直接的范数相对后向误差：计算结果是邻近右端项 $\mathbf b-\mathbf r$ 所对应的精确解。

若允许 $A$ 与 $\mathbf b$ 同时逐项扰动，常用的范数后向误差还会按 $\|A\|\,\|\widehat{\mathbf x}\|+\|\mathbf b\|$ 缩放。分量后向误差则逐项考虑数据尺度；当矩阵元素相差悬殊或含有结构性零元时，它可能比单一范数指标更有解释力。

计算残差本身也会舍入。若 $A\widehat{\mathbf x}$ 与 $\mathbf b$ 很接近，两者相减可能丢失有效数字；需要高可信残差时，可以用更高精度累加内积，或使用带误差控制的专门实现。

---

8. 舍入误差与主元增长

8.1 后向稳定性的含义

一个线性求解算法若产生 $\widehat{\mathbf x}$，并且可解释为邻近问题

$$(A+\Delta A)\widehat{\mathbf x}=\mathbf b+\Delta\mathbf b$$

的精确解，其中相对扰动与单位舍入误差同量级，就称其在相应扰动模型下后向稳定。后向稳定说明算法没有把输入改变很多；最终前向误差还要乘上问题的条件数才能估计。

对采用部分选主元的 Gauss 消元，经典误差分析给出形如

$$\frac{\|\Delta A\|}{\|A\|} \lesssim p(n)\,\rho,u$$

的界。这里 $u$ 是单位舍入误差，$p(n)$ 是随维数适度增长的因子，$\rho$ 是主元增长因子；符号 $\lesssim$ 省略了高阶项和依赖具体范数的常数。

8.2 主元增长因子

把消元过程中出现的所有系数记为 $a_{ij}^{(k)}$，定义

$$\rho= \frac{\max_{i,j,k}|a_{ij}^{(k)}|} {\max_{i,j}|a_{ij}|}.$$

增长因子衡量中间元素相对于输入最大元素放大了多少。即使最终结果规模适中，过大的中间量也会使绝对舍入误差增大。

部分选主元保证乘子不超过 $1$，但不能保证增长因子始终很小。对 $n$ 阶矩阵有最坏情形上界

$$\rho\le2^{n-1},$$

并且存在达到该数量级的矩阵。实际问题中增长通常温和，这也是列主元消元广泛使用的原因；“通常稳定”仍不是对每个输入都无条件成立的定理。

矩阵缩放不良也会影响绝对阈值、溢出风险和误差估计。按行、列平衡尺度有时能改善计算表现，但它不会神奇地消除所有方向上的病态性，条件数必须针对缩放后的变量解释。

---

9. 迭代改进

9.1 基本步骤

已有分解和初始计算解 $\widehat{\mathbf x}^{(0)}$ 后，迭代改进重复：

$$\mathbf r^{(k)}=\mathbf b-A\widehat{\mathbf x}^{(k)},$$

利用原来的分解近似求解

$$A\mathbf d^{(k)}=\mathbf r^{(k)},$$

再更新

$$\widehat{\mathbf x}^{(k+1)} =\widehat{\mathbf x}^{(k)}+\mathbf d^{(k)}.$$

不必重新分解 $A$，所以每轮主要增加一次残差计算和两次三角求解，成本为 $O(n^2)$。

9.2 它能够改进什么

若残差和校正都在精确算术中完成，一步就会得到精确解；浮点实现只能近似这一理想过程。残差使用比因子分解和三角求解更高的精度时，通常更能保留被相消隐藏的信息，这构成混合精度迭代改进的基本思路。

迭代改进不是无条件补救。其收敛和可达到的精度取决于矩阵条件数、因子分解的后向误差、校正求解精度、残差计算精度以及更新时的舍入。若问题严重病态到输入精度已经不能确定解，迭代改进也无法恢复输入中不存在的信息。

停止时应检查校正量和重新计算的残差，并设置迭代上限。残差不再下降或校正量无法改变当前浮点数，说明继续使用同一精度通常没有收益。

---

10. 方法选择与常见误区

10.1 直接解法对照

矩阵类型	推荐分解或算法	主要分解成本	核心条件
一般稠密方阵	列主元 Gauss 消元，即 $PA=LU$	约 $2n^3/3$	主元搜索后不存在零主元
多个右端项、同一一般矩阵	保存 $PA=LU$，重复三角求解	分解一次；每个右端项 $O(n^2)$	与一般 LU 相同
实对称正定矩阵	Cholesky 分解 $A=LL^{\mathsf T}$	约 $n^3/3$	对称且正定
满足无零修改主元条件的三对角矩阵	追赶法	$O(n)$	各步修改主元非零

表中的“直接”表示精确算术下经过有限步骤完成，不表示浮点结果必然精确，也不表示计算时间与维数无关。

10.2 常见误区

1. “矩阵可逆，所以无需换行的消元一定成功。” 可逆矩阵仍可能出现零的中间主元；顺序主子式条件比可逆性更强。
2. “选主元会改善矩阵条件数。” 行交换只改变计算次序，不改变原线性映射对数据扰动的固有敏感性。
3. “残差接近零，所以解一定准确。” 前向误差界还包含条件数；病态系统可以同时具有小残差和明显解误差。
4. “对称矩阵都能做 Cholesky 分解。” 还必须正定；不定矩阵需要其他分解策略。
5. “三对角矩阵都能安全使用不选主元的追赶法。” 仍需保证修改主元非零，并考虑舍入稳定性。
6. “求解线性系统应先算逆矩阵。” 实际计算通常直接分解并做三角求解；显式逆矩阵不是必要中间量。
7. “部分选主元保证绝对稳定。” 它通常表现可靠，但误差界仍受主元增长、维数和工作精度影响。

---

11. 本篇要点

• 诱导矩阵范数把向量误差、矩阵作用和条件数放进同一套语言中；范数选择必须前后一致。
• Gauss 消元把一般系统化为三角系统，列主元策略通过换行控制消元乘子，并形成 $PA=LU$。
• 分解成本通常高于三角求解，因此多个右端项应复用同一分解；矩阵结构则可能把成本从立方量级降到平方或线性量级。
• Cholesky 分解只适用于对称正定矩阵，追赶法需要无零修改主元等条件。
• 条件数描述问题对数据扰动的敏感性，后向误差描述计算结果对应的邻近问题；小残差只有结合条件数才足以判断解误差。
• 部分选主元的舍入误差与主元增长有关；迭代改进可以复用已有分解，但不能恢复输入精度中已经丢失的信息。

下一篇将研究线性方程组的迭代解法，比较 Jacobi、Gauss–Seidel、SOR 与共轭梯度法，并用谱半径和矩阵范数判断迭代是否收敛。

参考资料

1. Richard L. Burden、J. Douglas Faires、Annette M. Burden：《Numerical Analysis》，第 10 版，Cengage Learning，线性方程组的直接解法相关章节。
2. Tobin A. Driscoll、Richard J. Braun：《Fundamentals of Numerical Computation》，线性方程组、矩阵范数与条件数相关章节。
3. Gene H. Golub、Charles F. Van Loan：《Matrix Computations》，第 4 版，Johns Hopkins University Press，2013 年。
4. Nicholas J. Higham：《Accuracy and Stability of Numerical Algorithms》，第 2 版，Society for Industrial and Applied Mathematics，2002 年。
5. E. Anderson 等：《LAPACK Users’ Guide》，第 3 版，Society for Industrial and Applied Mathematics，1999 年。