MMA-Sim: Bit-Accurate Reference Model of Tensor Cores and Matrix Cores

发表时间: 2025-12

作者/机构: Peichen Xie, Yang Wang, Fan Yang, Mao Yang Microsoft Research

速读

一句话结论 本文针对现代 GPU 中矩阵乘法加速器因算术规范未公开导致的数值不精确与不一致问题，提出了首个覆盖 10 种 GPU 架构的比特级精确参考模型 MMA-Sim，成功实现了与真实硬件在比特级别上完全等价的算术行为模拟。

要解决什么问题 现代 GPU 广泛集成了矩阵乘法加速器（如 NVIDIA Tensor Cores 和 AMD Matrix Cores）来执行核心的矩阵乘加操作。然而，硬件厂商从未公开这些加速器在浮点矩阵乘法上的精确算术规范。这种黑盒状态导致了两个致命卡点：第一，不同架构在内部计算时存在精度和动态范围的暗中阉割，例如某些架构的 8 位浮点指令累积精度不足，或者在处理 16 位浮点时将次正规数直接冲刷为零，这会直接导致大语言模型在训练时出现不稳定的梯度截断或精度断崖式下降；第二，不同厂商甚至同一厂商的不同代际产品，对同一个矩阵乘加操作采用了完全不同的求和顺序、舍入模式和特殊值处理规则，导致相同的输入在不同硬件上产生不一致的数值输出。这种跨硬件的数值不一致性彻底破坏了深度神经网络训练和推理的可复现性，使得软件开发者在遇到数值异常时无从排查，也让新硬件的设计缺乏对齐标准。

怎么做的 核心思路是将封闭的硬件加速器视为黑盒，通过构造极具针对性的特定输入和随机输入，反向探测并逆向工程出硬件底层的算术行为，进而构建出 9 种白盒化的算术算法。首先，作者通过一致性检验，将复杂的矩阵乘加操作简化为对单一输出元素的点积加法建模，其核心计算式为： $$d = c + \sum_{k=0}^{K-1} a_k b_k$$ 为了探测这个公式在硬件中的具体执行细节，方法设计了四个关键的探测部件。第一是求和顺序探测，通过构造包含极大值、极小值及其相反数的输入，利用浮点数相加时的吸收效应，推导出硬件是采用二叉树成对求和还是多路融合求和。第二是累积精度探测，通过输入 $1.0+\epsilon$ 或 $-1.0+1.0+\epsilon$ 并不断将 $\epsilon$ 减半，直到硬件输出发生跳变，从而精准测出中间加法操作保留的小数位数。第三是舍入模式探测，利用最低有效位构造处于舍入边界的数值，观察硬件是向零舍入、向下舍入还是向最近偶数舍入。第四是特殊值探测，专门输入次正规数、负零和无穷大，观察硬件是否会发生冲刷为零或异常溢出。基于这些探测结果，作者为不同架构拼装出了具体的算术算法。以覆盖最广的融合点积加法（FDA）算法为例，它精确还原了大部分 Tensor Cores 的行为：首先计算无精度损失的乘积 $p_k = s_k \times 2^{e_k}$；接着将所有乘积和累加数 $c$ 对齐到最大指数 $e_{max}$，并将移出的位直接截断保留 $F$ 位小数；然后在这个定点数状态下完成无损求和；最后再将结果归一化并按照特定规则转换回目标浮点格式。通过这种积木式的算法重组，该模型能够精确模拟从输入到输出的每一次位移和舍入。

效果如何 实验通过对比模拟器输出与真实 GPU 硬件输出进行验证。测试硬件阵列极其庞大，涵盖了 NVIDIA 的 7 代架构（从 Volta 架构的 V100 一直到最新的 Blackwell 架构 B200）以及 AMD 的 2 代架构（CDNA2 的 MI250X 和 CDNA3 的 MI300X）。对比基线直接设定为这些真实物理硬件的实际运行路线。在针对每条指令输入超过 100 万组由随机比特流生成的极端和常规测试用例后，量化结果显示，除了部分指令在输出非规范 NaN 时的负载编码存在差异外，该模型在所有正常数值、次正规数和常规异常值上，均实现了与真实硬件 100% 的比特级完全一致。利用这个高精度的模型，作者还测出并证实了多个严重影响实际工程的硬件缺陷。例如，在 Hopper 和 Ada Lovelace 架构上，8 位浮点指令的底层累积精度实际上只有 13 个小数位，这直接解释了此前业界反馈的训练不稳定问题（该问题在 Blackwell 架构中被提升至 25 位才得以修复）。同时，测试发现 AMD CDNA3 架构在执行 TF32、FP16 和 BF16 指令时，采用了一种极其罕见的非对称向下舍入模式，这在特定计算下（如大数值相减后加上一个微小负数）会产生极其夸张的数值偏差，例如理论值为 -0.000001 的计算结果会被硬件直接输出为 -0.25。方法的局限性在于，目前硬件生成 NaN 时的内部负载编码规则依然是黑盒，模型尚无法做到 NaN 负载的完全对齐，但考虑到当前深度学习软件栈并不依赖 NaN 负载传递信息，这一代价在实际应用中完全可以接受。

A1 主要贡献

本文针对现代GPU中矩阵乘法加速器（MMA，如NVIDIA Tensor Cores和AMD Matrix Cores）因其浮点矩阵乘法运算规范不明确且未公开，导致数值不精确和不一致，进而影响深度神经网络（DNN）训练和推理的稳定性和可复现性的问题，提出了MMA-Sim，这是首个揭示了十种GPU架构（八种NVIDIA和两种AMD）中MMA详细算术行为的比特级精确参考模型。

