工程工具箱：误差诊断 · 算法选型 · 停止准则 · 函数地图

正文回答「为什么」，这一页回答「动手时先看哪一页、用哪个函数」。它是一张索引：拿到一个具体症状，先在本页定位到该算哪个量、该调哪个函数，再回正文读推导。本页不重述正文内容，也不重复复习层的公式卡与对照卡。

一、误差诊断：四个量、一张判据表

核心判据（把残差换算成误差）： $\dfrac{\lVert\boldsymbol{e}\rVert}{\lVert\boldsymbol{x}\rVert}\lesssim\mathrm{cond}(A)\cdot\dfrac{\lVert\boldsymbol{r}\rVert}{\lVert\boldsymbol{b}\rVert}$（定理 2.4.4）。求解器向后稳定只保证残差小；要谈解准不准，必须乘上条件数。丢掉几位有效数字用 $k\approx\log_{10}\mathrm{cond}(A)$ 估（2.4.2 换算表）。

「后向误差与条件数无关」的实测证据（同一批 Hilbert 矩阵，np.linalg.solve 求解、真解全为 $1$）：

读法：两行的后向误差都停在 $u\approx1.1\times10^{-16}$ 甚至更小——算法一直是向后稳定的；而前向误差从 $10^{-7}$ 涨到 $10^{-1}$，差了 6 个数量级。病态吃掉的是前向误差，与算法的稳定性无关（1.2.6 的 conditioning / stability 之分在这里落地）。

症状 → 先算哪个量 → 判据 → 处方：

可复算的两个锚点（本机 numpy 2.4.6，命令可直接跑）：8 阶 Hilbert 矩阵 $\mathrm{cond}_2=1.53\times10^{10}$，solve 最大误差 $1.2\times10^{-7}$；12 阶 $\mathrm{cond}_2=1.66\times10^{16}$，solve 最大误差 $0.52$（真值全为 $1$）——第二行说明解已经失去意义，而求解器并不报错。

二、算法选择树：解方程组与特征值

解 $A\boldsymbol{x}=\boldsymbol{b}$（方阵）

判定顺序（照这个顺序问四个问题）： ① 矩阵对称吗？正定吗？→ 决定能否用 Cholesky / LDL。② 稀疏吗（判据见第四节的经验法则）？→ 决定稠密库还是 scipy.sparse。③ 条件数量级多大？→ 决定要不要换路线（QR / SVD / 正则化 / 提高精度）。④ 同一个 $A$ 要解多少个右端？→ 多个右端先分解一次（lu_factor / cho_factor），再反复 *_solve。

超定最小二乘 $A\boldsymbol{x}\approx\boldsymbol{b}$（$m\gt n$）

特征值与奇异值问题

三、收敛判定与停止准则

三条停止准则，缺一不可： ① 相对残差（或相对增量）阈值 $\lVert\boldsymbol{r}_k\rVert/\lVert\boldsymbol{b}\rVert\le\varepsilon_{\mathrm{tol}}$；② 迭代次数上限（阈值没到也要停下来并报告）；③ 判据必须能算——用真解做判据等于没有判据。

实测收敛阶怎么估。 若误差满足 $e_{k+1}\approx C e_k^{p}$，则连续两个误差比给出

$$p\approx\frac{\ln!\left(e_{k+1}/e_k\right)}{\ln!\left(e_k/e_{k-1}\right)} .$$

可复算的例子（幂法求主特征值，矩阵与数据取自例 2.5.3：$A=\begin{pmatrix}5&2&0\2&4&1\0&1&3\end{pmatrix}$）：$\lvert\lambda_2/\lambda_1\rvert=0.5212$，实测 Rayleigh 商误差每步乘以 $0.2717$。两点读法：

• 比值的读法：$0.2717=0.5212^{2}$ 恰好说明「Rayleigh 商的误差按 $(\lambda_2/\lambda_1)^{2k}$ 衰减」——比只用 $\lVert A\boldsymbol{x}\rVert$ 估计（每步 $0.5212$）快一倍（定理 2.5.3）。
• 阶的读法：把实测误差代进上面的公式，得到 $p\approx0.98$——它告诉你的信息是「这个序列是几何衰减、比值 $0.27$」，而不是「一阶收敛」。报「阶」之前先确认序列真的落在 $e_{k+1}\approx Ce_k^{p}$（$p\gt1$）的形态里：割线法这类 $p\approx1.618$ 的序列（第 9 讲）才适合直接用三点公式估阶。

一份可以照抄的停止准则模板：

for k in range(maxit):                     # ② 上限永远写进循环
    x = step(x)                            # 迭代一步
    res = np.linalg.norm(b - A @ x) / bnorm # ① 相对残差
    if res < tol: break                    # 达标即停
else:
    print(f"未达标即触顶：k={maxit} 残差={res:.3e} "
          f"cond 估计={np.linalg.cond(A, 2):.2e}")  # ③ 触顶时报告诊断量

四、工具地图：任务 → 函数 → 后端

关于 numpy.linalg 的 LAPACK 后端。 NumPy 与 SciPy 调的是同一批 LAPACK / BLAS 例程（具体实现由发行版的 OpenBLAS、MKL 一类库提供），所以同一个矩阵在不同机器上可能出现末位差异——本教程各表因此都声明「逐位因平台的 BLAS 实现而异」。两处容易踩的差别：

