简介：《Accuracy and Stability of Numerical Algorithms》第二版由曼彻斯特大学数值计算权威 Nicholas J. Higham 撰写、SIAM 出版，面向科学计算、工程计算与数据分析领域的学习者和研究者，系统讲解有限精度环境下数值算法的准确性与稳定性评估方法。内容从相对误差、有效数字、误差来源、精度与准确性等基础概念切入，深入剖析前向与后向误差、条件数、消减现象，并以解二次方程、样本方差计算、线性方程组求解为例，对比 GEPP 与克拉默法则的稳定性与效率，还讨论了舍入误差的累积规律、非随机性及其有益作用，以及算法稳定性对问题的依赖性。资源包为 1 个 pdf 文件，压缩包约 13.68MB，完整呈现原书正文、目录与参考文献，便于按章节系统研读。目前已有 91 人学习下载，适合希望夯实数值分析理论基础、提升算法设计稳健性的读者参考。

1. 数值算法的精度与稳定性：为什么你的计算结果总在“漂移”

很多工程师第一次被数值稳定性“咬”到，是在一个看起来毫无风险的场景里：用双精度浮点累加一百万个量级相近的数，结果和理论值差了百分之几。代码逻辑没错，公式推导也没错，错的是我们默认了“浮点运算等于实数运算”。Accuracy and Stability of Numerical Algorithms 这个标题讲的正是这件事——算法的精度（算出来的结果离真值有多近）和稳定性（输入或中间步骤的微小扰动会不会被放大）是两套独立的评价体系，一个算法可以精度很高但极不稳定，也可以稳定但精度受限于问题本身的条件数。

这篇文章面向的是需要把数值方法落到生产代码里的从业者：做信号处理、做仿真、做机器学习底层算子、做金融定价、做图形几何的人都会撞上它。我会从误差来源讲起，落到可复现的代码实验，再给出参数选择和排查套路。读完你应该能判断自己手里的算法是“玄学漂移”还是“结构性问题”，并且知道怎么改。

2. 误差从哪来：浮点模型、条件数与算法稳定性的三角关系

2.1 浮点数的真实行为：不是精度不够，是分布不均匀

IEEE 754 双精度有 53 位有效二进制位，机器精度 eps 约为 2.22e-16。但关键不是这个数字有多小，而是浮点数在数轴上的分布是非均匀的：越靠近零越密，越大越稀疏。两个量级相差超过 2^53 的数相加，小的那个会被完全“吃掉”。这不是 bug，是浮点加法的定义。

我一般用下面这段代码让新人直观感受一下：

import numpy as np # 场景一：大数吃小数 a = 1.0 b = 1e-16 print(a + b == a) # True，b 被完全吞掉 # 场景二：累加顺序影响结果 vals = np.random.rand(1_000_000) s1 = 0.0 for v in vals: s1 += v # 顺序累加 s2 = np.sum(vals) # 成对求和 print(abs(s1 - s2)) # 通常能差到 1e-10 量级

逻辑说明：第一段验证浮点加法的吸收性，第二段对比朴素累加和 NumPy 的成对求和（pairwise summation）。参数上，1_000_000这个规模足以让误差累积到肉眼可见，但又不至于跑太久。如果你把np.sum换成math.fsum，会得到精确舍入的结果，代价是速度慢几倍。

提示：判断一个累加是否安全，先看数据量级跨度。如果 max/min 超过 1e15，朴素累加基本不可信。

2.2 条件数：问题本身的“病态程度”

条件数衡量的是输入扰动对输出的放大倍数。对于线性方程组 Ax=b，条件数 cond(A) = ||A||·||A⁻¹||。条件数为 10^k 意味着输入有 k 位有效数字的损失。这是问题固有的属性，跟你用什么算法无关。

import numpy as np # 构造一个 Hilbert 矩阵，经典病态例子 n = 12 H = np.array([[1.0/(i+j+1) for j in range(n)] for i in range(n)]) print("cond(H) =", np.linalg.cond(H)) # 量级约 1e16 b = H @ np.ones(n) x = np.linalg.solve(H, b) print("max error =", np.max(np.abs(x - 1.0)))