核心问题与研究目标
深度神经网络（DNN）日益增长的计算需求推动了专用浮点矩阵乘法加速器（MMA）的集成。然而，这些高性能MMA的快速发展也带来了新的挑战。由于缺乏公开的浮点矩阵乘法精确算术规范，不同的MMA在处理相同的矩阵乘法运算时，可能因精度和动态范围不足，导致DNN训练准确性下降，尤其是在大语言模型（LLM）中。同时，不同供应商甚至同一供应商不同代际产品之间的算术实现差异，导致了数值结果的不一致，损害了DNN训练和推理的可复现性。因此，本文的目标是构建一个比特级精确的参考模型MMA-Sim，以揭示这些MMA的详细算术行为，从而帮助软件开发者诊断和预防数值问题，并为硬件设计者提供设计数值一致性新MMA的参考。

创新点与主要贡献
本文通过一种基于测试的方法，结合精心设计的特定输入和随机输入，系统地剖析了MMA的算术行为，并在此基础上构建了MMA-Sim模型。主要贡献如下：

A3 背景知识/关键Observation/设计原则

MMA的计算密集性与数值问题。浮点矩阵乘法是现代深度神经网络（DNN），尤其是大语言模型（LLM）中计算最密集的操作。为了加速这一过程，现代GPU集成了专门的矩阵乘法加速器（MMA），如NVIDIA Tensor Cores和AMD Matrix Cores。这些MMA提供了多种矩阵乘法指令（MMI），用于加速一个或多个矩阵乘加操作：

$D_{M \times N} = A_{M \times K} \times B_{K \times N} + C_{M \times N},$

其中A、B、C是输入矩阵，D是输出矩阵。C和D作为累加器，通过迭代执行矩阵乘加操作来计算大型矩阵的乘法。例如，计算两个64×64矩阵的乘法需要进行(64/M)×(64/N)×(64/K)次矩阵乘加操作。

MMI复杂性与多样性带来的数值挑战。然而，MMI的复杂性和多样性引入了对DNN有深远影响的数值问题。目前没有规范来标准化MMI的设计和算术行为，因此硬件供应商可以自行决定实现方式。这导致不同供应商提供的MMI在形状和数据类型上各不相同（详见附录）。尽管同一供应商的几代MMA可能保留相同的形状和数据类型，但已有研究【[7] M. Fasi, N. J. Higham, M. Mikaitis, and S. Pranesh, “Numerical behavior of NVIDIA tensor cores,” PeerJ Computer Science, vol. 7, p. e330, 2021. [Online]. Available: https://doi.org/10.7717/peerj-cs.330; [8] X. Li, A. Li, B. Fang, K. Swirydowicz, I. Laguna, and G. Gopalakrishnan, “FTTN: Feature-Targeted Testing for Numerical Properties of NVIDIA & AMD Matrix Accelerators,” in International Symposium on Cluster, Cloud and Internet Computing (CCGRID). IEEE, 2024, pp. 39–46. [Online]. Available: https://doi.org/10.1109/CCGrid59990.2024.00014; [9] P. Xie, Y. Gao, Y. Wang, and J. Xue, “Revealing Floating-Point Accumulation Orders in Software/Hardware Implementations,” in USENIX Annual Technical Conference (USENIX ATC), 2025, pp. 1425–1440. [Online]. Available: https://www.usenix.org/conference/ atc25/presentation/xie】揭示，不同代的MMA对同一矩阵乘加操作（公式1）采用不同的计算方法。

未公开实现细节的影响。由于没有规范定义公式1的标准输出，硬件供应商可以自由地用不同的计算精度、动态范围、舍入模式、特殊值处理规则等来实现它。这些实现细节通常是未公开的，但它们对矩阵乘法的数值准确性和更高级别的DNN至关重要。例如，先前研究【[4] DeepSeek-AI, A. Liu, B. Feng, B. Xue, B. Wang, B. Wu, C. Lu, C. Zhao, C. Deng, C. Zhang, C. Ruan, D. Dai, D. Guo, D. Yang, D. Chen, D. Ji, E. Li, F. Lin, F. Dai, F. Luo, G. Hao, G. Chen, G. Li, H. Zhang, H. Bao, H. Xu, H. Wang, H. Zhang, H. Ding, H. Xin, H. Gao, H. Li, H. Qu, J. L. Cai, J. Liang, J. Guo, J. Ni, J. Li, J. Wang, J. Chen, J. Chen, J. Yuan, J. Qiu, J. Li, J. Song, K. Dong, K. Hu, K. Gao, K. Guan, K. Huang, K. Yu, L. Wang, L. Zhang, L. Xu, L. Xia, L. Zhao, L. Wang, L. Zhang, M. Li, M. Wang, M. Zhang, M. Zhang, M. Tang, M. Li, N. Tian, P. Huang, P. Wang, P. Zhang, Q. Wang, Q. Zhu, Q. Chen, Q. Du, R. J. Chen, R. L. Jin, R. Ge, R. Zhang, R. Pan, R. Wang, R. Xu, R. Zhang, R. Chen, S. S. Li, S. Lu, S. Zhou, S. Chen, S. Wu, S. Ye, S. Ma, S. Wang, S. Zhou, S. Yu, S. Zhou, S. Pan, T. Wang, T. Yun, T. Pei, T. Sun, W. L. Xiao, and W. Zeng, “DeepSeek-V3 Technical Report,” 2024, arXiv: 2412.19437. [Online]. Available: https://doi.org/10.48550/arXiv.2412.19437; [5] PyTorch Developers, “PyTorch Developer Notes - Numerical accuracy,” 2025. [Online]. Available: https://docs.pytorch.org/docs/stable/notes/ numerical accuracy.html】报道，NVIDIA Hopper架构上8位浮点（FP8）MMI的累积精度不足，以及CDNA2架构上16位浮点（FP16和BF16）MMI的动态范围不足，这些都可能导致DNN训练不稳定和训练精度显著下降。由于根本原因是硬件未公开的内部算术行为，软件开发者很难诊断和解决这些问题。