• np.linalg.cond(A, 2) 走 SVD；而 cond(A, 1) / cond(A, np.inf) 按「$\lVert A\rVert\cdot\lVert A^{-1}\rVert$」实现，会显式求逆。大矩阵上想省这一步，用 SciPy 暴露的 gecon（先 LU 再做 $O(n^{2})$ 范数估计），调用签名随版本略变，用前查当前版本文档。
• np.linalg.svd 默认走 gesdd（分治），np.linalg.qr 走 geqrf；要换成别的例程（如 gesvd）得用 SciPy 的 svd(lapack_driver='gesvd')。
• 最小二乘的底层是 LAPACK 的 SVD 型求解器（gelsd 系列）。本机 numpy 2.4.6 的 lstsq 文档只承诺「用 SVD」，未点名例程——具体例程这一项需核实（要精确核对就查当前版本文档或源码）。

何时该用 scipy.sparse： 判据是「非零元占比」与「能不能只做矩阵–向量乘」两条。经验法则（社区惯例，不是定理）：$\mathrm{nnz}/n^{2}$ 低于约 $10$、且矩阵来自网格、图或文本这类天然稀疏结构时，换稀疏表示通常值得。三条实测提醒：① 稀疏格式要转对——spsolve 用 CSC（A.tocsc()）；② 直接法（splu）有 fill-in，分解后非零元可能比原矩阵多几个量级，内存是第一个瓶颈；③ 迭代法（cg / gmres）不收敛时先加预条件子（spilu 做不完全 LU），而不是一味加迭代次数。

五、三条「不要做」的事

1. 不要用显式求逆解方程组。

x = np.linalg.inv(A) @ b        # 不要这样写
x = np.linalg.solve(A, b)       # 正确的写法

• 代价账：solve 是「LU（$\frac23n^{3}$）+ 两次三角回代（$2n^{2}$）」；inv 还要由 LU 显式构造 $A^{-1}$（再加约 $n^{3}$），总计约 2.5 倍工作量。

• 精度账（可复算，$n$ 阶 Hilbert 矩阵、真解全为 $1$）：

  $n$	$\mathrm{cond}_2$	solve 最大误差	inv @ b 最大误差
  8	$1.53\times10^{10}$	$1.2\times10^{-7}$	$7.2\times10^{-7}$
  12	$1.66\times10^{16}$	$5.2\times10^{-1}$	$1.9\times10^{1}$

  显式求逆多引入一层舍入，误差被条件数再放大一次（1.3.6 的矩阵条件数是这件事的定量依据）。

• 同一 $A$ 多个右端时，正确做法不是 inv，而是分解一次、回代多次：lu_factor(A) + 循环 lu_solve。

2. 不要把 $A^{\mathrm T}A$ 送给特征值分解。

数学上 $A^{\mathrm T}A$ 的特征值确实给出奇异值的平方，但 $\mathrm{cond}(A^{\mathrm T}A)=\mathrm{cond}(A)^{2}$——有效数字直接掉一半。可复算的账：取 $12\times10$ 的 Vandermonde 设计矩阵（$x$ 均匀取自 $[0,1]$，np.vander(x, 10, increasing=True)），实测

$$\mathrm{cond}(A)=6.84\times10^{6},\qquad \mathrm{cond}(A^{\mathrm T}A)=4.68\times10^{13},\qquad \text{比值}=6.84\times10^{6}=\mathrm{cond}(A).$$

所以：算 SVD 用 np.linalg.svd(A)；解最小二乘用 lstsq / QR；算条件数用 cond(A, 2)。eigh(A.T @ A) 只适合演示「条件数为什么会平方」（2.6.7 常见错误 2、2.8.2 陷阱 5）。

3. 不要把默认阈值当数据的分辨率。

• np.linalg.matrix_rank 的默认阈值是 $\sigma_1\cdot\max(m,n)\cdot\mathrm{eps}$（本机 numpy 2.4.6 文档如此），它只适合「精确矩阵、数值上降秩」的场景。带噪数据的实测：给一个数值秩 3 的矩阵加 $10^{-9}$ 量级噪声后，默认阈值判出 $\operatorname{rank}=20$，按噪声水平截断得 $3$（实验三 (d)）。
• 同一条也适用于伪逆：np.linalg.pinv 的默认 rcond 是 $10^{-15}$，不截断的伪逆会把 $1/\sigma_{\min}$ 放大 $\mathrm{cond}(A)$ 倍，与直接求逆同样危险（2.6.4 思考）。正确做法：按数据噪声水平给 rcond（或 tol），并报告所选阈值。

六、验收清单（提交数值结果前过一遍）

七、症状 → 先看哪一页 / 用哪个函数

依据第一章 + 补充篇正文的公式、量级与实验整理 · 本页所有函数名与默认值在本机 numpy 2.4.6 / scipy 1.18.0 上核对过；各版本默认阈值会变，动手前查当前版本文档 · 与先修与自测的分工：那一页管基础，本页管操作 · 整理日期 2026-10-20。