简介《Accuracy and Stability of Numerical Algorithms》是数值计算领域权威学者Nicholas J. Higham的经典著作第二版由SIAM在费城出版面向科学计算、工程分析与数据科学读者系统解答有限精度计算下“算法为什么出错、如何保证可靠”的核心问题。全书从相对误差、有效数字谈起梳理误差来源与前/后向误差的区别用条件数刻画问题敏感性并以二次方程求解、样本方差计算、线性方程组求解等实例说明消减、舍入误差积累及不稳定性带来的陷阱给出稳定算法设计和提高精度的改进路径。资源为1个PDF电子书文件约13.68MB可直接在电脑、平板等设备阅读目前已有91人学习。书中包含大量实例、练习与文献指引章节自有限精度计算原理延伸至浮点算术、误差分析、矩阵计算等主题目录清晰既适合入门者建立数值算法的误差分析框架也适合进阶者按需查阅具体结论是科学计算领域难得的常备参考书。1. 数值算法的精度与稳定性为什么公式没问题算出来却不敢信做控制系统辨识时我曾把一张六阶矩阵交给同事的求解器行列式不为零、残差只有 1e-10参数却让仿真直接发散。问题不在建模而在“数值算法的精度与稳定性”这个步骤被完全忽略了。有限位浮点数下任何一个加减乘除都可能让误差放大几个数量级精度决定单步误差有多大稳定性决定误差在传播中是逐渐衰减还是彻底爆炸。这篇文章我按实际排查的顺序写先找误差来源再给量化手段最后落到矩阵分解重构、参数调整与避坑适合正在写数值核心、做控制和仿真的工程师照着走一遍。2. 精度误差从哪里来机器精度、舍入误差与灾难性消去2.1 浮点数不是实数机器精度与有效数字的边界在数学上1 和 1e-20 相加还是 1.000…0001在计算机里double 的有效数字只有 1516 位1 1e-20 会被直接舍入成 1。这个基本事实划定了数值算法精度的天花板也解释了为什么“公式正确”不等于“结果正确”。IEEE 754 下几种常用格式的关键参数如下类型尾数有效位机器精度 eps最小正规数十进制有效位数float3223 bit约 1.19e-7约 1.18e-38约 7 位float6452 bit约 2.22e-16约 2.23e-308约 1516 位用np.finfo可以亲手验证import numpy as np # 看 double 和 float 的机器精度 print(np.finfo(np.float64).eps) # 2.220446049250313e-16 print(np.finfo(np.float32).eps) # 1.1920929e-07逻辑说明机器精度 eps 指“1 与下一个可表示的浮点数之间的相对距离”。单次浮点运算通常最多引入 eps / 2 的相对误差这个量叫单位舍入 u。工程上常说的“double 能保证 15 位有效数字、float 只能保证 7 位”就是从这里来的。参数说明如果中间量本身量级接近 1e8而最终结果又需要保留到 1e-8你实际上已经快把 double 的全部有效位用完了。对这类问题要么换变量让量级归一化要么接受误差并提前设计稳定结构。2.2 复现灾难性消去一元二次方程“标准公式”翻车教科书求根公式 x (-b ± sqrt(b^2 - 4ac)) / 2a 在 b^2 远大于 4ac 时会触发经典数值事故灾难性消去catastrophic cancellation。拿具体方程举例x^2 - 1634.0 x 0.5 0两根相差超过 1600 倍。用两个版本计算小根import numpy as np S 1634.0 # 两根之和 P 0.5 # 两根之积 disc np.sqrt(S * S - 4.0 * P) # 版本1直接套教科书公式求得小根 x2_naive (S - disc) / 2.0 print(x2_naive) # 约 3.0517578125e-05真实小根约 3.06156e-04相对误差超过 90% # 版本2先用加法算大根再用韦达定理反推小根 x1 (S disc) / 2.0 x2_stable P / x1 print(x2_stable) # 约 3.06159e-04误差回到正常水平逻辑说明版本1里 S 和 disc 都约等于 1634两者相减只保留最高位的几个有效数字其余全是舍入噪声小根相对误差被放大若干个数量级。版本2中加法没有消去效应先算大根再用两根乘积等于 P 的性质反推小根避开减法。参数说明当 disc 与 S 的相对差小于约 sqrt(eps)double 下约 1e-8时直接相减就会开始丢有效位。实际代码里不要迷信“标准公式”先看这两个数靠得多近越近越要找稳定等价形式。2.3 用条件数量化“问题有多敏感”numpy 一行拿到 cond(A)精度误差还有一个更系统化的度量条件数。它衡量问题本身不是算法对输入扰动的敏感程度。对 Ax bcond(A) 是误差放大的上界参数解 x 的相对误差约等于 cond(A) 乘上输入数据的相对误差。动手看一眼import numpy as np A np.array([[1.0, 1.0], [1.0, 1.0 1e-8]]) b np.array([2.0, 2.0 1e-8]) print(np.linalg.cond(A)) # 约 4e8严重病态 x np.linalg.solve(A, b) print(x)逻辑说明这个矩阵只比“两行相等”偏离了 1e-8但两行几乎线性相关cond(A) 高达 4e8。在 double 下cond(A) × eps 约等于 1e-8相当于解的有效位数快被耗尽结果自然不稳定。参数说明我的判断习惯是cond(A) 1 / eps 时大多数情况可用cond(A) 1 / sqrt(eps) 时已经危险任何输入扰动都会被放大到影响有效位这时候要么换变量归一化要么改走更稳定的分解。误差来源主要有两条路一条是浮点运算本身的舍入另一条是问题本身对扰动的高敏感度。先分清这两条路各占多少再谈稳定算法否则后面所有参数都是盲调。3. 稳定性怎么量化前向误差、后向误差与扰动检验3.1 前向误差与后向误差误差到底跑到哪里去了稳定性回答一个具体问题误差在算法推进过程中会被放大多少倍。数值误差理论把它分成两句话。前向误差forward error是近似解和真实解的差比如 ||x_calc - x_true|| / ||x_true||。这个指标最直观但真实解通常根本不知道否则就不需要算了。后向误差backward error换个角度近似解 x_calc 不一定精确但能不能找到一个尽量小的输入扰动让“扰动过的方程”精确解正好是 x_calc换句话说误差是被算进了解里还是被算进了矩阵数据里。后向稳定的算法有一个重要性质它产生的解相当于对“稍微改过的原始输入”求精确解只要输入本身也是测量值这个扰动通常可以接受。算法稳定性看后向误差问题本身的敏感性看条件数前向误差 ≈ 条件数 × 后向误差。这句话是整套检验方法的核心实际工程里不需要推导只需要把两边分别算出来对比。代码可以这样贴近import numpy as np A np.array([[2.0, 1.0], [1.0, 3.0]]) b np.array([1.0, 2.0]) x_calc np.linalg.solve(A, b) # 相对残差后向误差的工程近似 res np.linalg.norm(b - A x_calc) / (np.linalg.norm(A) * np.linalg.norm(x_calc) 1e-300) print(relative residual:, res) # 一般在 1e-16 级别的量级逻辑说明相对残差是后向误差最容易拿到的上界估计。如果 res 在 1e-15 附近说明算法本身没引入显著额外误差如果 res 在 1e-8 甚至更高算法已经在帮输入数据“加戏”了。参数说明分母里加 1e-300 只是防除零的保底不要用绝对残差做判断因为绝对残差会随矩阵量级一起膨胀掩盖真实问题。3.2 残差检查不够用跟残差配合的扰动测试残差小真的代表解对吗对病态矩阵答案是否定的残差可以小到 1e-14解却可能错到 10%。原因是病态矩阵把数据扰动放大但扰动方向恰好在零空间附近时残差依然很小。我一般会常备一个扰动测试把输入数据打乱一点看解的变化幅度是否和 cond(A) 对得上import numpy as np def stability_probe(A, b, pert_scale1e-8, trials50, seed42): rng np.random.default_rng(seed) x0 np.linalg.solve(A, b) rel_delta [] for _ in range(trials): dA A * (1.0 pert_scale * rng.standard_normal(A.shape)) db b * (1.0 pert_scale * rng.standard_normal(b.shape)) x1 np.linalg.solve(dA, db) rel_delta.append(np.linalg.norm(x1 - x0) / np.linalg.norm(x0)) return np.percentile(rel_delta, 95) A np.array([[1.0, 1.0], [1.0, 1.0 1e-8]]) b np.array([2.0, 2.0 1e-8]) print(stability_probe(A, b))逻辑说明95 百分位扰动幅度如果和 cond(A) × pert_scale 同一个量级说明算法本身没有额外放大误差问题源头是条件数如果这个数字远大于该量级说明算法实现里存在额外的失稳点需要回头查消元、约分或求和顺序。参数说明pert_scale 取 1e-8 是折中太小会被舍入噪声盖住太大会掩盖小扰动行为。trials 取 50 是为了拿到统计分布而非单个偶然结果固定 seed 方便在回归测试里复现。3.3 扰动幅度怎么设把稳定性测试做成自动验收扰动测试不是一次性调试工具。我建议把它做成独立函数写进测试套件任何矩阵求解、拟合、积分模块在提交前都跑一遍。验收标准可以这样定如果 cond(A) 小于 1e8则 95 百分位的解扰动应小于扰动幅度的 100 倍如果 cond(A) 大于 1e12直接挂红牌要求先做归一化或换算法。def assert_stable(A, b, tol1e-4, pert_scale1e-8): cond np.linalg.cond(A) p95 stability_probe(A, b, pert_scale) if cond 1e12: raise RuntimeError(matrix too ill-conditioned, cond {:.2e}.format(cond)) if p95 / max(pert_scale, 1e-15) tol: raise RuntimeError(p95 perturbation too large: {:.2e}.format(p95)) return cond, p95逻辑说明这个函数的作用是把“这个算法稳不稳”从一个玄学问题变成一行输出。矩阵条件数、扰动下的解的波动幅度、残差三者一起看才能区分问题出在数据上还是出在算法实现上。参数说明tol 按项目对解精度的实际需求设定控制类工程我建议从严到 1e-4纯仿真可以先放到 1e-2。跑一次低成本、高回报尤其是后面要接优化器或反馈控制时这个检查能省掉大量半夜排查时间。4. 把不稳定算法改造成稳定版本主元、QR 分解与 SVD 的选型路线4.1 高斯消元不加主元一个能复现的翻车示例高斯消元本身没问题问题出在“不加主元”的偷懒实现。当某一步对角线元素比同一列其他元素小很多时用它做除数会放大舍入误差。import numpy as np def gauss_no_pivot(A, b): A A.copy().astype(np.float64) b b.copy().astype(np.float64) n len(b) for k in range(n - 1): factor A[k1:, k] / A[k, k] A[k1:, k:] - np.outer(factor, A[k, k:]) b[k1:] - factor * b[k] x np.zeros(n) for i in range(n - 1, -1, -1): x[i] (b[i] - A[i, i1:] x[i1:]) / A[i, i] return x A np.array([[1e-10, 1.0], [1.0, 1.0]]) b np.array([1.0, 2.0]) print(gauss_no_pivot(A, b)) # 结果明显偏离 [1, 1] print(np.linalg.solve(A, b)) # numpy 默认做部分主元结果接近 [1, 1]逻辑说明第一行第一列元素是 1e-10用它做第一次消元的除数会让第二行运算整体放大 1e10 倍舍入误差跟着放大。numpy 的solve默认带部分主元所以结果正常真实项目里如果自己写了简化版消元或复用了老代码很容易把这类问题带进来。参数说明判断消元是否安全的快速办法是看每步主元与所在行最大元素的比值比值小于 1e-12 就要敲警钟。更稳妥的做法是让主元始终取当前列绝对值最大的元素这就是部分主元策略。4.2 三种分解的稳定性对比LU、QR、SVD同一个病态矩阵LU、QR、SVD 给出的解稳定性不同。SVD 最稳因为它把所有方向按奇异值显式切成独立通道直接绕开了“用近奇异矩阵当除数”的问题。QR 不需要显式求法方程也比普通 LU 稳代价是计算量高一些。LU 最快在条件数可控时足够用。用一段代码把三种算法放在同一张矩阵上比较import numpy as np A np.array([[1.0, 1.0], [1.0, 1.0 1e-10]]) b np.array([2.0, 2.0 1e-10]) x_lu np.linalg.solve(A, b) # numpy 对 double 标准路径带部分主元 Q, R np.linalg.qr(A) x_qr np.linalg.solve(R, Q.T b) # QR 分解后回代 U, s, Vt np.linalg.svd(A) x_svd Vt.T (U.T b / s) # 直接按奇异值求逆 print(x_lu) print(x_qr) print(x_svd) print(np.linalg.cond(A))逻辑说明这个矩阵的 cond 已超过 1e10明显跨过危险线。LU 在这种场合会丢掉大部分有效位QR 因为用正交变换保持数值稳定效果稍好SVD 则把近零奇异值直接暴露出来改动或剔除它后就不会强行放大噪声。参数说明Q.T b是把右端项投影到 Q 的空间再回代 R x Q^T b。SVD 里U.T b / s是把对角系统拆开解后续再加截断阈值即可。4.3 病态程度三档决定选型参数表与阈值选型别靠感觉按条件数分三档我一般这样定条件数范围推荐算法关键参数与注意点相对代价cond 1e8LU 部分主元默认选择注意输入归一化低1e8 ~ 1e12QR 或带截断的 SVD精度优先时tol max(m, n) * eps * max(s)中cond 1e12先归一化或换变量再用 SVD手动截断奇异值考虑降阶或换基高谨慎U, s, Vt np.linalg.svd(A) tol max(A.shape) * np.finfo(float).eps * s.max() # 常见截断阈值 keep s tol x Vt[:keep.sum(), :].T ((U[:, :keep.sum()].T b) / s[keep])逻辑说明截断阈值把小于 max(m, n) * eps * max(s) 的奇异值视为噪声方向避免用近奇异分量去放大输入扰动。注意先归一化数据再把截断做在归一化后的矩阵上否则阈值没有意义。参数说明max(m, n) 是矩阵行数列数的较大值eps 是浮点机器精度max(s) 是最大奇异值。这个阈值是很多数值库 pinv 的默认取法实际项目中如果奇异值谱本身有明显断点也可以用“断点右侧全部丢弃”的策略比固定阈值更符合业务语义。5. 精度稳定性避坑5 个真实翻车现场与成因5.1 浮点相等判断翻车 vs np.isclose现象两个计算结果理论上相等代码里用if x y:判断结果时而 true 时而 false甚至对同一个输入两次运行结果都不同被当成“玄学”。原因浮点运算不满足数学上的精确相等1/3 在二进制里是无限循环小数舍入后和人类预期的值天然有差异。解决统一改用近似比较import numpy as np expect (0.1 0.2) if np.isclose(expect, 0.3, rtol1e-9, atol1e-12): print(pass)参数说明rtol 是相对容差atol 是绝对容差两者同时生效。数据量级接近 0 时靠 atol 兜底数据量级大时靠 rtol 兜底。不要只给 atol否则数值量级一变测试全红。5.2 大数加小数接龙求和误差悄悄累积成常数现象循环里反复执行s v单步误差看起来很小一百万个周期后总误差大到连整数部分都不对收敛曲线尾部飘出明显毛刺。原因浮点加法执行时会把两个操作数对齐到统一指数小数部分的低位数被截掉。大量小量级项被逐步喂进一个很大的和中每一步都在丢低位。解决优先用按量级合并的求和方式比如math.fsum或者把正负项分组后再加。Kahan 补偿求和也是一种手段import math vals [1e-8] * 1_000_000 [1.0] print(sum(vals) - 1.0) # 输出不是 0约 1.0109509e-08 print(math.fsum(vals) - 1.0) # 输出接近 0逻辑说明math.fsum内部维护一个补偿寄存器把每步截掉的低位收集起来融进下一轮整体误差能控制在几倍 eps 以内。这个场景在标定、积分、滤波里反复出现值得作为默认工具。参数说明如果必须手写高性能累计器就按块排序再求和让量级接近的数先加。绝对量级相差超过 1e10 的数据在求和前先重新归一化。5.3 残差小却解不对把条件数纳入验收标准现象线性方程组结果用于控制器参数标定时残差只有 1e-12但闭环仿真却完全不对把矩阵打印出来各行肉眼看起来也差不多。原因矩阵病态严重残差被矩阵方向性掩盖解的误差仍可能很大。此时 cond(A) 远大于 1 / eps。解决验收工具里强制打印条件数并把条件数与残差一起存档cond np.linalg.cond(A) res np.linalg.norm(b - A x) / (np.linalg.norm(A) * np.linalg.norm(x)) if cond * max(res, np.finfo(float).eps) 1e-8: raise RuntimeError(unreliable solve: cond too high)逻辑说明判断继续用“残差小 正确”会漏掉病态问题补上 cond 后能直接暴露“数据本身不可靠”的情况。条件数是问题属性不是算法的锅。5.4 单精度双精度混用接口层没守住现象内部计算在 double 下正常换成 GPU 推理或 C 接口后精度断崖式变差有人把锅甩给硬件但同样代码在 CPU double 下仍正常。原因接口或中间存储层误用了 float32或者把 double 数据传入只接受np.float32的函数后悄悄做了截断。浮点位数每减一半有效数字直接降到 7 位条件数稍高就全军覆没。解决在接口处强制做类型检查并对跨底层传递的数据显式声明类型def solve_double(A, b): if A.dtype ! np.float64: A A.astype(np.float64) if b.dtype ! np.float64: b b.astype(np.float64) return np.linalg.solve(A, b)逻辑说明类型转换应该发生在传给求解器之前而不是之后。显式检查能防止“局部变量是 double传参后被自动压成 float”这类黑匣子行为。参数说明单精度并非不能用前提是 cond(A) * eps 小于你可接受的误差。模型推理时常用 float16 提速但一旦回传数值参与控制决策必须恢复 double 并重新跑一遍扰动测试。5.5 只测标准矩阵真实数据直接失稳现象测试套件用单位矩阵、对角占优矩阵、随机正态矩阵全部通过上线后遇到真实采集数据线性求解立刻爆炸或者迭代不收敛。原因标准矩阵的条件数普遍适中真实数据往往存在强相关性、不同物理量纲混合导致条件数达到 1e12 甚至更高。测试集没覆盖这类坏形状。解决在测试数据里加入三类特殊矩阵接近成比例的列、量纲相差巨大的行、带噪声测量值。同时把真实数据直接合入回归测试固定 seed 后跑扰动测试作为回归底线。6. 精度、稳定性与性能怎么权衡我的默认检查清单数值算法的精度与稳定性不是二选一需要在性能约束下综合取舍。我自己的默认检查清单如下每个新工程开始前都会过一遍第一先测条件数再选算法。np.linalg.cond一行就能拿到结果条件数在 1e8 以内直接走 LU超过就上 QR 或 SVD。先算条件数再调参比等到残差异常再回头查要快很多。第二所有验收都附带扰动测试。把输入矩阵加一个 1e-8 量级的随机扰动看解的变化幅度是否被限制在合理倍数内。这是整套稳定性判断里最可靠的一个实际动作。第三在关键接口显式标注精度类型。函数参数直接声明np.float64跨语言、跨硬件边界时增加 dtype 断言避免单双精度混用引发的失稳。性能敏感时可以在验证模块用 float32 提前模拟但最终控制参数一律走 double。第四用相对残差而非绝对残差验收。同时把条件数写入日志存档把条件数和残差当作基本观测点长期监控结果变化。对稳定性的养成我在自己代码里保留这条习惯每次求解前先打印cond求解后再打印relative residual两者一旦出现量级异常立即隔离问题源。第五优先改数据表示而不是调大容差。真实项目中很多看似“精度不足”的问题根源是数据量纲和尺度混乱。先把变量归一化到 1 附近再去折腾更复杂的算法往往事半功倍。有一次我把一个病态矩阵前先做列缩放把最大奇异值和最小奇异值之比从 1e14 压到 1e6结果原先要截断到只剩两个奇异值的系统直接能用了。这件事让我确信精度稳定性问题的第一反应不该是换更贵的算法而是先检查输入数据的量纲、存储位数和条件数。这几个动作成本极低、收益却能覆盖大部分实际事故。希望对你也有帮助。本文还有配套的精品资源点击获取