跨硬件的数值不一致性。此外，由于MMI在形状和算术行为上存在差异，不同的MMA对相同的浮点矩阵乘法会产生不一致的输出【[6] X. Li, A. Li, B. Fang, K. Swirydowicz, I. Laguna, and G. Gopalakrishnan, “Discovery of Floating-Point Differences Between NVIDIA and AMD GPUs,” in International Symposium on Cluster, Cloud and Internet Computing (CCGRID), 2024, pp. 663–666. [Online]. Available: https://doi.org/10.1109/CCGrid59990.2024.00083】。这损害了DNN的可复现性，使得在不同MMA上保持一致的DNN行为变得非常困难。不了解现有MMA的算术行为，新的硬件供应商也很难确保与它们的数值一致性 。

构建比特级精确参考模型的必要性。因此，完整而准确地建模MMA的算术行为对于解决这些问题至关重要。为此，本文旨在构建一个MMA的比特级精确参考模型，该模型能提供对每个MMI算术行为的彻底剖析和比特级精确模拟。

A2 方法细节

A. 工作流程概述

比特级精确参考模型的定义。比特级精确参考模型是一个可执行的行为模型，其算术行为与相应硬件在比特层面等价。对于每个矩阵乘法加速器（MMA）和每个矩阵乘法指令（MMI），该模型提供一个算术算法$f$，使得对于由MMI指定的形状和数据类型的任何输入矩阵A、B和C，都满足：

$f(A, B, C) = \text{MMI}_{\text{MMA}}(A, B, C),$

基于测试的建模方法。为了构建这个比特级精确参考模型，我们采用一种基于测试的非侵入式方法，将硬件行为视为一个黑盒。具体来说，对于每一代Tensor Cores和Matrix Cores，我们检查每个MMI，并遵循以下工作流程为该MMI构建一个算术算法$f$：

  1. 首先，我们构造精心设计的输入来测试MMI，并根据其在MMA上的输出来检测MMI的特性。这些特性包括求和顺序、累积精度、舍入模式、特殊值处理规则等。
  2. 基于检测到的特性，我们构造一个算法来计算公式1，并用Python实现该算法。
  3. 然后，我们通过对比测试来验证算法的正确性。给定超过一百万组随机输入，如果我们的算法输出与硬件输出在比特级别上一致，我们就在此停止，并检查下一个MMI。
  4. 否则，我们分析输出的差异并相应地优化我们的算法以修复不匹配。测试和优化的过程会重复进行，直到验证成功。

在本节的其余部分，我们首先展示了MMI的矩阵乘加操作建模可以简化为点积加法操作的建模。然后，我们介绍我们测试的主要特性以及我们如何为它们设计特殊输入。

B. 问题简化

输出元素计算的一致性检验。我们首先检验输出矩阵D的每个元素是否都以相同的方式计算——即，点积加法计算

$$d_{ij} = c_{ij} + \sum_{k=0}^{K-1} a_{ik}b_{kj}$$

是否依赖于索引i和j。为了测试这一点，我们首先设置i=0和j=0，并为k=0, 1, ..., K-1生成随机输入$a_{0k}$、$b_{k0}$和$c_{00}$。然后，我们为所有其他的i和j设置$a_{ik} = a_{0k}$，$b_{kj} = b_{0j}$，以及$c_{ij} = c_{00}$。通过将该输入发送到MMA，我们看到输出中的每个元素$d_{ij}$在比特级别上都是相同的。我们在大规模随机测试中观察到了这种一致的行为。因此，D中的每个输出元素都是以相同的方式计算的。

简化为点积加法算法。在本文的后续部分，我们将公式3的表示法简化为

$$d = c + \sum_{k=0}^{K-1} a_k b_k$$

$$= f(c, a_0, a_1, \ldots, a_{K-1}, b_0, b_1, \ldots, b_{K-1}).$$

因此，我们的目标被简化为为每个MMI的点积加法计算构建一个(2K+1)输入的算法$f(c, a, b)$。

C. 求和顺序

使用特殊输入检测求和顺序。我们使用精心设计的输入来检测公式4计算中的求和顺序。输入的构造方式是使得K+1个加数

$$\begin{aligned} p_k = \begin{cases} a_k b_k & \text{if } 0 \leq k < K \\ c & \text{if } k = K \end{cases} \end{aligned}$$

由一个非常大的2的幂$X$、其负数$-X$以及K-1个值为非常小的2的幂$y$的相同数字组成。当将这样的输入发送到MMA时，硬件输出除以$y$的结果表示在$p_i = X$和$p_j = -X$中和（$p_i+p_j=0$）之后相加的加数数量，这样我们就可以确定$p_i$和$p_j$之间求和操作的优先级。我们在不同位置放置$X$和$-X$来运行这个数值实验，从而可以根据优先级信息构建求和树。我们利用【[9] P. Xie, Y. Gao, Y. Wang, and J. Xue, “Revealing Floating-Point Accumulation Orders in Software/Hardware Implementations,” in USENIX Annual Technical Conference (USENIX ATC), 2025, pp. 1425–1440. [Online]. Available: https://www.usenix.org/conference/ atc25/presentation/xie】中的树构建算法来构造求和树。

典型的求和树结构。通过这种方法，我们可以生成一个代表求和顺序的求和树。例如，图1展示了我们为所有代际的Tensor Cores和Matrix Cores上的MMI构建的四种典型求和树。前两种是二叉树，表示求和由不同顺序的二元加法操作构成。图1(a)展示了Ampere架构上DMMA.884指令的顺序求和树。图1(b)展示了CDNA2架构上v mfma f32 32x32x8 f16指令的求和树，表示在四元组内进行成对求和，然后在c和各组和之间进行顺序求和。后两种是多路树，t元树表示求和基于t项融合求和，即t个项在中间求和结果被归一化为浮点数之前相加。图1(c)展示了Ampere架构上HMMA.1684.F32.TF32指令的求和树，表示使用了5项融合求和。图1(d)展示了Ampere架构上HMMA.1688.F32.TF32指令的求和树，表示迭代使用了两个五项融合求和。

