数值算法精度与稳定性实战：从浮点误差到混合精度

简介《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) # Trueb 被完全吞掉 # 场景二累加顺序影响结果 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 条件数问题本身的“病态程度”条件数衡量的是输入扰动对输出的放大倍数。对于线性方程组 Axb条件数 cond(A) ||A||·||A⁻¹||。条件数为 10^k 意味着输入有 k 位有效数字的损失。这是问题固有的属性跟你用什么算法无关。import numpy as np # 构造一个 Hilbert 矩阵经典病态例子 n 12 H np.array([[1.0/(ij1) 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 指数增长n12 时已经接近双精度极限。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 内部用分块和向量化累加误差远小于朴素三重循环。参数上 n2000 时手动循环会跑几十秒实际调试可以降到 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对对称矩阵也能算但数值路径不同可能产生微小虚部。参数上n300 时差异还小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_THREADS1编译时去掉-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, dtypetorch.float32) B torch.randn(n, n, dtypetorch.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 量级。参数上n4096 是 Tensor Core 的典型工作尺寸太小则加速比不明显。迭代精化是另一个技巧解出 x 后算残差 r b - Ax再用 r 解一次修正量更新 x。对于条件数在 1e8 以内的矩阵两到三轮迭代精化能把前向误差降到机器精度。import numpy as np def iterative_refinement(A, b, max_iter3): 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_iter3对多数问题够用1e-14是收敛阈值。注意残差要用高精度算否则修正量本身就不准。我自己的习惯是任何涉及矩阵分解的代码先跑一遍条件数估计再决定用不用迭代精化。条件数超过 1e12 就直接上高精度库别在双精度里硬扛。数值算法这件事后悔药很少前期多花十分钟做误差分析比后期调一周的“玄学 bug”划算得多。希望帮到你。本文还有配套的精品资源点击获取