07.GPU 与 NPU 精度误差分析#
Author: 陈春昀
!!!!!!!!!!!! 文章的开篇介绍，本节要介绍的内容
硬件浮点运算如何引入误差#
浮点表示#
!!!!!!!!!!!! 加入 GPU 和 NPU 当前用的浮点数介绍，BF16、FP16、FP8、FP6、FP4。
十进制中我们使用科学计数法表示一个浮点数，例如
\((-4.5)_{10}\times10^{0} = (-1)\times(4\times10^{0}+5\times10^{-1})\)
同样的数用二进制表示为
\((-1.001)_{2}\times2^{2}=(-1)\times(1\times2^2+0\times2^1+0\times2^0+1\times2^{-1})\)
他们的值都是-4.5，是同一个数。
从上面可以看出，表示一个浮点数，我们只需要三个部分：
符号位表示该数为正数还是负数，如上述例子中的-1；
小数点左边只有一位数的小数部分，如上述例子中十进制的 4.5，或者二进制的 1.001；
指数部分，如上述例子中等号左边十进制的指数 0，或二进制的指数 2。
计算机硬件也是使用这三个部分表示一个浮点数：符号（sign，S），尾数（fraction，F），指数（exponent，E）。
我们对浮点数的表示做一些改进：
二进制下浮点数总是可以写成 \(±1.xxxx\times2^{yyyyy}\) 的形式，最左边的数字总是 1，可以省去，这样可以为硬件节省一位的空间；
如果需要比较两个浮点数的大小，指数部分大的浮点数一定比指数部分小的浮点数大，可以先比较两个浮点数的指数部分。我们把浮点数的指数加上一个 bias 得到指数 E，从而保证 E 为非负数（无符号数），相较于比较有符号数，硬件可以更快地比较两个无符号数的大小。假设有 n 位指数，则 \(bias=2^{n-1}-1\)。
对于一个（S，F，E）组合而成的浮点数，其对应的浮点数的值为：
但上述表示存在一个问题：当 \((E-bias)=0\) 时，正浮点数的值为 \(1.F\)，最小只能表示到 1.0，无法表示 (0,1)之间的数，负浮点数也是一样，无法表示 (-1,0)。这种表示方法也无法表示 0。为了解决这一问题，IEEE 754 标准使用一种非规格化方法，约定当 \((E-bias)=0\) 时，不再假定尾数为 \(1.F\) 的规格化形式，而是 \(0.F\) 的非规格化形式，且指数值被固定为 \(1-bias\) ，即
此外，IEEE 还规定了 0、无穷、NaN 的表示，这里不做重点介绍。
为简化理解，我们使用 1 位符号、 2 位指数（此时 bias=1）、2 位尾数表示一个浮点数，且 S=0（即只考虑非负浮点数）。一共可以表示的数如下表所示：
E |
指数值 |
F |
S=0 时，表示的数字 |
备注 |
|---|---|---|---|---|
00 |
1-bias = 0 |
00 |
0 |
非规格化 |
00 |
1-bias = 0 |
01 |
\((0.01)_2\times2^{0}=1\times2^{-2}=0.25\) |
非规格化 |
00 |
1-bias = 0 |
10 |
\((0.10)_2\times2^{0}=2\times2^{-2}=0.5\) |
非规格化 |
00 |
1-bias = 0 |
11 |
\((0.11)_2\times2^{0}=3\times2^{-2}=0.75\) |
非规格化 |
01 |
E - bias = 0 |
00 |
\((1.00)_2\times2^{0}=1\times2^{0}=1\) |
|
01 |
E - bias = 0 |
01 |
\((1.01)_2\times2^{0}=1\times2^{0}+1\times2^{-2}=1.25\) |
|
01 |
E - bias = 0 |
10 |
\((1.10)_2\times2^{0}=1\times2^{0}+1\times2^{-1}=1.5\) |
|
01 |
E - bias = 0 |
11 |
\((1.11)_2\times2^{0}=1\times2^{0}+1\times2^{-1}+1\times2^{-2}=1.75\) |
|
10 |
E - bias = 1 |
00 |
\((1.00)_2\times2^{1}=1\times2^{1}=2\) |
|
10 |
E - bias = 1 |
01 |
\((1.01)_2\times2^{1}=1\times2^{1}+1\times2^{-1}=2.5\) |
|
10 |
E - bias = 1 |
10 |
\((1.10)_2\times2^{1}=1\times2^{1}+1\times2^{0}=3\) |
|
10 |
E - bias = 1 |
11 |
\((1.11)_2\times2^{1}=1\times2^{1}+1\times2^{0}+1\times2^{-1}=3.5\) |
|
11 |
00 |
特殊值无穷大，∞ |
||
11 |
非 0 |
不是一个数，NaN |
上表一共表示了从 0 到 3.5 之间的 12 个数，如下图所示。
然而，在 0 到 3.5 之间存在无数个实数，我们使用 5 位的浮点数只能准确表示其中的 12 个数。如果要表示其他的浮点数，比如表示 2.25 这个数，怎么办？要么增加更多位尾数，比如增加一位表示尾数 F， F 从 2 位增加到 3 位，F=001，使用 \((1.001)\times2^1\) 表示 2.25；要么就只能进行舍入操作，使用距离 2.25 比较近的两个浮点数之一，即 2 或者 2.5 来表示。总结就是，当一个浮点数需要太多位尾数才能准确表示，就会发生舍入，而舍入会引入精度误差。
IEEE 754 规定了五种舍入规则：
向最接近可表示值舍入 (Round to nearest, ties to even)
向正无穷方向舍入 (Round toward +∞)
向负无穷方向舍入 (Round toward -∞)
向零舍入 (Round toward zero)
向最接近可表示值舍入，中间值远离零方向 (Round to nearest, ties away from zero)
硬件厂商可以选择其中任意一种舍入模式，并没有强制规定必须使用哪一种，其中最常用且默认的是向最接近可表示值舍入模式。如果硬件厂商选择的舍入模式不同，且计算时涉及舍入操作，那么不同产品得到的结果很可能是不同的。
为了量化了浮点数表示的不确定性和舍入误差，我们使用 ULP（Unit in the Last Place）衡量误差，ULP 表示在给定指数下两个相邻可表示浮点数之间的最小间隔。一种更直观的理解是：ULP 是浮点数末尾 1 个 bit（最低有效位）所代表的值。例如在上述例子中， 2.25 的两个相邻可表示浮点数分别为 2 和 2.5 ，则 \(ULP=2.5-2=0.5\)。
一次舍入操作带来的精度误差不会超过 1ULP，似乎并不大，但使用舍入之后的浮点数再进行算术运算，可能会发生误差累计或者放大误差，比较典型的就是除法运算，且除数接近 0，例如用上面的浮点数计算 \(0.25/0.125\) ，精确的结果应该是 2 ，但 0.125 会被舍入到 0.25 或者 0，如果舍入到 0.25，得到的结果是 1，是精确结果的一半；如果舍入到 0，会发生除 0，得到的结果是无穷，误差更大。
浮点运算#
我们以浮点数的加法为例，说明浮点运算操作如何带来误差。为了便于理解，我们使用两个只有 4 位有效位，指数为 2 位的十进制数作为示例：
**步骤一：**指数对齐。对指数较小的数进行调整，使其指数等于较大的数的指数。这里会使用到右移操作，由于只有 4 位有效位，最后一位 1 会被丢弃。
**步骤二：**有效位相加。
**步骤三：**规格化，并检查结果中的指数位是否发生上溢或者下溢，若发生，抛出异常。
**步骤四：**若没有发生上溢或者下溢，对结果进行舍入。
**步骤五：**检查步骤四的结果是否为规格化数，若不是，返回步骤三，若是，得出结果。
从上面的例子可以看出，加法操作有两个地方会引入误差：步骤二的移位操作，步骤四的舍入操作。
移位操作：假设两个浮点数的指数部分相差为 n，则较小的浮点数需要右移 n 位与较大的浮点数指数对齐，较小的浮点数会丧失其最低 n 位的精度。这在两个浮点数的大小相差非常大的时候会很明显。例如：计算 \(a+b+(-a)\) 的值，准确的答案应该是 b，但如果 a 和 b 的指数相差太大，b 会因为移位操作丧失最低 n 位的数值，使 b 变成一个更靠近 0 的数，最后的结果不等于 b。
a = 1000
b = 0.0000000000000125
a+(-a)+b
#输出 1.25e-14
a+b+(-a)
#输出 0.0
假如我们需要计算多个数的和，比如 \(a+b+b+b+...+b\) ，一共 \(10^{14}\) 个 b， \(10^{14}\times b=1.25\) ，正确的结果应该是 1001.25，但由于计算顺序为 \(a+b+b+b+...+b\)，对 b 进行移位后变成了 0，计算机不停地执行 \(a+b=a\) ，最后结果为 1000。
那么应该如何尽可能避免这种误差呢？一种缓解这种情况的操作是，将加法中的元素按照升序排序之后再相加。较小的浮点数累加成较大的浮点数后，再与更大的浮点数相加，可以减少两个浮点数的指数之间的差距，从而减少右移的位数。
在并行计算中，常会将数字分成组，并行计算每一组内的数字的和，类似的，将这些数字进行预先排序，使较为接近的数字分在一组，可以减少移位操作带来的误差。
舍入操作：上一节已经介绍过，当一个浮点数需要太多位尾数才能准确表示，就会发生舍入，而舍入会引入精度误差。这里不再赘述。
上面介绍了加法操作，乘法操作可以转化为加法操作，除法操作和超越函数（如三角函数）则更为复杂，往往通过多项式逼近的算法实现，具体精度与硬件实现方式有关，可以查看硬件供应商是否提供精度保证，例如英伟达对其数学操作的舍入模式和精度保证进行了说明（https://developer.nvidia.com/zh-cn/blog/cuda-math-method-cn/）。
实际应用中的浮点误差#
交换律和结合律#
a = 0.1
b = 0.2
c = 0.3
d = 0.4
e = 0.5
f = 0.6
x = (a+b) + (c+d) + (e+f)
y = (a+d) + (b+e) + (c+f)
#计算的结果为
#x = 2.1
#y = 2.0999999999996
先看上面的例子。从数学上来说，上面两个累加的结果 x 和 y 应该相等才对，那么到底计算机得出的结果中哪个才是准确的呢？
其实两个结果都不是准确的！
计算机是不能准确表示 0.1 这个浮点数的，0.1 在二进制中是一个无限循环小数（0.0001100110011...），因此无法精确表示，必须进行舍入。
而 \(2.1=2+0.1=2^1+0.1\) ，因此二进制也不能准确表示 2.1 这个数，之所以显示 x=2.1，是由于 Python 将浮点数转换为字符串进行打印时，它会进行一些处理以提供“友好”的输出。如果我们打印查看更多更多位数，变量 x 就“原形毕露”了，这就是为什么说这两个结果都不准确：
print("{:.20f}".format(x))
#输出 2.10000000000000008882
print("{:.20f}".format(y))
#输出 2.09999999999999964473
但是，相对而言，x 比 y 更准确，计算 x 的时候，累加的元素是按照升序顺序计算的，涉及的右移操作更少，比如计算 0.1+0.2 时，只需要将 0.1 右移一位就能与 0.2 进行指数对齐，最多带来 1ULP 的误差，而计算 0.1+0.4 时需要将 0.1 右移两位，最多带来 3ULP 的误差。
因此，在进行计算的时候，要特别注意计算的顺序，浮点运算是不满足数学上的交换率和结合率的，同一硬件上不同的计算顺序可能得出不同的结果。
硬件实现差异#
!!!!!!!!!!!! 这里需要再深入介绍技术点，加上配图，要深入去理解技术点。不能比视频浅。
浮点运算单元差异：硬件计算过程中，中间结果的保留位数会影响精度。例如进行 16 位的乘加操作 \(a\times b+c\) ，如果硬件先进行乘法操作 \(a\times b\)，得到一个中间结果 tmp，理论上，tmp 必须是 32 位的才能保证没有精度损失，如果 tmp 只有 16 位，中间结果会经历一次舍入。
然后 tmp 与 c 相加，可能还会经历一次舍入。如果硬件上使用融合乘加（FMA）单元计算 \(a\times b+c\)，由于 FMA 会使用 32 位保留中间结果，可以避免一次中间结果的舍入误差。
数值近似优化：不同硬件对低精度计算（FP16）中的非规格化数的处理可能不同。例如在 Tensor Cores ( Volta 及以后)：矩阵乘法（如 FP16 MMA），输入非规格化数通常被默认置零（ Flush-To-Zero，FTZ），从而引发精度损失，但这是速度优化的一部分，如果想避免，可以手动关闭 FTZ。
浮点数格式：计算中使用不同的浮点数格式，所具备的精度不同。深度学习可能采用 16 位或者更少位的浮点格式（如下图），以进行更快的计算。
下面以 FP16、FP8-E5M2、FP8-E4M3 三种格式为例，说明使用浮点表示 0.1 时产生的误差。
import torch
# 检查设备支持
torch.cuda.is_available()
device = torch.device("cuda")
# 创建变量 a、b、c 分别等于 0.1 (格式为 FP16、FP8-E5M2、FP8-E4M3)
a = torch.tensor(0.1, dtype=torch.float16, device=device)
b = torch.tensor(0.1, dtype=torch.float8_e5m2, device=device)
c = torch.tensor(0.1, dtype=torch.float8_e4m3fn, device=device)
#结果
a.item(),b.item(),c.item()
#(0.0999755859375, 0.09375, 0.1015625)
#与 0.1 相比分别对应约 0.02%，6.25%，1.56%的相对误差
软件与算法实现差异#
!!!!!!!!!!!! 这里需要再深入介绍技术点，加上配图，要深入去理解技术点。不能比视频浅。
数学库的底层实现：NV cuBLAS 与昇腾 CANN 算子库在矩阵乘法的分块（Tiling）、循环展开（Loop Unrolling）、内存访问的实现上存在差异，以及不同库对矩阵乘法的并行拆分策略（GEMM 分块大小）差异，但这些的核心原因依然是引起了浮点数加法计算顺序的改变。在进行浮点运算时，凡是涉及计算顺序的改变都要多加小心，计算结果可能不一致。例如进行卷积运算时，不同的数学库，使用维诺格拉德（Winograd）卷积与使用朴素卷积相比结果很可能不一样。
编译器优化：编译器对指令重排（Instruction Reordering）、融合乘加（FMA）指令使用影响中间结果的精度。指令重排可能会影响浮点数计算的顺序。融合乘加对中间结果的保留相比于一次乘法后再进行一次加法，减少了一次中间结果的舍入，更加精确。
其他差异#
!!!!!!!!!!!! 这里需要再深入介绍技术点，加上配图，要深入去理解技术点。不能比视频浅。
线程调度与计算顺序不同：并行计算中，线程块的执行顺序和线程间的数据同步存在非确定性。例如，对 CUDA 多个线程计算的浮点数结果使用原子操作进行加和，先计算出结果的线程先进行加和运算，由于加法的顺序是不确定的，同样的算法同样的设备多次运算可能得出不同的结果，建议最好使用规约（Reduction）算法。
总结与思考#
!!!!!!!!!!!! 总结成一段话，或者把上面最核心的 2/3 个点一句话列出来。
浮点表示不一定是精确的，浮点运算的结果也不一定是精确的；
浮点运算时要特别留意舍入和移位操作带来的误差；
浮点加法计算顺序的改变可能通过影响移位操作的舍入误差积累的方向进而影响到计算结果；
对于浮点计算，应该关注的是计算误差在一个可接受的范围内，而不是追求结果的一致。