图1: Tensor Cores和Matrix Cores上的四种典型求和顺序。#i代表乘积a_i*b_i。
图1: Tensor Cores和Matrix Cores上的四种典型求和顺序。#i代表乘积a_i*b_i。

D. 累积精度

二元加法精度的检测。我们使用精心设计的输入来检测公式4求和中的精度。对于每个基于二元加法运算的求和，我们使用$1.0 + \epsilon$来测试其精度。具体来说，对于$p_i$和$p_j$的求和，我们首先设置$p_i = 1.0$和$p_j = \epsilon = 1.0$，并将其他加数设为零。接下来，我们反复将$\epsilon$减半，直到硬件输出满足$d \neq 1.0 + \epsilon$。然后，我们得出结论，该操作中$p_j$的累积精度在小数点后为$1 - \log_2 \epsilon$位。通过对不同的i和j索引运行此数值实验，我们可以确定所有加法操作中所有加数的精度。

t项融合求和精度的检测。对于每个t项融合求和，我们使用$-1.0 + 1.0 + \epsilon$来测试操作的精度。类似地，对于包含$p_i, p_j, p_k$的求和操作，我们首先设置$p_i = -1.0, p_j = 1.0, p_k = \epsilon = 1.0$，并将其他加数设为零。接下来，我们不断将$\epsilon$减半，直到硬件输出$d \neq \epsilon$。然后，我们得出结论，$p_k$的精度在小数点后为$1 - \log_2 \epsilon$位。通过对不同的i, j, k索引运行此数值实验，我们可以确定所有求和操作中所有加数的精度。

E. 舍入模式

舍入模式的检测。在检测到累积精度后，我们通过精心设计的输入来分析每个求和操作的舍入模式。我们首先设置最低有效位的单位$u = 2\epsilon$。然后，我们构造输入，使未舍入的和分别等于$1.0 + 0.75u, 1.0 + 0.25u, -1.0 - 0.75u$和$-1.0 - 0.25u$。根据舍入后的结果，我们参考表I来识别舍入模式，可能是向上舍入（RU）、向下舍入（RD）、向零舍入（RZ）、远离零舍入（RA）或向最近舍入（RN）。

向最近舍入的平局规则检测。如果舍入模式属于RN，我们进一步测试$1.0+0.5u, 1.0+1.5u, -1.0-0.5u$和$-1.0-1.5u$。根据舍入后的结果，我们参考表II来识别RN的平局决胜规则，可能是舍入平局向上（RNU）、舍入平局向下（RND）、舍入平局向零（RNZ）、舍入平局远离零（RNA）、舍入平局到偶数（RNE）或舍入平局到奇数（RNO）。

表I：根据四个输入的舍入结果检测舍入模式。
图片描述

表II：根据四个输入的舍入结果检测向最近舍入模式的平局决胜规则。
图片描述

F. 特殊值处理

特殊值的定义。次正规数、负零、无穷大和NaN是浮点算术中的特殊值。我们构造边界条件下的输入来检测MMI的特殊值处理规则。

次正规数的处理。对于次正规数，我们首先将c设置为$0.5\epsilon$（其中$\epsilon$是最小的正规数），并将其他输入设为零。如果硬件输出为零，则表明c的次正规输入被冲刷为零。类似地，我们将一个输入$a_k$（或$b_k$）设为$0.5\epsilon$，将其对应的$b_k$（或$a_k$）设为1，并将其他输入设为零。如果硬件输出为零，则表明$a_k$（或$b_k$）的次正规输入被冲刷为零。接下来，我们通过设置输入$a_k = 0.5$，其对应的$b_k = \epsilon$，并将其他输入设为零，来检测次正规值是否在乘法后被冲刷。最后，我们通过设置一个加数$p_i = 1.5\epsilon$，另一个加数$p_j = -\epsilon$，并将其他输入设为零，来检测次正规值是否在加法后被冲刷。

负零的处理。对于负零，我们通过将输入设置为$-0.0 + \sum_{k=0}^{K-1} (+0.0) \times (-0.0)$来检测是否可以生成负零。硬件输出是否等于-0.0表明是否可以生成负零。

无穷大的处理。对于无穷大，我们通过为某个k设置$a_k = 2.0$，设置$b_k = 2^E, c = -2^E$，并将其他输入设为零（其中E是输入浮点格式的最大指数）来检测乘法是否会溢出。我们还通过测试$2^E + 2^E - 2^E$来检测中间求和是否会溢出。然后，我们通过测试边界条件$(2-u)\times2^E +0.75u\times2^E, (2-u)\times2^E +0.5u\times2^E, (2-u)\times2^E +0.25u\times2^E, (-2+u)\times2^E -0.75u\times2^E, (-2+u)\times2^E -0.5u\times2^E$和$(-2+u)\times2^E -0.25u\times2^E$来检测溢出规则。

NaNs的处理。对于NaNs，我们为某个k设置$p_k = 0.0 \times \infty$，并将其他输入设为零；然后我们为某i和j设置$p_i = -\infty$和$p_j = +\infty$，并将其他输入设为零。我们发现在所有测试的MMA中，这两种情况都可以生成NaNs。

比特级精确参考模型

九种算术算法的构建。基于我们对十个GPU架构（包括所有八个带Tensor Cores的NVIDIA GPU架构，从Volta到RTX Blackwell，以及最新的两个带Matrix Cores的AMD架构，CDNA2和CDNA3）的矩阵乘法指令（MMI）进行的数值实验，我们总共构建了九种算术算法。表III列出了我们为Tensor Core指令构建的算术算法，表IV列出了为Matrix Core指令构建的算术算法。

