浮点精度与数值稳定性：IEEE 754、FP32/FP64/BF16/TF32 与混合精度

引言

浮点数是科学计算的底座，也是性能与正确性交锋最激烈的地方。把 FP64 换成 FP32 往往能直接省掉一半内存带宽，把 TF32 或 BF16 交给张量核则能带来数倍吞吐，但代价是误差悄悄进入结果——很多「加速成功」的案例最后是在收敛性变差或结果对不上参考文献时才暴露问题。

本文按「表示 → 舍入 → 误差 → 稳定性 → 类型版图 → 混合精度 → 张量核 → 验证 → 迁移」讲解浮点精度与数值稳定性：IEEE 754 的位布局与舍入语义、灾难性消去与条件数、补偿求和与迭代精化、FP32/FP64/TF32/BF16/FP16/FP8 的取舍边界，以及一套可落地的精度验证与可复现性方法。

前置：/hpc-roofline-model/（Roofline 判定计算还是访存瓶颈）、/hpc-gpu-kernel-optimization/（GPU 内核与张量核优化）、/hpc-ai-hpc-convergence/（AI 与 HPC 的混合负载）。

目录

• 1. IEEE 754 基础：位布局与表示范围
• 2. 舍入模式与特殊值
• 3. 误差来源与条件数
• 4. 数值稳定性：求和与补偿算法
• 5. 数据类型版图：从 FP64 到 FP8
• 6. 混合精度计算：迭代精化与残差校正
• 7. 张量核与混合精度实践
• 8. 精度验证与可复现性
• 9. 实战：从 FP64 到混合精度的迁移
• 10. 速查表与一句话记忆
• 延伸阅读

1. IEEE 754 基础：位布局与表示范围

1.1 三种字段

一个浮点数由三部分组成：符号位、指数位（决定动态范围）、尾数位（决定精度）。指数采用偏置表示，尾数隐含一个前导 1（正规数）：

值 = (-1)^sign × 1.mantissa × 2^(exponent - bias)

FP64: 1 + 11 + 52 位   bias = 1023
FP32: 1 +  8 + 23 位   bias = 127
FP16: 1 +  5 + 10 位   bias = 15
BF16: 1 +  8 +  7 位   bias = 127

1.2 关键指标对比

注意 BF16 与 FP32 的指数位相同，因此动态范围一致、只是精度低；而 FP16 的指数位只有 5 位，动态范围窄得多，这也是为什么深度学习里 BF16 比 FP16 更不容易上溢。

1.3 机器 epsilon 的含义

机器 epsilon 是 1.0 与下一个可表示浮点数之间的距离。它决定了任何单次运算的相对误差下界：一次加减乘除的结果相对误差不超过 u（单位舍入误差）。这个数字是所有误差分析的基本单位。

2. 舍入模式与特殊值

2.1 四种舍入模式

IEEE 754 规定了四种舍入方向，默认是「就近舍入、逢半取偶」（round to nearest, ties to even）：

RN  就近取偶（默认）        误差最小，可复现
RZ  向零舍入                截断，误差有偏
RU  向正无穷舍入            区间算术的上界
RD  向负无穷舍入            区间算术的下界

编译器默认假设 RN 且不读取舍入模式寄存器。若代码里用 fesetround 改变舍入方向，必须加 -frounding-math（GCC）或 -ffp-model=strict（Clang），否则常量折叠会算出错误结果。

2.2 特殊值

±0        有符号零，0.0 == -0.0 为真但 1/0.0 与 1/-0.0 符号不同
±Inf      溢出或除以零的结果，参与运算仍为 Inf
NaN       无效运算结果（0/0、Inf-Inf），与任何值比较均为假
Subnormal 非正规数，指数全零，用于表示比最小正规数更小的值

非正规数对性能有实际影响：部分 GPU 与加速器对非正规数走慢路径甚至直接刷零（flush to zero）。在需要极小量级累加的内核里，这个差异会改变结果。

2.3 FMA 与融合运算

FMA（Fused Multiply-Add）把 a*b + c 用一条指令完成，中间乘积不产生舍入，既更快也更准。但它改变结果的最后几位，导致开启与关闭 FMA 时结果不一致：

# GCC 默认 -ffp-contract=fast，允许跨语句收缩为 FMA
gcc -O3 -ffp-contract=off -o app app.c    # 关闭收缩，便于对拍

3. 误差来源与条件数

3.1 三类误差

表示误差   输入数据无法被浮点精确表示（如 0.1）
舍入误差   每次运算引入不超过 u 的相对误差
截断误差   离散化与迭代截断引入的方法误差

前两者是浮点固有的，第三者属于数值方法本身。混合精度优化只影响前两者，绝不能拿它去掩盖截断误差。

3.2 条件数与误差放大