逻辑说明：Hilbert 矩阵的条件数随 n 指数增长，n=12 时已经接近双精度极限。np.linalg.solve用的是 LU 分解，算法本身是稳定的（后向稳定），但问题病态导致前向误差巨大。参数上，n 超过 13 基本就没法用双精度直接解了，需要换正则化或高精度库。

2.3 稳定性：算法会不会放大已有的误差

稳定性分前向稳定和后向稳定。后向稳定是说：算法给出的解，是某个微小扰动后问题的精确解。LU 分解、QR 分解、Householder 变换都是后向稳定的。而经典 Gram-Schmidt 正交化只是前向稳定，改进版（MGS）才后向稳定。

选型理由很简单：优先选后向稳定算法，再考虑条件数。如果问题本身病态，任何算法都救不了前向精度，但后向稳定算法至少保证你不会雪上加霜。

3. 把稳定性落到代码：四个可复现的实验与修法

3.1 实验一：矩阵乘法的累加顺序

矩阵乘法 C = A·B，每个元素是 A 的一行和 B 的一列做点积。点积的累加顺序直接影响误差。

import numpy as np np.random.seed(0) n = 2000 A = np.random.randn(n, n) B = np.random.randn(n, n) # 方法一：NumPy 默认（底层 BLAS，分块累加） C1 = A @ B # 方法二：手动三重循环（朴素累加） C2 = np.zeros((n, n)) for i in range(n): for j in range(n): s = 0.0 for k in range(n): s += A[i, k] * B[k, j] C2[i, j] = s print("max diff =", np.max(np.abs(C1 - C2)))

逻辑说明：BLAS 内部用分块和向量化累加，误差远小于朴素三重循环。参数上 n=2000 时手动循环会跑几十秒，实际调试可以降到 200。这个实验的意义是：不要自己写矩阵乘法，除非你在做 BLAS 本身。

3.2 实验二：求解线性方程组的三种姿势

import numpy as np from scipy.linalg import lu, qr np.random.seed(1) n = 500 A = np.random.randn(n, n) + n * np.eye(n) # 对角占优，条件数可控 b = np.random.randn(n) # 姿势一：直接求逆（最差） x_inv = np.linalg.inv(A) @ b # 姿势二：LU 分解求解 x_lu = np.linalg.solve(A, b) # 姿势三：QR 分解求解 Q, R = np.linalg.qr(A) x_qr = np.linalg.solve(R, Q.T @ b) for name, x in [("inv", x_inv), ("lu", x_lu), ("qr", x_qr)]: print(name, "residual =", np.linalg.norm(A @ x - b))

逻辑说明：求逆再乘的残差通常最大，因为显式求逆引入了额外误差。LU 和 QR 残差接近，但 QR 在条件数较大时更稳。参数上，n * np.eye(n)里的系数控制对角占优程度，系数越小条件数越大，三种方法的差距越明显。

注意：np.linalg.solve内部会根据矩阵结构选择 LU 或 Cholesky，不要假设它一定用 LU。

3.3 实验三：特征值分解的对称性陷阱

对称矩阵用eigh，非对称用eig，混用会导致结果不可复现。

import numpy as np np.random.seed(2) n = 300 S = np.random.randn(n, n) S = (S + S.T) / 2 # 强制对称 # 正确：对称专用 w1, v1 = np.linalg.eigh(S) # 错误：通用算法 w2, v2 = np.linalg.eig(S) print("eigh 虚部 max =", np.max(np.abs(np.imag(v1)))) print("eig 虚部 max =", np.max(np.abs(np.imag(v2)))) print("特征值差异 =", np.max(np.abs(np.sort(w1) - np.sort(np.real(w2)))))

逻辑说明：eigh利用对称性，保证实特征值和正交特征向量，速度快一倍以上。eig对对称矩阵也能算，但数值路径不同，可能产生微小虚部。参数上，n=300 时差异还小，n 上千后eig的虚部会明显起来。

3.4 实验四：求和与方差的稳定算法

方差计算用E[x²] - E[x]²是经典翻车点，当均值远大于标准差时会灾难性抵消。

import numpy as np np.random.seed(3) x = 1e8 + np.random.randn(1_000_000) # 均值 1e8，标准差 1 # 不稳定：直接公式 var_naive = np.mean(x**2) - np.mean(x)**2 # 稳定：Welford 在线算法 mean = 0.0 M2 = 0.0 for i, v in enumerate(x): delta = v - mean mean += delta / (i + 1) M2 += delta * (v - mean) var_welford = M2 / len(x) print("naive =", var_naive) print("welford =", var_welford) print("numpy =", np.var(x))