MMA-Sim的实现。基于这些算术算法，我们构建了一个名为MMA-Sim的模拟器，它作为我们的比特级精确参考模型。使用MMA-Sim，用户可以通过提供输入浮点矩阵A、B和C来模拟任何GPU架构的任何MMI，并获得输出矩阵D，使得$d_{ij} = f(c_{i,j}, a_{i0}, a_{i1}, ..., b_{0j}, b_{1j}, ...)$，其中f是该MMI的算术算法。

A. 顺序融合乘加（SFMA）算法

SFMA算法描述。顺序融合乘加（SFMA）算法是为Tensor Core的DMMA指令以及Matrix Core的双精度（FP64）和单精度（FP32）MFMA指令构建的。在此算法中，$f(c, a, b)$通过顺序应用标准的融合乘加（FMA）操作来计算，如算法1所示。

算法1 顺序融合乘加（SFMA）算法
算法1 顺序融合乘加（SFMA）算法

B. 分组成对求和（GPS）算法

GPS算法描述。分组成对求和（GPS）算法是为CDNA2架构上的Matrix Core 16位浮点（FP16和BF16）MFMA指令构建的。在此算法中，$f(c, a, b)$是标准浮点加法和乘法操作与次正规数冲刷的组合。具体来说，计算包括以下四个步骤：

  1. 冲刷输入次正规数：如果A、B或C中的任何数字是次正规数，它将被冲刷为正零（+0.0）。
  2. 乘法：每个$a_k$与$b_k$使用标准的FP32乘法操作相乘，生成FP32乘积$p_k$。然后，如果$p_k$是次正规数，它将被冲刷为零，并保留其符号。
  3. 成对求和：每G个连续的乘积使用标准的FP32加法操作成对求和，并进行符号保留的次正规数冲刷，如算法2所示。此步骤生成K/G个FP32组和，表示为$\gamma_0, \gamma_1, ..., \gamma_{K/G-1}$。
  4. 顺序求和：最后，c和K/G个组和使用标准的FP32加法操作进行顺序求和，并进行符号保留的次正规数冲刷。

GPS算法实现。算法3详细说明了GPS算法。不同的指令使用不同的G值，如下所示：

算法2 成对求和函数

要求：$p_0, p_1, ..., p_{G-1}$
确保：$\gamma = \sum p_k$

function PAIRWISESUM(p_0, p_1, ..., p_{n-1})
    if n = 1 then
        return p_0
    else
        γ ← PAIRWISESUM(p_0, p_1, ..., p_{n/2-1})
        γ ← γ + PAIRWISESUM(p_{n/2}, p_{n/2+1}, ..., p_{n-1})
        γ ← γ × 0.0 if |γ| < 2⁻¹²⁶ ▷ 符号保留地冲刷为零
        return γ

算法3 分组成对求和（GPS）算法

图片描述
图片描述

C. 融合点积加法（FDA）算法

FDA算法描述。融合点积加法（FDA）算法是为Tensor Cores上所有的HMMA、QMMA、HGMMA、HQMMA、UTCHMMA和UTCQMMA指令构建的，涵盖了从TF32到FP4的浮点格式。在此算法中，一个FDA操作的计算包括五个步骤：

  1. 检查特殊值：即无穷大和NaN。如果输入中有任何NaN，或者存在$0 \times \infty$的情况，输出d为NaN（FP32输出编码为0x7FFFFFFF，FP16输出为0x7FFF，即NVIDIA的规范NaNs）。如果乘积和c中同时存在+∞和-∞，输出为NaN。如果只存在一种无穷大，输出就等于它。否则，不存在特殊值，进入下一步。
  2. 计算乘积：使用$a_k$和$b_k$的有效数和指数计算乘积。具体来说，我们用$s_{ak} \times 2^{e_{ak}}$和$s_{b} \times 2^{e_{bk}}$表示$a_k$和$b_k$，其中对于正规数$1 \le |s| < 2$，对于次正规数$|s| < 1$，指数从它们的原始浮点表示中获得。乘积$p_k = a_k b_k = s_k \times 2^{e_k}$通过$s_k = s_{ak} \times s_{bk}$和$e_k = e_{ak} + e_{bk}$计算。此步骤中乘积不进行归一化，有效数以定点数形式计算和存储，没有精度损失、溢出或下溢。
  3. 对齐与截断：将乘积和c对齐到它们的最大指数$e_{max} = \max\{e_c, e_0, e_1, e_2, ...\}$。每个有效数$s_k$（包括$s_c, s_0, s_1, ...$）相应地右移$e_{max} - e_k$位。结果只有F个小数位，移出的位被截断（实际上是向零舍入）。F的值取决于架构和指令：
    • Volta: HMMA为F=23。
    • Turing和Ampere: HMMA为F=24。
    • Ada Lovelace: HMMA为F=24，QMMA为F=13。
    • Hopper: HMMA和HGMMA为F=25，QGMMA为F=13。
    • Blackwell: HMMA、UTCHMMA和UTCQMMA为F=25。
    • RTX Blackwell: HMMA和QMMA为F=25。
  1. 求和：由于在第3步中进行了对齐和截断，此步骤只需使用定点算术（F个小数位）对移位后的有效数求和。因此，求和结果是精确的，且与求和顺序无关。
  2. 归一化：将和归一化为输出的浮点格式，即FP32或FP16。第4步后的和表示为$s_{sum} \times 2^{e_{max}}$。对于FP32输出，该数被归一化到FP32的指数范围。如果归一化结果的绝对值大于或等于$2^{128}$，则变为无穷大。否则，其有效数在第13个小数位（Ada Lovelace QMMA或Hopper QGMMA指令）或第23个小数位（其他指令）处向零舍入（RZ）。对于FP16输出，该数被归一化到FP16的指数范围，然后其有效数在第10个小数位使用向最近偶数舍入（RNE）规则进行舍入。如果舍入结果的绝对值大于或等于$(2-2^{-11}) \times 2^{15}$，则变为无穷大。否则，再次遵循RNE规则将其归一化为FP16。