问题的条件数决定了输入扰动被放大多少倍。对线性方程组 Ax = b，条件数 κ(A) = ‖A‖·‖A⁻¹‖：

前向误差 ≤ 条件数 × 后向误差

后向误差很小（算法稳定）但条件数极大（问题病态）时，结果依然不可信。这就是为什么残差小不等于误差小：迭代法里 ‖b - Ax_k‖ 可以很小，而 ‖x_k - x‖ 依然很大，/hpc-petsc-solver/ 中的迭代求解器同样受此约束。

3.3 灾难性消去

两个量级相近的数相减，有效数字会被大幅抵消：

/* 病态写法：b^2 远大于 4ac 时，一个根会发生灾难性消去 */
double x1 = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);

/* 稳定写法：用韦达定理求另一个根，避免相近数相减 */
double x1 = (-b - copysign(sqrt(b * b - 4.0 * a * c), b)) / (2.0 * a);
double x2 = c / (a * x1);

判别式 b*b - 4*a*c 本身也可能因 a、b、c 的表示误差而失去意义，此时问题已经病态，换公式也救不回来。

4. 数值稳定性：求和与补偿算法

4.1 朴素求和的误差增长

把 n 个数依次累加，误差随 n 线性增长，最坏情况相对误差约为 n·u·Σ‖x_i‖ / ‖Σx_i‖。当正负项相消时，分母很小而分子不小，误差被放大几个数量级。

4.2 Kahan 补偿求和

double kahan_sum(const double *x, size_t n) {
    double sum = 0.0, c = 0.0;
    for (size_t i = 0; i < n; ++i) {
        double y = x[i] - c;        /* 减去上一轮丢失的低位 */
        double t = sum + y;
        c = (t - sum) - y;          /* 本轮丢失的低位 */
        sum = t;
    }
    return sum;
}

Kahan 的代价是破坏了结合律，无法直接向量化或并行化。在 GPU 上更常用的是分块成对求和：每个线程块内用 Kahan 或成对求和，块间用树形归约。这也解释了为什么 OpenMP 的 reduction 与手写累加结果不同——归约顺序不同。

4.3 内积与矩阵乘的误差界

内积 x·y 的朴素计算误差约为 n·u·Σ‖x_i·y_i‖。这也是高精度 BLAS 存在的理由：某些库提供补偿内积（如 ddot 的扩展精度版本），代价是吞吐下降数倍。在矩阵乘里，累加精度比乘法精度更关键——用 BF16 输入但 FP32 累加，误差远小于用 BF16 累加。

5. 数据类型版图：从 FP64 到 FP8

5.1 选型的第一原则

选类型时先问两个问题：动态范围够不够（会不会上溢或下溢）与精度够不够（误差能否被问题容忍）。动态范围由指数位决定，精度由尾数位决定，两者互不替代。

FP64 → 科学计算默认：CFD、分子动力学、谱方法
FP32 → 中等精度需求：图像处理、部分机器学习、单精度求解器
TF32 → 张量核上的「FP32 加速档」：范围同 FP32，精度约 10 位尾数
BF16 → 深度学习训练：范围同 FP32，精度低但不易上溢
FP16 → 推理与混合精度训练：精度尚可但范围窄，需 loss scaling
FP8  → 大模型推理与训练：需要 per-tensor 或 per-block scaling

5.2 带宽视角

混合精度最大的直接收益往往不是算力而是带宽。同样规模的数组，FP32 比 FP64 少占一半内存、少传一半数据；在 Roofline 的访存受限区，这意味着接近 2 倍的加速。这也是为什么很多内核在算术强度不变的情况下，仅换数据类型就能大幅提速。

5.3 算力视角

在张量核上，各类型的峰值吞吐差异巨大（以某代数据中心 GPU 为例的典型比例）：

这些比例是数量级参考，具体取决于芯片代际与指令形态。

6. 混合精度计算：迭代精化与残差校正

6.1 核心思想

混合精度不是「把所有 FP64 换成 FP32」，而是把便宜的操作放在低精度、把关键的操作留在高精度。最经典的框架是迭代精化（Iterative Refinement）：

求解 Ax = b：
  1. 用低精度（FP32 或 FP16）做 LU 分解，得到近似解 x_0
  2. 用高精度（FP64）计算残差 r = b - A x_k
  3. 用低精度分解求解 A d = r，更新 x_{k+1} = x_k + d
  4. 重复直到残差满足要求

6.2 为什么有效

分解是 O(n³) 的主要开销，放在低精度能带来数倍加速；残差计算是 O(n²)，留在 FP64 代价很小却能保证最终精度接近 FP64。收敛条件是迭代矩阵的谱半径小于 1，通常几次迭代即可收敛。