逻辑说明：朴素公式里mean(x**2)和mean(x)**2都是 1e16 量级，相减后有效数字几乎全丢。Welford 算法每步只处理增量，避免了大数相减。参数上，1e8这个均值让抵消效应非常明显，实际数据里均值/标准差比值超过 1e8 就要警惕。

4. 避坑与排查：数值不稳定问题的五个典型现场

4.1 现象：结果每次运行都不一样，但代码没有随机性

原因：多线程 BLAS 的归约顺序不确定，或者编译器开了-ffast-math重排了浮点运算。解决：固定线程数（OMP_NUM_THREADS=1），编译时去掉-ffast-math，或者接受这种非确定性并在容差内比较。

4.2 现象：小规模数据正常，大规模数据结果发散

原因：误差随问题规模累积，比如迭代法的迭代次数不够，或者矩阵条件数随 n 增长。解决：先算条件数，再看残差随 n 的变化曲线。如果残差斜率异常，换后向稳定算法或加预条件。

4.3 现象：GPU 上和 CPU 上结果对不上

原因：GPU 的归约顺序、FMA 指令使用、以及 float32 默认精度都和 CPU 不同。解决：统一用 float64 做验证，或者把容差放宽到 float32 的 eps 量级（约 1e-7）。不要指望 bitwise 一致。

4.4 现象：稀疏矩阵求解比稠密还慢

原因：填充（fill-in）导致稀疏结构被破坏，或者用了不适合的预条件子。解决：换 AMD/COLAMD 重排序，或者改用迭代法（CG、GMRES）配对角预条件。参数上，先看nnz(LU)/nnz(A)这个比值，超过 10 就要考虑换方法。

4.5 现象：特征值出现负的或复数，但矩阵明明是对称正定的

原因：数值误差让最小的特征值跨过零，或者用了eig而不是eigh。解决：对称矩阵一律eigh，正定性检查用 Cholesky 分解是否成功，而不是看特征值符号。

5. 进阶：用混合精度和迭代精化把精度“买”回来

混合精度是现在 GPU 上的主流做法：用 float16 或 bfloat16 做矩阵乘法，用 float32 或 float64 做累加。核心思想是把误差控制在累加环节，因为乘法本身的相对误差是固定的。

import numpy as np import torch # 模拟混合精度矩阵乘法 n = 4096 A = torch.randn(n, n, dtype=torch.float32) B = torch.randn(n, n, dtype=torch.float32) # 纯 float32 C32 = A @ B # 混合精度：转 float16 乘，float32 累加 A16 = A.half() B16 = B.half() C16 = (A16 @ B16).float() print("相对误差 =", torch.norm(C32 - C16) / torch.norm(C32))

逻辑说明：A16 @ B16在 GPU 上走 Tensor Core，速度快但精度低。转回 float32 后和纯 float32 结果比较，相对误差通常在 1e-3 量级。参数上，n=4096 是 Tensor Core 的典型工作尺寸，太小则加速比不明显。

迭代精化是另一个技巧：解出 x 后算残差 r = b - Ax，再用 r 解一次修正量，更新 x。对于条件数在 1e8 以内的矩阵，两到三轮迭代精化能把前向误差降到机器精度。

import numpy as np def iterative_refinement(A, b, max_iter=3): x = np.linalg.solve(A, b) for _ in range(max_iter): r = b - A @ x dx = np.linalg.solve(A, r) x = x + dx if np.linalg.norm(dx) < 1e-14 * np.linalg.norm(x): break return x

逻辑说明：每轮用当前残差解一个修正方程，逐步逼近。参数上，max_iter=3对多数问题够用，1e-14是收敛阈值。注意残差要用高精度算，否则修正量本身就不准。

我自己的习惯是：任何涉及矩阵分解的代码，先跑一遍条件数估计，再决定用不用迭代精化。条件数超过 1e12 就直接上高精度库，别在双精度里硬扛。数值算法这件事，后悔药很少，前期多花十分钟做误差分析，比后期调一周的“玄学 bug”划算得多。希望帮到你。

本文还有配套的精品资源，点击获取