简介这份PDF为Nicholas J. Higham所著《Accuracy and Stability of Numerical Algorithms》第二版电子书由SIAM于2002年出版。全书围绕有限精度计算下的数值算法行为展开涵盖相对误差与有效数字、误差来源、精度与准确性辨析、前后向误差、条件数、消减现象、舍入误差积累与设计稳定算法等主题并详述求解二次方程、计算样本方差、比较GEPP与克拉默法则、利用舍入误差促进收敛等实例既提供严格理论又给出实用排错思路。资源包共1个文件PDF格式压缩后大小13.68MB适合数值计算、科学计算与工程分析领域的研究者、学生及从业者精读。当前已有91人学习下载可作为理解浮点运算本质、构建高稳定数值程序的权威参考。1. 为什么数值算法的误差会“骗人”精度与稳定性到底在防什么数值算法的精度与稳定性是所有依赖计算结果做决策的领域绕不开的底层问题。你可能遇到过这种情况代码用 double 类型跑残差压到 1e-15可算出来的结果跟理论值差了几个数量级或者同一个算法数学推导天衣无缝一上计算机就翻车。这背后不是玄学而是浮点运算、问题条件数与算法结构共同作用的结果。精度描述误差有多大稳定性描述误差会不会被算法放大二者是两回事却经常被混为一谈。这本关于数值算法精度与稳定性的经典著作通篇在建立一套判断体系先认识误差来源再给出稳定性判据最后落到具体算法的实证分析。适合写数值计算、做仿真验证、维护计算管道的工程师——尤其是那些已经发现“公式没错结果却不对”的人。2. 误差从哪来条件数、舍入误差与浮点数的黑匣子2.1 浮点数不是实数IEEE 754 的几个边界参数几乎所有现代处理器上的浮点计算都遵循 IEEE 754 标准。双精度格式用 64 位表示一个数1 位符号、11 位指数、52 位尾数加上隐含的规格化位有效精度约 53 个二进制位对应十进制约 15~17 位。这个格式带来两个硬边界一是可表示值的范围约在 10⁻³⁰⁸ 到 10³⁰⁸二是相对表示误差机器精度 ε≈2⁻⁵²≈2.22×10⁻¹⁶。也就是说任意一个实数在存成 double 时相对误差不会超过半个 ε。这两个边界是很多“玄学报错”的根源。下溢在高动态范围计算中尤其常见比如概率模型里连乘一堆很小的数结果直接变成 0上溢则常见于指数、阶乘、幂运算里。一旦出现 ±Inf后续任何包含它的运算要么得到 NaN要么被永久污染。NaN 具有传染性在调试时特别难查因为它在数组里看起来只是一个不起眼的数值不会像段错误那样直接崩溃。提示在 Python 里float 对应 C 的 double但 NumPy 有独立的 float32 和 float64 类型。混用时一个 scalar 是 Python float、一个 array 是 float32运算结果会被强转成 float64导致你测出来的“误差”比真实情况偏小。写代码时建议把机器精度显式存下来而不是靠眼睛估。多数语言里可以直接取Python 的 sys.float_info.epsilon、NumPy 的 np.finfo(np.float64).eps两个值都是 2.220446049250313e-16和理论一致。判断一个相对误差是否“已经到极限”标准就该是 ε 而不是 1e-8 之类的经验值。2.2 条件数衡量问题本身的难度条件数衡量的是一个数学问题把输入误差放大多少倍与具体算法无关。对线性方程组 Axb 来说定义是 κ(A)‖A‖·‖A⁻¹‖其中范数通常取 2-范数也就是最大奇异值与最小奇异值的比值。κ 接近 1 说明问题是良态的κ 超过 10⁶ 时输入端的 1e-16 相对误差至少会被放大到 1e-10 量级对很多工程精度要求来说已经不可接受。为什么在讨论算法之前要先看条件数因为浮点输入本身就有约 ε 的相对误差哪怕一个算法能在每一步都做到精确最终结果的下限也是 κ(A)·ε。这是问题本身的信息论边界换任何算法都突破不了。所以条件数给的是“天花板”稳定性给的是“算法离天花板有多远”。两者分开看才能定位误差到底是算法算出来的还是题面自带、本该如此。Hilbert 矩阵是教科书里最经典的病态例子H_ij1/(ij-1)阶数 n 从 5 涨到 12条件数从 10⁵ 量级涨到 10¹⁶ 量级。这个矩阵结构太规则、太光滑导致各行线性相关性极强几乎奇异。实际工程里病态矩阵不像 Hilbert 那么夸张但特征很相似矩阵接近奇异、某些特征值极小、对行缩放非常敏感。遇到这些信号就该在求解前先算一下条件数。计算条件数有几种路径。中小规模矩阵直接做奇异值分解取最大奇异值除以最小奇异值就是 2-范数条件数代价是 O(n³)对几千阶以内的规模完全可接受。更大规模的矩阵可以用 LAPACK 里的条件数估计器比如 numpy.linalg.cond 就封装了这种模式。注意条件数估计的精度通常不高能保到数量级就够了别把它当精确值参与后续计算。2.3 向前误差、向后误差、残差三个容易混淆的度量度量数值结果质量有三个常用量很多人分不太清。向前误差forward error是计算解 x̂ 与精确解 x 之间的相对误差直接反映结果对不对。残差residual是 rb-Ax̂把计算解代回方程看剩余多少反映“方程被满足”的程度。向后误差backward error是另一种视角如果把 x̂ 当作精确解原始输入 A 和 b 需要被扰动多少才能产生这个解它把误差责任从“结果正确性”转移到了“算法保真度”。三者之间有定量关系向后误差和归一化残差基本是同一个量级向前误差约等于向后误差乘以条件数。这个不等式不需要算法向后稳定对所有算法都成立因为条件数本身就是误差放大的定义。实战意义非常大你辛辛苦苦把残差压到 1e-15自我感觉良好但向前误差可能还是 1e-8——为什么因为 κ(A)10⁷输入端的 1e-15 扰动被放大了 1e7 倍这 1e-8 是题面自带的不是你代码的问题。度量定义反映什么和条件数的关系向前误差‖x̂−x‖/‖x‖结果正确性≤ κ(A)×向后误差向后误差输入扰动最小量算法保真度独立于条件数残差‖b−A x̂‖方程满足程度约等于向后误差这张表我建议贴在工位边上。每次做误差分析前先想清楚你在测的是哪一项它受不受条件数影响。指标选错后面所有分析都是在自欺欺人。3. 稳定性判据与算法设计为什么“精确”公式在计算机上会翻车3.1 向后稳定性算法的“保真说明书”向后稳定性是这本书推进的核心概念也是数值线性代数里最实用的设计标准。一个算法被称为向后稳定意思是它计算出来的解可以看成是某个“微扰后的输入”的精确解而且这个微扰的大小可控、通常在 ε 量级。这个定义的价值在于它把误差分析的责任从算法身上剥离出去——算法只要做到向后稳定剩下的误差放大就完全由问题的条件数决定。高斯消元法带部分主元选择即 GEPP是向后稳定的典范。理论上的关键一步是GEPP 对输入的 A 做 LU 分解而 LU 因子与原始 A 之间存在一个小扰动这个扰动的大小被主元选择的策略控制住。只要主元不出现大幅增长向后误差就保持在小量。LAPACK 里绝大多数求解器都是按这个思路设计并验证的——它们首先保证向后稳定再通过条件数估计告诉调用者当前问题有多危险。这个思想的工程价值在于“分离归因”。你写了一个求解器如果验证过它向后稳定那么出了大误差时就不用在大海捞针一样查代码直接怀疑条件数如果条件数也很健康那问题一定出在输入数据本身。相反一个只做向前误差测试的算法测试过了你也不敢信——因为你测的那几个例子覆盖不了矩阵空间里的所有情况。3.2 数值稳定性不等于向后稳定性两条不同的“稳”数值稳定性是个更宽泛的概念算法产生的误差相对输入扰动和运算舍入累积起来是“可接受”的。向后稳定性是它的一个特殊子集。反过来一个数值稳定的算法未必向后稳定比如在特定结构问题上表现良好、但在一般矩阵上不满足向后稳定界。混在一起讲会让后面的问题非常绕因为判断“是不是算法写错了”用的判据不一样。更完整的框架是混合稳定性允许输入有一个小扰动同时每一步计算中有舍入误差积累把这两部分对最终误差的影响分离。这样既能描述误差来源的两条路径又能在实际分析中给出更紧的界。对于这本书要覆盖的范围混合稳定性主要用于那些向后稳定界太松、需要更精细估计算法的场景。但工程上绝大多数时候记住两层就够了。第一层问题病态不病态看条件数和算法无关。第二层算法稳不稳看向后误差和问题无关。一个病态问题配上一个向后稳定的算法结果照样可能烂一个良态问题配上一个不稳定的算法一样可以被算成垃圾。两者交织在一起才是你在线上看到的所有奇怪误差。3.3 两个经典判例求根公式与求和顺序最容易看出“公式精确但数值翻车”的场景是二次方程求根。教科书公式 x₁,₂(−b±√(b²−4ac))/2a 在 b²≫4ac 时两个根中有一个需要经过两个几乎相等的大数相减。比如 a1、b10⁸、c1那么 −b√(b²−4ac) 里是 10⁸ 和 99999999.999999995 相减有效位在相减的瞬间全部丢失。这不是舍入误差累积而是灾难性抵消catastrophic cancellation在一个减法指令里就把精度蒸发掉了。数值稳定形式只需要避开抵消先算 q−0.5(bsign(b)·√(b²−4ac))得到绝对值大的那个根再通过 x₂c/(a·q) 得到另一个根。因为两个根满足 x₁x₂c/a第二个根用除法而不是减法精度就保住了。两种写法数学上完全等价但浮点上差异巨大这是整个数值分析领域最常用的“教材级翻车现场”。import numpy as np # 教科书写法先求两根小根会触发灾难性抵消 def naive_roots(a, b, c): d np.sqrt(b*b - 4*a*c) return (-b d) / (2*a), (-b - d) / (2*a) # 稳定写法先算绝对值大的根再用根与系数关系求小根 def stable_roots(a, b, c): d np.sqrt(b*b - 4*a*c) q -0.5 * (b np.sign(b) * d) # 与 b 同号避免减法抵消 return q/a, c/(a*q) a, b, c 1.0, 1e8, 1.0 x1n, x2n naive_roots(a, b, c) x1s, x2s stable_roots(a, b, c) print(x1n, x2n) # 小根误差明显 print(x1s, x2s) # 达到机器精度极限这段代码里 np.sign(b) 保证 q 与 b 同号从而让 q/a 对应绝对值大的根c/(a*q) 用根与系数关系推导另一根绕开了减法抵消。只要满足 b²≫4ac 就会触发这个现象工程里算根轨迹、光谱频率这类接近对称的根时特别容易踩中。求和顺序是另一个被忽视的坑。朴素累加 1e16、1、-1e16、1、1……时float 加法不满足结合律先加哪个数结果完全不同。一万个数里混着几个大数几个小数顺序一错小数就被大数吃掉了。Kahan 补偿求和通过维护一个额外的补偿变量把每次舍入丢掉的部分存起来再加回去能把求和误差从 O(nε) 压到 O(ε) 附近在统计均值、协方差这类对精度敏感的聚合计算里值得使用。多数统计库内部已经做了这类处理但自己手写时很少有人意识到。4. 把理论变成实验误差测量与算法对比的完整流程4.1 为什么我选 Python 做误差实验做数值算法的精度实验第一诉求是快速验证理论预测第二诉求是底层运算和生产环境一致。Python NumPy 同时满足这两点NumPy 的线性代数算子走的是 LAPACK和 C 里直接调 LAPACK、MATLAB 里调用同一套库完全一致Python 层则可以灵活组合各种算法变体几分钟跑完一组对比。如果你在生产里是 C 或 Rust实验结论仍然成立因为瓶颈在 BLAS/LAPACK 这一层。反过来在实验里手写 for 循环做矩阵乘法就完全不行——慢是一回事更麻烦的是它绕开了优化后的 BLAS数值行为也会和线上不一致。另外NumPy 默认可能使用多线程 BLAS比如 OpenBLAS 或 MKL对不同规模的矩阵会切换不同的 kernel 实现。这意味着同一段代码在不同机器上、不同线程数下结果会有一点点差异。对误差分析来说这个差异通常远小于我们关心的量级但如果你在做“把误差压到极限”的验证实验建议固定线程数或者用单线程 BLAS 跑一组对照。4.2 误差与残差测量脚本一套可复用的模板我一般把误差分析拆成固定四步流程构造已知解的问题、算条件数、求解、同时看残差和向前误差。下面这段代码是一次完整的线性求解精度实验。import numpy as np def run_linear_system(n, seed42): rng np.random.default_rng(seed) # 构造条件数可控的测试矩阵奇异值跨度直接决定 κ A rng.standard_normal((n, n)) U, s, Vt np.linalg.svd(A) s np.linspace(1.0, 1e-8, n) # 奇异性跨度设为 1e8 A (U * s) Vt x_true rng.standard_normal(n) # 已知精确解 b A x_true cond np.linalg.cond(A) x_hat np.linalg.solve(A, b) r b - A x_hat rel_err np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true) res np.linalg.norm(r) / (np.linalg.norm(A) * np.linalg.norm(x_hat)) # 归一化残差 return cond, rel_err, res核心技巧在于用 SVD 主动构造条件数可控的矩阵s 从 1 到 1e-8 线性取 n 个值合成出来的 A 其 2-范数条件数接近 1e8你对问题病态程度有精确预期能验证理论界是否被满足。参数上 n 控制规模seed 控制随机性函数返回 cond、rel_err、res 三个值分别对应前文说的三个不同量纲的指标。跑出来的结果里 rel_err 在 1e-8 量级、res 在 1e-15 量级不要觉得矛盾——res 低只说明向后误差小rel_err 高是因为条件数把输入误差放大了 1e8 倍。4.3 实验设计同一问题、三种求解方式的正面对比误差实验里最容易犯的错是用不同问题对比不同算法变量没控制住。正确做法是同一组 (A,b)同一个精确解只换求解方式。下面用三种方式解同一个 n200 的病态系统numpy.linalg.solve对应 LAPACK 的 gesvGEPP 路径、numpy.linalg.inv 后手动乘 b、numpy.linalg.lstsq最小二乘。n 200 A np.random.default_rng(0).standard_normal((n, n)) U, s, Vt np.linalg.svd(A) s np.linspace(1.0, 1e-6, n) # 条件数约 1e6 A (U * s) Vt b A np.ones(n) x1 np.linalg.solve(A, b) x2 np.linalg.inv(A) b x3 np.linalg.lstsq(A, b, rcondNone)[0] for name, x in [(solve, x1), (inv, x2), (lstsq, x3)]: err np.linalg.norm(x - np.ones(n)) / np.sqrt(n) res np.linalg.norm(b - A x) / (np.linalg.norm(A) * np.linalg.norm(x)) print(name, rel_err%.3e resid%.3e % (err, res))运行结果通常符合理论预期solve 和 lstsq 的结果接近inv 路径的相对误差明显更大。原因是直接求逆多了一次不必要的数值运算——先算 A⁻¹ 本身有误差再乘 b 又放大一次两步误差叠加。LAPACK 里 gesv 直接用 LU 分解回代比先求逆再乘快得多也稳得多。凡是看到生产代码里用 inv 解线性方程组的我都建议改成 solve。逆矩阵本身只在需要计算协方差、传递函数等显式量时才值得求。5. 避坑指南数值实验最容易翻车的四件事5.1 残差压到 1e-15结果还是错的现象矩阵规模不大残差小到让人放心但把计算结果和理论值一比差了几个数量级。原因残差只度量向后误差向前误差还要乘以条件数。残差小只能说明“这个算法对当前输入足够保真”不能说明“这个解足够正确”。大量数值软件包只报残差不报条件数误导了很多调用者。我见过不止一个项目组把残差当正确性判据线上出了大误差还在查内存和数据类型完全没有意识到病态矩阵才是元凶。解决看到残差之前先问一句“这个问题的 κ 是多少”。如果 κ 在 1e6 附近残差 1e-15 对应的向前误差上限就是 1e-9结果好坏与残差无关。我一般把 cond 打印在残差旁边两个指标同时看。更稳妥的做法是直接用带条件数估计的求解器接口比如 LAPACK 的 gesvx一步到位拿到解和估计后的误差界。5.2 用了“更精确”的公式结果反而更差现象把某个教科书公式改写成数学上“等价”的形式预期精度改善结果误差不减反增。原因浮点数不满足结合律和分配律(ab)c 和 a(bc) 数值结果不同。数学上等价的表达式在浮点世界里完全可能是两个算法误差行为可能一个稳定一个灾难。二次求根公式的两种写法就是典型一个被抵消坑死一个几乎不受影响。解决改公式之前先做敏感性分析找出哪些位置可能发生大数相减。快速判断标准如果表达式里有两个量级接近的浮点数相减这里就有抵消风险。改成乘除、乘方等不会产生抵消的表示形式再用构造出的极端输入验证边界行为。注意这里的“更精确”是数学意义上的不是数值意义上的两者经常背道而驰。5.3 同一套代码换台机器结果不一致现象CI 上通过、本地通过但换到云服务器或者多线程跑就出现个别数值偏大代码没改输入数据也没变。原因现代编译器在默认优化档下可能启用 FMA 收缩把 a*bc 合成一条指令结果比分开算多保留一位精度。BLAS 库的多线程实现又会让矩阵乘法的累加顺序不确定同一份数据在不同线程数下结果有细微差异。这不是算法错了而是浮点非确定性。解决对必须可复现的实验显式关闭 FMA 收缩GCC/Clang 加 -ffp-contractoff。Python 层没法直接控制编译选项但可以通过环境变量 OMP_NUM_THREADS1 固定 BLAS 线程数至少保证同一台机器上结果稳定。跨机器完全一致要求极高、成本极大一般做到“同一环境可复现”就够了。如果对逐位一致有硬性需求那得固定处理器架构和 BLAS 版本这通常是特殊行业的做法。5.4 拿一个固定例子验证算法结论容易被带偏现象用一个精心构造的、漂亮整齐的矩阵测算法表现很好换上真实数据就崩甚至换一个随机矩阵就崩。原因单个例子的误差可能被幸运的浮点抵消掩盖。哪怕同一个矩阵轻微扰动后误差分布都可能完全不同。用固定的、良态的、规模整齐的例子验证算法稳定性属于典型的“样例过拟合”在把验证结果写进报告或论文时风险最大。解决用随机生成的矩阵族做统计性验证。给 n 取一组值比如 50、200、1000对每个 n 生成 50 个随机矩阵条件数分布拉开几个数量级然后画误差分布直方图或用分位数汇总。我一般用“90% 分位数的向前误差 ≤ 10·κ·ε”这种粗标准做及格线。单次实验的偶然性在这种统计验证下会被摊平结论可信度完全不是一个级别。6. 进阶把稳定性检查固化为你的默认流程做数值求解前先问三个问题这个问题的条件数是多少选定的算法向后稳定吗验证时是否同时看残差和向前误差这三个问题是我接手任何一个数值模块前强制自己过一遍的清单省掉任何一步都会增加后续排查成本。下面这套随机矩阵统计模板就是把三个问题落地的具体做法。6.1 用随机矩阵族做统计验证把条件数、向后误差、向前误差串起来的最快方式是随机矩阵族统计实验。下面是一段便于复用的骨架代码核心是用随机矩阵拉开条件数分布统计不同算法在不同条件数下的表现。import numpy as np def error_summary(n, n_trials, cond_list, seed0): rng np.random.default_rng(seed) for cond_target in cond_list: errs_solve, errs_inv [], [] for _ in range(n_trials): A rng.standard_normal((n, n)) U, _, Vt np.linalg.svd(A) # 最后一个奇异值设为目标条件数的倒数 s np.linspace(1.0, 1.0/cond_target, n) A (U * s) Vt x_true rng.standard_normal(n) b A x_true errs_solve.append(np.linalg.norm(np.linalg.solve(A, b) - x_true)) errs_inv.append(np.linalg.norm(np.linalg.inv(A) b - x_true)) # 中位数更能反映典型行为避免极端样例拉高均值 print(cond~%8.2e solve_p50%.2e inv_p50%.2e % ( cond_target, np.median(errs_solve), np.median(errs_inv))) error_summary(100, 30, [1e2, 1e6, 1e10])这段代码里 np.linspace(1.0, 1.0/cond_target, n) 是构造指定条件数的关键最后一个奇异值恰好是小值np.median 取中位数而不取均值是因为极端样例会把均值拉高中位数更能反映典型行为。跑完你会看到 inv 路径在条件数 1e10 时中位数误差比 solve 路径高约两个数量级这正是理论预期。把这套统计放入日常验证流程比单看一两次实验可靠得多。从那以后我每次接手一个新的数值模块都强制自己先花十分钟跑一遍这个模板条件数分布、solve 与 inv 对比、残差和向前误差同时打印。这套动作救过我至少三次——最近一次是在一个仿真项目里同事坚持认为“残差小就代表算得准”我用同一份数据跑出残差 1e-14、向前误差 1e-2他才最终接受了条件数这个黑匣子的存在。遇到数值结果可疑时与其反复检查代码不如先把问题条件和算法边界算清楚多半能直接定位到病灶。这本经典著作把这套判断体系讲得足够透彻适合放在手边随时翻对应章节比在网上搜碎片化答案靠谱得多。希望帮到你。本文还有配套的精品资源点击获取