TF32截断。虽然TF32数编码为E8M10，只有19个有效位，但它存储在32位中。我们发现Tensor Cores总是将13个最低有效位视为零进行计算。这导致一个有趣的后果：NaN值可能被转换为无穷大。例如，如果输入元素是0x7F800001（NaN），它将被转换为无穷大（0x7F800000）。

比例因子。QMMA.SF和UTCQMMA指令为MXFP8、MXFP6和MXFP4格式引入了UE8M0比例因子。为了乘以这些比例因子，只需修改第1步和第2步。在第1步，如果任何比例因子是NaN，结果应为NaN。在第2步，乘积$p_k = a_kb_k = s_k \times 2^{e_k}$通过$s_k = s_{ak} \times s_{bk}$和$e_k = e_{ak} + e_{bk} + e_{aSF\lfloor k/32 \rfloor} + e_{bSF\lfloor k/32 \rfloor}$计算。

表III：从Volta到RTX Blackwell的八个架构上NVIDIA Tensor Core指令的比特级精确算术算法。每个单元格显示了精确模拟该指令的算法。N/A表示该指令在该架构上不可用。
图片描述

表IV：CDNA2和CDNA3架构上AMD Matrix Core指令的比特级精确算术算法。指令前缀“v mfma”已省略。除非特别说明，我们使用CDNA3上重新格式化的指令名称。
图片描述

D. FDA链（CoFDA）算法

CoFDA算法描述。对于大多数Tensor Core指令，$d = g(c, a_0, a_1, ..., a_{K-1}, b_0, b_1, ..., b_{K-1})$，其中g(c, a, b)代表一个FDA操作。然而，对于Ampere和Ada Lovelace上的HMMA.1688.F32.TF32、HMMA.16816.F32、HMMA.16816.F16和HMMA.16816.F32.BF16，以及Ada Lovelace上的QMMA.16832.F32.f8type.f8type和QMMA.16832.F16.f8type.f8type，一个FDA操作只能计算K/2对被乘数，d是使用两个FDA操作链式计算的，即：

$d=g\left(g\left(c, a_0, a_1, \ldots, a_{K/2-1}, b_0, b_1, \ldots, b_{K/2-1}\right), \right. \\ \left. a_{K/2}, a_{K/2+1}, \ldots, a_{K-1}, b_{K/2}, b_{K/2+1}, \ldots, b_{K-1}\right).$

E. 分组点积融合求和（GDFS）算法

GDFS算法描述。分组点积融合求和（GDFS）算法是为Tensor Cores上的所有OMMA和UTCOMMA指令构建的，涵盖了MXFP4和NVFP4数据类型。此算法中，一个GDFS操作的计算包括七个步骤：

  1. 检查特殊值：由于FP4数没有无穷大或NaN，UE4M3和UE8M0数没有无穷大，此步骤相对简单。如果任何比例因子是NaN或c是NaN，输出d为NaN（编码为0x7FFFFFFF）。如果c是无穷大，d=c。否则，进入下一步。
  2. 计算乘积：使用定点算术计算$a_k$和$b_k$的乘积$p_k=a_kb_k$。此步骤无精度损失。
  3. 分组求和：每16个连续的乘积使用定点算术求和，生成4个组点积，表示为$\sigma_0, \sigma_1, \sigma_2, \sigma_3$。
  4. 乘以比例因子：将组点积乘以各自的比例因子。结果表示为有效数和指数$\gamma_k = \sigma_k a_{SF\lfloor 16k/S \rfloor} b_{SF\lfloor 16k/S \rfloor} = s_k \times 2^{e_k}$。有效数$s_k$以F=35个小数位的定点数存储，到目前为止没有精度损失。
  5. 对齐与截断：将$s_0, s_1, s_2, s_3$和c对齐到它们的最大指数$e_{max}$，并使用向零舍入（RZ）规则截断移出F=35个小数位的比特。
  6. 求和：使用定点算术对第5步的结果求和。
  7. 归一化：将和归一化为FP32。

UE4M3截断。虽然UE4M3数只有7个有效位，但它存储在8位中。Tensor Cores总是将其最高有效位视为零。

F. 融合点积向下舍入加法（FDRDA）算法

FDRDA算法描述。融合点积向下舍入加法（FDRDA）算法是为CDNA3架构上的Matrix Core TF32、FP16和BF16 MFMA指令构建的。此算法中，一个FDRDA操作的计算包括八个步骤：

  1. 检查特殊值：与FDA算法的第1步相同。
  2. 计算乘积：与FDA算法的第2步基本相同，但乘积可能溢出。如果$p_k \ge 2^{128}$，溢出为+∞；如果$p_k \le -2^{128}$，溢出为-∞。
  3. 再次检查特殊值：如果乘积中同时存在+∞和-∞，输出d为NaN。如果只有一种无穷大，输出等于它。否则，进入下一步。
  4. 对齐乘积：将乘积对齐到它们的最大指数$e_{dot} = \max\{e_0, e_1, e_2, ...\}$，截断移出的位。此步骤与FDA算法的第3步类似，但不包括c。此算法中，F=24。
  5. 求和乘积：使用定点算术对乘积求和，不包括c。结果表示为$s_{dot} \times 2^{e_{dot}}$。
  6. 对齐点积结果与c：将点积结果与c对齐到它们的最大指数$e_{max} = \max\{e_{dot}, e_c\}$。有效数$s_{dot}$右移$e_{max}-e_{dot}$位，移出第31个小数位的部分被向下舍入。有效数$s_c$右移$e_{max}-e_c$位，移出第24个小数位的部分被向下舍入。
  7. 求和：使用定点算术对两个舍入后的有效数求和，生成$s_{sum}$。
  8. 归一化：使用向最近偶数舍入（RNE）规则将结果$s_{sum} \times 2^{e_{max}}$归一化为FP32。