/* 伪代码：FP32 分解 + FP64 残差精化 */
lu_factor_fp32(A, LU);                 /* 一次性开销，走 FP32 */
x = lu_solve_fp32(LU, (float *)b);
for (int k = 0; k < max_iter; ++k) {
    r = b - matvec_fp64(A, x);         /* 关键：FP64 残差 */
    if (norm2_fp64(r) < tol) break;
    d = lu_solve_fp32(LU, (float *)r); /* 修正量走 FP32 */
    x = x + d;
}

6.3 更现代的变体

当矩阵条件数较大时，纯精化可能不收敛，此时用 Krylov 方法做修正步（GMRES-IR）能显著扩大适用范围。这类方法把低精度分解当作预条件子，而不是精确求解器，因此对精度损失更宽容。

7. 张量核与混合精度实践

7.1 累加必须留在高精度

张量核的基本模式是「低精度输入、高精度累加」：FP16 或 BF16 的乘法结果在 FP32 累加器里相加。这是混合精度能保持可接受误差的关键——乘法精度影响单次误差，累加精度影响误差随规模的累积。

推荐：  A(fp16) × B(fp16) → 累加器 fp32
避免：  A(fp16) × B(fp16) → 累加器 fp16   （误差随 K 维度线性增长）

7.2 数值范围问题

FP16 的最大值是 65504，梯度或激活值稍大就上溢为 Inf。工程上的对策：

• Loss scaling：训练时把 loss 放大若干倍，反向传播后再缩回，让梯度落在 FP16 可表示范围内。
• Per-tensor 或 per-block scaling：FP8 场景下按张量或分块统计最大值，动态选择缩放因子。
• BF16 优先：范围与 FP32 相同，多数训练场景可直接替换而不需 scaling。

7.3 与 Roofline 的结合

混合精度让计算吞吐提升，但访存带宽往往不变。因此在访存受限的内核里，换精度必须先减少数据搬运（例如数据本身就用低精度存储），否则算力提升会被内存墙吃掉。判断方法仍然是 Roofline 模型。

8. 精度验证与可复现性

8.1 与高精度参考对拍

最可靠的验证方式是构造一个已知解析解或高精度参考解的算例，然后统计误差：

from mpmath import mp
mp.dps = 50                      # 50 位十进制精度
ref = mp.quad(lambda t: mp.exp(-t * t), [0, 1])

对大规模问题无法求解析解时，退而求其次的做法是：用 FP128 或双倍-双倍（double-double）跑一遍，比较低精度结果的相对误差是否落在理论界内。

8.2 编译器带来的非确定性

同一份代码在不同优化级别下结果不同，原因通常是：

FMA 收缩       -ffp-contract=fast 把 a*b+c 融合，少一次舍入
向量化重排     向量归约改变求和顺序
OpenMP 归约    线程数与调度策略改变归约树形状
快速数学       -ffast-math 允许重新结合

若要位级可复现，需要同时固定：-ffp-contract=off、固定线程数与调度、禁用重排类优化，并用确定性归约。

8.3 有界误差而非位相等

更现实的工程目标是有界误差：给定输入，结果与参考解的相对误差不超过某个阈值。把阈值写进 CI，任何优化改动只要越界就失败。这比追求位级一致更可维护。

9. 实战：从 FP64 到混合精度的迁移

9.1 迁移步骤

□ 建立 FP64 基准解与误差度量（相对误差、迭代次数、收敛残差）
□ 识别可降精度的环节：矩阵组装、分解、矩阵向量乘、预处理
□ 先只降分解精度，残差保留 FP64，做迭代精化
□ 逐项开启低精度路径，每步都跑误差回归
□ 记录迭代次数变化，收敛变慢说明精度已到临界
□ 对最终结果做一次 FP64 复算校验

9.2 一个稀疏求解器的实测

以某三维 CFD 隐式求解器为例（单节点，双路 CPU）：

结论很典型：全低精度不可用，混合精度可用且加速明显；精化迭代次数随分解精度下降而增加，这部分开销必须计入总时间。

9.3 常见坑与对策

10. 速查表与一句话记忆

一句话记忆：浮点优化的口诀是「范围看指数位、精度看尾数位，误差看条件数；把便宜的大头（分解、矩阵乘）放到低精度，把关键的残差与累加留在 FP64/FP32，用迭代精化把精度补回来，每一步都要有误差回归与高精度参考对拍」。

延伸阅读

• /hpc-roofline-model/ — 判断计算受限还是访存受限
• /hpc-gpu-kernel-optimization/ — 张量核与内核精度取舍
• /hpc-ai-hpc-convergence/ — AI 负载与 HPC 的融合
• /hpc-amd-rocm/ — AMD GPU 上的混合精度实践
• /hpc-molecular-dynamics/ — 分子动力学中的精度与截断
• 高性能计算专题 — 高性能计算专题