TF32截断。与Tensor Cores类似，Matrix Cores总是将TF32输入的13个最低有效位视为零，这可能将NaN转换为无穷大。

G. FDRDA链（CoFDRDA）算法

CoFDRDA算法描述。对于CDNA3上的大多数TF32、FP16和BF16 MFMA指令，$d = g(c, a_0, ..., a_{K-1}, b_0, ..., b_{K-1})$，其中g(c, a, b)代表一个FDRDA操作。然而，对于v_mfma_f32_16x16x8_xf32、v_mfma_f32_16x16x16_f16和v_mfma_f32_16x16x16_bf16，一个FDRDA操作只能计算K/2对被乘数，d是使用两个FDRDA操作链式计算的，其形式如公式6所示。

H. 分组融合点积向下舍入加法（GFDRDA）算法

GFDRDA算法描述。分组融合点积向下舍入加法（GFDRDA）算法是为CDNA3架构上的Matrix Core FP8 MFMA指令构建的。此算法中，一个GFDRDA操作的计算包括十个步骤：
1.、2.、3. 特殊值检查与乘积计算：与FDRDA算法的第1、2、3步相同。

  1. 分组对齐乘积：将偶数索引的乘积对齐到它们的最大指数$e_{even} = \max{e_0, e_2, ...}$，将奇数索引的乘积对齐到它们的最大指数$e_{odd} = \max{e_1, e_3, ...}$。移出F=24个小数位的比特被向零舍入（RZ）截断。
  2. 分组求和：分别使用定点算术对两组乘积求和，结果表示为$s_{even} \times 2^{e_{even}}$和$s_{odd} \times 2^{e_{odd}}$。
  3. 对齐组和：将两个点积对齐到它们的最大指数$e_{dot} = \max{e_{even}, e_{odd}}$。移出F=24个小数位的比特被向下舍入（RD）。
  4. 求和：使用定点算术对两个舍入后的有效数求和，生成$s_{dot}$。
  5. 对齐点积结果与c：将点积结果与c对齐到它们的最大指数$e_{max} = \max{e_{dot}, e_c}$。$s_{dot}$在第31个小数位向下舍入。$s_c$在第24个小数位舍入，但有一个额外检查：如果$e_c < e_{max} - 25$，则$s_c$向零舍入，否则向下舍入。

9.、10. 最终求和与归一化：与FDRDA算法的第7、8步相同。

I. GFDRDA链（CoGFDRDA）算法

CoGFDRDA算法描述。一个GFDRDA操作只能计算16对被乘数。因此，对于K=32的FP8 MFMA指令，d是使用两个GFDRDA操作链式计算的，其形式如公式6所示。

A4 实验环境

软件配置

硬件配置
实验在以下GPU上运行，涵盖了所有配备Tensor Cores的NVIDIA架构以及最新的两种配备Matrix Cores的AMD架构：

A4 实验结果

正确性验证

主要发现与洞察
通过MMA-Sim，本文确认了一些已知问题并发现了新的未公开特性：

  1. 累积精度的降低与提升

    • 发现: 先前报道的NVIDIA Hopper架构FP8指令累积精度不足导致LLM训练不稳定的问题【[4] DeepSeek-V3 Technical Report, 2024, arXiv】，在Ada Lovelace架构上也同样存在。两种架构的QMMA/QGMMA指令均使用13个小数位的累积精度，而非[4]中报道的14位尾数位。
    • 发现: NVIDIA在其最新的Blackwell和RTX Blackwell架构中，已将FP8矩阵乘法指令的累积精度提升至25个小数位，与TF32/FP16/BF16指令的精度相同，从而解决了此问题。AMD CDNA3的FP8指令也使用了24个小数位。
    • 架构建议: 未来MMA设计应至少使用23个小数位（即FP32的尾数位数）的累积精度。
    • 软件建议: 在Hopper和Ada Lovelace架构上，建议开发者使用FP16或BF16指令进行FP8矩阵乘法，以牺牲效率换取更高精度。

  2. 有限的动态范围

    • 发现: AMD CDNA2架构中将次正规数冲刷为零的行为会影响DNN训练稳定性【[5] PyTorch Developer Notes - Numerical accuracy, 2025】，此问题在CDNA3架构中通过完全支持次正规数得到解决。
    • 发现: CDNA3上的中间乘积在绝对值超过$2^{128}$时仍会溢出为无穷大，而NVIDIA Tensor Cores则不会。
    • 架构建议: 未来MMA设计应完全支持FP32的动态范围，包括次正规数，并避免中间结果溢出。
    • 软件建议: 在CDNA2上，建议对FP16矩阵乘法的输入进行缩放或归一化以避免次正规数，或如[5]所建议，使用BF16进行计算。

  3. CDNA3上的非对称舍入

    • 发现: CDNA3架构的TF32、FP16和BF16 MFMA指令使用了非对称的向下舍入（RD）模式，而大多数其他MMA使用对称的向零舍入（RZ）或向最近偶数舍入（RNE）。
    • 影响: 这种非对称舍入在特定条件下（如A、B矩阵数值大，C矩阵数值小且为负）会引入不可忽略的数值误差和偏差。例如，计算2048*2048 - 2048*2048 - 0.000001的结果为-0.25，远偏离真实值。图2的分布图显示了RD模式导致的负向偏差，而假设的RZ模式则表现出对称性。
    • 发现: AMD在CDNA3的FP8 MFMA指令中增加了一个补丁，当累加数c的指数远小于最大指数时，会将其向零舍入而非向下舍入，从而缓解了上述问题。
    • 架构建议: 建议未来的MMA设计使用与浮点数编码一致的符号-数值算术，并选择RZ作为最自然的舍入模式。
    • 软件建议: 在CDNA3上，建议仅使用Matrix Cores进行矩阵乘法（设置C=0），并使用通用计算单元进行全FP32精度的矩阵累加，以减轻TF32/FP16/BF16指令的此问题。
      图2: 使用向下舍入（RD）模式的CDNA3 FP16 MFMA指令的数值偏差δRD，以及使用假设的向零舍入（RZ）模式的FP16 MFMA指令的数值偏差δRZ的分布。

A7 补充细节

现有MMA算术行为建模研究。虽然许多研究尝试对现有MMA的算术行为进行建模，但没有一个被验证为比特级精确。Raihan等人【[15] Modeling Deep Learning Accelerator Enabled GPUs, 2019, ISPASS】直观地将Tensor Core的点积操作建模为乘法和加法操作的完整二叉树，但Hickmann等人【[16] Experimental Analysis of Matrix Multiplication Functional Units, 2019, ARITH】通过数值实验证明这是不正确的，并指出了Volta架构Tensor Core的几个特性。Fasi等人【[7] Numerical behavior of NVIDIA tensor cores, 2021, PeerJ Computer Science】将数值实验扩展到Turing和Ampere架构的Tensor Cores，Li等人【[8] FTTN: Feature-Targeted Testing for Numerical Properties of NVIDIA & AMD Matrix Accelerators, 2024, CCGRID】进一步扩展到Hopper架构的Tensor Core和CDNA1/CDNA2架构的Matrix Cores。Xie等人【[9] Revealing Floating-Point Accumulation Orders in Software/Hardware Implementations, 2025, USENIX ATC】发明了一种算法来揭示软硬件实现中的浮点累积顺序。然而，这些研究只覆盖了加速器算术行为的一部分，不足以构成一个比特级精确的参考模型。此外，Valpey等人【[17] An SMT Formalization of Mixed-Precision Matrix Multiplication: Modeling Three Generations of Tensor Cores, 2025, arXiv】报告称其中一些结果是不正确的。由于缺乏比特级的验证，这些研究的结果可靠性不足。

MMA导致的AI研发挑战。由于缺乏标准和一致的算术行为，矩阵乘法加速器（MMA）给AI研究和开发带来了挑战。Deepseek开发者【[4] DeepSeek-V3 Technical Report, 2024, arXiv】报告称，由于Hopper架构FP8矩阵乘法指令的累积精度不足，他们的LLM训练精度显著下降。为缓解此问题，他们使用额外的CUDA Cores进行全FP32精度的矩阵累加。PyTorch开发者【[5] PyTorch Developer Notes - Numerical accuracy, 2025】报告称，CDNA2架构的FP16和BF16 MFMA指令会将次正规数冲刷为零，而次正规值在反向传播的梯度计算中频繁出现，可能导致FP16 DNN训练无法收敛。为缓解此问题，他们在反向传播中将FP16转换为BF16进行计算，这增加了动态范围但牺牲了精度。

数值不一致性与可复现性问题。此外，先前的研究报告了由于浮点算术的复杂特性，在不同硬件【[10] Estimation of numerical ´ reproducibility on CPU and GPU, 2015, FedCSIS】、GPU供应商【[11] Testing GPU Numerics: Finding Numerical Differences Between NVIDIA and AMD GPUs, 2024, SC24-W】以及不同代MMA【[6] Discovery of Floating-Point Differences Between NVIDIA and AMD GPUs, 2024, CCGRID】之间存在数值不一致问题。DNN的可复现性是AI研究人员的一个重要目标【[12] Problems and Opportunities in Training Deep Learning Software Systems: An Analysis of Variance, 2020, ASE; [13] Towards Training Reproducible Deep Learning Models, 2022, ICSE; [14] Defeating Nondeterminism in LLM Inference, 2025, Thinking Machines Lab】，但与硬件相关的数值不一致性使其变得相当困难【[9] Revealing Floating-Point Accumulation Orders in Software/Hardware Implementations, 2025, USENIX ATC】。

A5 结论

本文介绍了MMA-Sim，这是首个揭示NVIDIA和AMD GPU中矩阵乘法加速器（MMA）详细算术行为的比特级精确参考模型。通过系统性的测试和大规模验证，MMA-Sim实现了与真实硬件的比特级等价，并揭示了许多先前未公开的算术行为。

MMA-Sim将作为开源项目发布，为透明的数值分析、跨供应商比较、算术行为复现以及未来加速器和AI系统的精度感知协同设计提供一个公开可用的工具。

A6 附录

A. 用于矩阵乘法的浮点格式

传统与低精度格式。浮点矩阵乘法是数值计算、高性能计算和深度学习中最基本和主导的操作之一。传统上，IEEE-754双精度（FP64）和单精度（FP32）格式最为流行，其矩阵乘法遵循标准的IEEE-754算术【[18] IEEE Standard for Floating-Point Arithmetic, 2019】。随着深度神经网络（DNN）的普及，尤其是在大语言模型（LLM）中，为了实现高效率和可扩展性，低精度浮点格式被广泛采用。

各类低精度格式。

B. Tensor Core指令

指令族与命名规则。NVIDIA GPU提供Warp级的矩阵乘加指令（DMMA, HMMA, QMMA, OMMA）以使用Tensor Cores执行矩阵乘加操作。此外，Hopper架构提供Warp组级指令（HGMMA, QGMMA），Blackwell架构提供Warp组级/Warp组对级指令（UTCHMMA, UTCQMMA, UTCOMMA）。这些指令及其兼容性列于表III。指令名称通常指示数据类型和形状（MNK或MxNxK），例如“16816”表示M=16, N=8, K=16。

各指令族的数据类型。

$$d_{ij} = c_{ij} + \sum_{k=0}^{K-1} a_{ik}a_{i[k/S]}^{SF} b_{kj}b_{[k/S]j}^{SF}$$


其中S=32是块大小。

C. Matrix Core指令

指令族与命名规则。近期的AMD GPU提供Warp级的矩阵融合乘加（MFMA）指令以使用Matrix Cores执行矩阵乘加操作。本文覆盖了CDNA2和CDNA3架构，其Matrix Core指令列于表IV。