Ozaki Scheme II: A GEMM-oriented emulation of floating-point matrix multiplication using an integer modular technique

• 作者/机构: Katsuhisa Ozaki1, Yuki Uchino2, and Toshiyuki Imamura2

A1 主要贡献

本文提出了一种新颖的、基于通用矩阵乘法（GEMM）的矩阵乘法仿真算法，该算法利用中国剩余定理（Chinese Remainder Theorem）进行高精度计算。

• 核心问题：标准浮点运算在进行矩阵乘法时可能精度不足，而传统的高精度计算库或方法（如Ozaki方案）在利用现代硬件（尤其是低精度矩阵引擎）方面存在性能或效率上的限制。
• 研究目标：开发一种新的矩阵乘法仿真方法，该方法既能利用硬件（如GPU上的INT8张量核心）的高性能，又能提供高精度计算结果，同时比现有方法（如传统的Ozaki方案）更高效。
• 创新点与主要贡献：
  1. 提出Ozaki方案II：这是一种基于中国剩余定理的新型矩阵乘法仿真方法。与依赖于将矩阵分解为多个“切片”（slices）的传统Ozaki方案（本文称为Ozaki方案I）不同，新方法通过模运算将高精度计算分解为多个独立的、低精度（如INT8）的矩阵乘法。
  2. 高效利用低精度硬件：该方法专为利用低精度矩阵引擎（如NVIDIA GPU的INT8张量核心）而设计，从而继承了高度优化的GEMM例程的计算效率。
  3. 灵活的精度控制与更高的计算效率：与传统Ozaki方案相比，新方法显著减少了达到同等精度所需的矩阵乘法次数。用户可以通过调整模数（moduli）的数量来直接控制计算量和精度。
  4. 卓越的性能表现：
     • 在GPU上，使用INT8张量核心仿真FP64运算，在NVIDIA RTX 4090上实现了7.4至9.8 TFLOPS的性能，在NVIDIA GH200上实现了56.6至80.2 TFLOPS的性能，后者的性能甚至超过了其原生的FP64运算性能。
     • 在CPU上，使用FP64运算仿真四倍精度（quadruple-precision）运算，与传统的Ozaki方案相比，实现了高达2.3倍的加速。

下表展示了NVIDIA GPU在不同精度下的浮点和整数性能，凸显了FP16和INT8张量核心的巨大吞吐量优势，这是本研究方法所利用的基础。

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

2.1 符号表示

符号定义。本文定义了若干符号和运算。$F$表示由IEEE计算机协会（2019）定义的二进制浮点数集合。常量$u$定义为单位舍入误差，例如FP64的$u = 2^{-53}$。对于$k \in N$，$Z_k$是一个整数集合，其中$a \in Z_k$意味着$|a| \leq 2^k$。符号$fl(\cdot)$表示使用浮点算术计算的结果，并假设计算中不发生上溢或下溢。对于矩阵$A = (a_{ij})$，符号$|A|$表示其元素为$|a_{ij}|$的矩阵。$gcd(a, b)$代表两个整数$a$和$b$的最大公约数。

对称模运算的定义。模运算$a \pmod m$（其中$a \in Z$，$m \in N$）可以有多种定义方式。本文采用对称模，其中余数$r$满足以下条件：
$$-\frac{m}{2} \leq r \leq \frac{m}{2}.$$
形式上，$r = a \pmod m$ 定义为：
$r = a - m \cdot \left\lfloor \frac{a}{m} + \frac{1}{2} \right\rfloor.$
该定义确保余数是在与$a$模$m$同余的整数中离零最近的那个。对于一个矩阵$A$，表达式$A \pmod m$指的是对矩阵的每个元素应用模运算，从而得到一个与$A$维度相同的新矩阵。

2.2 传统的Ozaki方案

矩阵分解。对于矩阵$A \in F^{p \times q}$和$B \in F^{q \times r}$，目标是获得$AB$的近似值。Ozaki方案【15, Ozaki K, Ogita T, Oishi S and Rump SM (2012) Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications. Numerical Algorithms 59(1): 95–118.】【16, Ozaki K, Ogita T, Oishi S and Rump SM (2013) Generalization of error-free transformation for matrix multiplication and its application. Nonlinear Theory and Its Applications, IEICE 4(1): 2–11.】将矩阵分解为：
$$A = A_1 + A_2 + \cdots + A_{k-1} + \underline{A_k},$$
$$B = B_1 + B_2 + \cdots + B_{\ell-1} + \underline{B_\ell},$$
其中
$A_i \in \mathbb{F}^{p \times q}, \quad B_j \in \mathbb{F}^{q \times r}$
对于$1 \leq i \leq k - 1$ 和 $1 \leq j \leq \ell - 1$，并且
$$ \underline{A}_k := A - \sum_{i=1}^{k-1} A_i \in \mathbb{F}^{p \times q}, \quad \underline{B}_l := B - \sum_{i=1}^{l-1} B_i \in \mathbb{F}^{q \times r}. $$
对于公式(2)，我们称$k$为$A$的切片数，$\ell$为$B$的切片数。这里，我们根据Ozaki等人(2012)的研究设置$k = \ell$，但也可以像Ozaki等人(2013)讨论的那样设置$k \neq \ell$。然后，$AB$被转换为：
$$AB = \sum_{i+j \leq k} A_i B_j + \sum_{i=1}^{k-1} A_i \underline{B_{k+1-i}} + \underline{A_k} B.$$
当$i + j \leq k$时，矩阵乘积$A_iB_j$可以使用浮点运算无舍入误差地计算。值得注意的是，Ozaki等人(2025)也提出了$AB$的另一种形式，但本文不予考虑。

Ozaki方案的构成与特点。当$k = \ell$时，Ozaki方案由以下三部分组成：
- 第1部分：如(2)所示的矩阵分解。
- 第2部分：如(3)所示的$k(k + 1)/2$个矩阵乘积的计算。
- 第3部分：对计算出的矩阵乘积进行归约。
计算成本方面，第1部分为$O(pq) + O(qr)$，第2部分为$k(k + 1)pqr + O(pr)$，第3部分为$O(pr)$。因此，对于足够大的$p, q, r$，第2部分主导了总成本。Ozaki方案的一个优点是，可以在这个计算密集的部分应用优化的BLAS例程。然而，一个缺点是需要大量的内存来存储矩阵，因为它们在(2)中被分解为多个加数。在第2部分中，根据矩阵的结构选择适当的BLAS例程。例如，如果$A$是三角矩阵，$B$是通用矩阵，则使用TRMM。如果$A$和$B$没有特殊结构，则使用GEMM。

多组分格式下的Ozaki方案。即使当$A$和$B$以多组分格式表示时，类似的方法也可以实现高精度的计算结果。例如，假设$A$和$B$以双字格式表示：$A := A_h + A_\ell$ 和 $B := B_h + B_\ell$，其中$A_h, A_\ell \in F^{p \times q}$，$B_h, B_\ell \in F^{q \times r}$，使得：
$$f1(A_h + A_l) = A_h, \quad f1(B_h + B_l) = B_h.$$
同样，我们将$A_h + A_\ell$和$B_h + B_\ell$分解为未求和的浮点矩阵之和，使得：
$$A_h + A_\ell \approx A_1 + A_2 + \cdots + A_{e-1} + \underline{A_e},$$
$$B_h + B_\ell \approx B_1 + B_2 + \cdots + B_{f-1} + \underline{B_f},$$
其中$A_i \in F^{p \times q}$，$B_j \in F^{q \times r}$，并且对于$1 \leq i \leq e - 1$和$1 \leq j \leq f - 1$，$fl(A_iB_j) = A_iB_j$。此外，
$$ \mathbb{F}^{p \times q} \ni \underline{A_e} \approx A_h + A_\ell - \sum_{i=1}^{e-1} A_i, \\ \mathbb{F}^{q \times r} \ni \underline{B_f} \approx B_h + B_\ell - \sum_{i=1}^{f-1} B_i. $$
然后，我们基于(3)得到计算结果。请注意，在这种情况下，第3部分的计算需要高精度计算，例如双字算术。

利用矩阵引擎进行仿真。接下来，我们介绍一种使用近期GPU中可用的矩阵引擎来仿真矩阵乘法的方法，典型例子是NVIDIA张量核心（TCs）。设非奇异对角矩阵$D_i \in F^{p \times p}$和$E_i \in F^{r \times r}$的对角元素是2的幂，目的是实现无误差的对角缩放。Mukunoki等人【8, Mukunoki D, Ozaki K, Ogita T and Imamura T (2020) Dgemm using tensor cores, and its accurate and reproducible versions. In: Sadayappan P, Chamberlain BL, Juckeland G and Ltaief H (eds.) High Performance Computing. Cham: Springer International Publishing, pp. 230–248.】设置：
$$ A \approx D_1^{-1} D_1 A_1+D_2^{-1} D_2 A_2+\cdots+D_k^{-1} D_k A_k, \\ B \approx B_1 E_1 E_1^{-1}+B_2 E_2 E_2^{-1}+\cdots+B_k E_k E_k^{-1}, $$
其中对于所有$1 \leq i \leq k$，$D_iA_i$和$B_iE_i$中的所有元素都可以在FP16中表示。乘积$D_iA_i$和$B_iE_i$的项表现得好像它们的有效数字具有
$$\left\lfloor \frac{24 - \log_2 n}{2} \right\rfloor$$
位的精度。然后，$(D_iA_i)(B_jE_j)$可以使用FP16 TCs无舍入误差地计算。需要注意的是，FP16 TCs可以将结果存储为FP16或FP32，但在这种情况下，计算结果存储在FP32中。然后，通过以下方式获得近似值$\tilde{C}$：
$$\tilde{C} := \text{f1}\left( \sum_{i+j \leq k+1} D_i^{-1} (((D_i A_i)(B_j E_j))) E_j^{-1} \right)$$
这涉及$k(k + 1)/2$次矩阵乘法。在计算矩阵乘法$(D_iA_i)(B_jE_j)$之后，结果从FP32转换为FP64，并在FP64中进行求和。

使用INT8张量核心。Ootomo等人【13, Ootomo H, Ozaki K and Yokota R (2024) Dgemm on integer matrix multiplication unit. The International Journal of High Performance Computing Applications in Press. DOI:10.1177/ 10943420241239588.】类似地使用了(5)的形式，其中所有$i$的$D_iA_i$和$B_iE_i$的所有元素都存储在INT8中。考虑到INT8 TCs操作的结果存储在INT32中，他们利用这一特性，在使用INT8 TCs时，当$q \leq 2^{29}$时，实现了乘积$(D_iA_i)(B_jE_j)$的无误差计算。使用INT8 TCs的第一个优点是，理论上它的速度是FP16 TCs的2倍或4倍，如表1所示。第二个优点是，如果$q < 2^{17}$，分裂矩阵的有效数字保持7位，与$q$无关。如果$q > 2^{10}$，则(6)的数量小于7。

Ozaki方案I及其精度趋势。Uchino等人【19, Uchino Y, Ozaki K and Imamura T (2025) Performance enhancement of the ozaki scheme on integer matrix multiplication unit. The International Journal of High Performance Computing Applications DOI:10.1177/10943420241313064. URL https://doi.org/10.1177/10943420241313064.】讨论了ozIMMU的加速和舍入误差分析。他们通过降低第3部分归约的成本实现了加速。具体来说，他们减少了使用FP64执行的求和计算次数。在本文中，我们将这类算法称为Ozaki方案I。使 用INT8 TCs的Ozaki方案I包括以下三个部分：
- 第1部分：如(5)所示的矩阵分解。
- 第2部分：使用INT8 TCs计算$k(k + 1)/2$个矩阵乘积。
- 第3部分：在FP64中对计算出的矩阵乘积进行归约。
在此，我们介绍使用INT8 TCs的Ozaki方案I的精度趋势。在这种情况下，(5)中所有$i$的$D_iA_i$和$B_iE_i$都由INT8表示。我们使用MATLAB生成了两个矩阵：

其中randn(n)生成一个均值为0、方差为1的正态分布随机数$n \times n$矩阵。表2显示了每个切片所需的矩阵乘法次数。图1显示了Ozaki方案I计算(7)的结果相对于切片数（左）和矩阵乘法次数（右）的最大相对误差。从图1可以看出，Ozaki方案I的精度与切片数成正比提高。当切片数从$k$增加到$k + 1$时，需要额外进行$k + 1$次矩阵乘法。因此，如图1所示，当在对数尺度上绘制时，精度相对于矩阵乘法次数并未呈现线性增长。

2.3 中国剩余定理

定理概述。此处简要介绍中国剩余定理，因为它对所提出的方法至关重要。
令$m_1, m_2, \ldots, m_s$是两两互质的正整数，即对于所有$i \neq j$，$gcd(m_i, m_j) = 1$。对于任何给定的整数$a_1, a_2, \ldots, a_s$和

存在一个唯一的整数$x$同时满足以下同余方程组：

那么，$x$在模$M$的意义下是唯一确定的。解可以明确地构造为：

其中$M_i = M/m_i$，而$y_i$是$M_i$模$m_i$的模乘法逆元，即$M_iy_i \equiv 1 \pmod{m_i}$。

在矩阵乘法中的应用。我们将此技术应用于矩阵乘法。设$A' \in Z^{p \times q}$和$B' \in Z^{q \times r}$。令

那么，我们有

在本文中，我们使用这种方法，即使用中国剩余定理进行直接重构。另一种方法是使用Garner方法。

A2 方法细节

Ozaki方案II的提出。本节中，我们提出了一种基于GEMM的方法，该方法结合了中国剩余定理和无误差矩阵乘法，我们称之为Ozaki方案II。设$A \in F^{p \times q}$和$B \in F^{q \times r}$。我们准备两个对角矩阵$D$和$E$，其对角元素为2的幂次方，我们有

其中

算法流程。我们首先设定$s$，即矩阵乘法的次数。接下来，我们从一个表中选取$m_i \geq 2$（对于$1 \leq i \leq s$），其中$m_1, m_2, \ldots, m_s$两两互质，并根据(8)设定$M$。矩阵$A' \in Z^{p \times q}$和$B' \in Z^{q \times r}$生成如下：

对于$1 \leq t \leq s$。如果我们设置合适的$m_i$，那么我们可以使用BLAS例程无舍入误差地计算$A'_iB'_t$。这一点将在第3.1节和第3.2节中以INT8 TCs和FP64为例进行详细讨论。

结果重构。令$C_i \equiv A'_iB'_i \pmod{m_i}$和$Y = \sum_{i=1}^{s} C_iM_iy_i$。那么，我们有
$C' \equiv A'B' \mod M = Y \mod M.$
因此，$c'_{ij}$的候选项是

结果的唯一性。假设对于所有的$(i, j)$对

如果$c_{max} - c_{min} < M$，我们可以在(13)中找到唯一的$C'$。为简单起见，令

如果

被满足，我们可以找到矩阵$X \in Z^{p \times r}$使得

如果$M$小于或等于$2c_{max}$，我们可能会找到多个结果的候选项（见图2）。尽管为简单起见，$c_{max}$被定义为所有元素的最大值，但也可以进行逐元素的定义。

缩放因子的确定。在(11)中的常数$k_A$和$k_B$对于计算结果的准确性非常重要，因为它与截断误差密切相关。这里我们解释一种找到$k_A$和$k_B$的简单方法。从(11)出发，

根据(16)，我们将(11)中的$k_A + k_B$设置为

然后，我们可以在(11)中设置对角矩阵$D$和$E$。如果$A$和$B$具有相同的精度，我们假设$k_A = k_B$，从而得到以下表达式：

缩放因子优化的说明。如果矩阵$A$和$B$的精度不同，例如，$A$是FP32，$B$是FP64，我们假设$k_A < k_B$。如果$s$增加，$M$也增加，结果是$k_A + k_B$变得更大。需要注意的是：
- 如果矩阵是稀疏的，$q$在(16)中可以减少。
- $q2^{k_A+k_B}$在(16)中是对$c_{max}$上限的高估。替代方法是使用柯西-施瓦茨不等式或使用低精度计算（见第3.1节）。
- 如果$A'B'$的范围已知，即$c_{min}$和$c_{max}$，我们可以改进(17)。

Ozaki方案II的整体流程。总的来说，Ozaki方案II包括以下四个部分：
- 第1部分：根据给定的$s$、$A$和$B$计算$k_A$和$k_B$。然后如(11)所示获得$A'$和$B'$。
- 第2部分：对$i = 1, \ldots, s$重复以下步骤
- 第2-a部分：如(12)所示确定$A'_i$和$B'_i$。
- 第2-b部分：使用BLAS例程计算$C_i := A'_iB'_i$。
- 第2-c部分：计算$Z := Z + C_iM_iy_i$。
- 第3部分：从矩阵$Z$中确定唯一的候选项$X$。
- 第4部分：应用(10)中的逆缩放$D^{-1}XE^{-1}$。
请注意，$M_iy_i$和$m_i$是预先计算并存储在一个表中的。在第2-b部分，我们根据矩阵$A$和$B$的结构使用适当的函数，如GEMM、TRMM或SYRK。在第3部分，我们通过范围缩减来找到(15)中的$x_{ij}$，因此在第2-c部分高精度计算至关重要。在Ozaki方案I中，我们需要保留$A_1, \ldots, A_k$和$B_k, \ldots, B_\ell$，这会消耗大量内存。然而，Ozaki方案II在获得$C_i$后立即丢弃矩阵$A'_i$和$B'_i$。

3.1 使用INT8 TCs进行矩阵乘法

模数的选择。如果我们在GPU上使用INT8 TCs，我们设置$2 \leq m_i \leq 256$（对于$1 \leq i \leq s$），并且$m_1, m_2, \ldots, m_s$两两互质。在这些条件下，我们准备尽可能大的$m_1, m_2, \ldots, m_s$。例如，当$s = 16$时，我们将$m_i$设置为：

这些数字和(9)中的$M_iy_i$被存储在一个表中，用于$s = 2, 3, \ldots$。与任何数模$m_i$的模运算结果都落在-128到127的范围内，这可以用INT8表示。在模256运算结果为128的情况下，INT8类型的回绕行为将此值映射为-128，从而避免了任何问题。然后，在使用cublasGemmEx时，对于$q < 2^{17}$，$A'_iB'_i$中不会出现错误，因为结果由INT32存储。如果$q \geq 2^{17}$，可以应用分块矩阵乘法来实现无误差的矩阵乘法。该代码已由Uchino (2025)开源【18, Uchino Y (2025) GEMMul8: GEMM emulation using int8 matrix engines based on the Ozaki Scheme II. R-CCS github repositry. URL https://github.com/RIKEN-RCCS/GEMMul8.】 。

精度与模数数量的关系。图3显示了对于$p = q = r \in \{1024, 4096, 16384\}$时，(18)中的$k := k_A = k_B$。对于$s = 16$，预计可以实现$k = 53$，从而允许用16次矩阵乘法进行FP64仿真。这里的期望是成立的，前提是矩阵元素之间的绝对值没有显著差异。为了达到与双精度算术相当的结果，Ozaki方案I需要7到8个切片，对应28到35次矩阵乘法。这一观察表明了Ozaki方案II的潜在优势。

算法1概述。以下是使用INT8张量核心的所提方法的概述。函数trunc(·)和round(·)分别以向零取整和四舍五入到最近偶数模式将输入舍入为整数，$E$是适当大小的全1矩阵。F64是FP64中的浮点数集合。

算法1：使用INT8张量核心的所提方法概述

# 将FP64转换为INT8：
1: 确定移位值 D ∈ F_p×p^64 和 E ∈ F_r×r^64
2: A' := trunc(DA)  {A' ∈ F_p×q^64 ∩ Z_p×q}
3: B' := trunc(BE)  {B' ∈ F_q×r^64 ∩ Z_q×r}
4: (a'_t)_ij := a'_ij mod m_t (1 ≤ t ≤ s)
5: (b'_t)_ij := b'_ij mod m_t (1 ≤ t ≤ s)

# 使用INT8 TCs进行矩阵乘法：
6: C'_t := A'_t B'_t (1 ≤ t ≤ s)

# 将矩阵乘积转换为UINT8：
7: (c''_t)_ij := (c'_t)_ij − ⌊(c'_t)_ij / m_t⌋ * m_t (1 ≤ t ≤ s)

# 累加矩阵乘积并进行逆缩放：
8: C''' := Σ_{t=1 to s} C''_t * M_y_t / m_t
9: C''' := C''' mod M
10: C := D⁻¹C'''E⁻¹

3.2 使用FP64进行矩阵乘法

模数的选择。我们将素数$m_1, m_2, \ldots, m_s$设置为满足

对于(12)中的$A'_t \in Z^{p \times q}$和$B'_t \in Z^{q \times r}$，使用(20)，我们有

对于所有$(i, j)$对和$1 \leq t \leq s$。因此，在使用BLAS中的GEMM时，$1 \leq i \leq s$的$A'_i B'_i$中不会出现舍入误差。类似地，考虑(16)，$k_A + k_B$的获取方式如(17)所示。

模数选择的进一步讨论。需要注意的是，$1 \leq i \leq s$的$m_i$不必是素数，只要它们两两互质即可。然而，与INT8 TCs的情况不同，当使用FP64时，为了确保能轻易找到足够大的素数，它们通常被选择为素数。例如，对于$s = 16$和$n = 2^{10}$，我们从(20)中设置$m_i$为：

原因是从(20)出发，

因此

注意(22)中的$m$满足

M的增长率分析。预期

因为$M$与$s$成正比，所以$k_A + k_B$也与$s$成正比。相反，(19)中的$m$满足

$M$的增长率随着$s$的增加而减小。

多字算术的应用。在此，我们处理FP64在多字算术中的应用。Ozaki方案II可以应用于任何多字格式；我们假设矩阵表示为一个由$v$个浮点矩阵组成的未求和的总和，例如

其中对于$1 \leq i \leq v - 1$，

同样，我们设置对角矩阵$D$和$E$使得

多字算术的具体计算。在第2-a部分，我们计算

其中对于所有的$i$，$A' \in Z^{p \times q}$和$B' \in Z^{q \times r}$。其余部分与原始方法相同。

多字算术下的精度。图4展示了常数$k := k_A = k_B$作为模数$s$的函数的行为。当$s$在16到20的范围内时，$k$的值大约为160，这对应于与三字算术获得的精度水平相当。当$s$增加到21到25的范围时，$k$增加到大约210，表明精度水平可与四字算术相媲美。

扩展到多矩阵乘积。到目前为止，我们已经考虑了两个矩阵的乘积；然而，同样的方法可以扩展到三个或更多矩阵的乘积。例如，在处理三个矩阵$A、B$和$C$的乘积时，我们将它们转换为

其中矩阵$D、E$和$F$是对角矩阵，其对角线元素是2的幂，并且$k_A、k_B$和$k_C$被确定为满足

然后，可以应用与(16)和(17)中类似的过程。如果$k_A = k_B = k_C$，我们设置

A4 实验环境

本研究在GPU和CPU两个平台上进行了数值实验，以评估所提方法的准确性和计算性能。所有实验均假设输入矩阵具有相同的精度水平，即$k_A=k_B$。

• GPU实验环境 (使用INT8 TCs仿真FP64)

  • 硬件配置:
    • NVIDIA GH200 Grace Hopper Superchip
    • NVIDIA GeForce RTX 4090 GPU
  • 软件配置:
    • NVIDIA CUDA Toolkit 12.8.61
  • 数据集:
    • 测试矩阵$A \in F^{p \times q}$ 和 $B \in F^{q \times r}$通过以下公式生成：
      
      其中$\phi$控制指数分布，rand是(0, 1]内的均匀分布随机数，randn是标准正态分布随机数。$\phi=0.5$时经验上与HPL基准测试中的指数分布相当。
  • 对比方法:
    • DGEMM: 使用cublasGemmEx的原生FP64矩阵乘法。
    • OS II-fast-s: 本文提出的算法1，使用$s$个模数，并通过柯西-施瓦茨不等式估计上界以保证结果唯一。
    • OS II-accu-s: 本文提出的算法1，使用$s$个模数，并通过低精度（INT8）cublasGemmEx直接计算上界以保证结果唯一。
    • ozIMMU EF-S: 传统的Ozaki方案I（由Uchino (2024)实现【17, Uchino Y (2024) Accelerator for ozIMMU. R-CCS github repositry. URL https://github.com/RIKEN-RCCS/accelerator_for_ozIMMU.】），使用$S$个切片 。

• CPU实验环境 (使用FP64仿真四倍字长算术)

  • 硬件配置:
    • Intel® Core™ i7-8665U 处理器 (4核)
    • Intel® Core™ i9-10980XE 处理器 (18核)
  • 软件配置:
    • MATLAB 2024b, Windows 10
    • 代码实现为MATLAB可执行文件(MEX)，使用Visual Studio 2022编译。
    • 编译标志: Core i7使用/openmp和/arch:AVX2；Core i9使用/openmp和/arch:AVX512。
  • 数据集:
    • 输入矩阵以四倍字长（$v=4$）格式表示：
      
      其中$A^{(1)}$和$B^{(1)}$由MATLAB的randn函数生成，并满足多字格式的非重叠条件。
  • 内部计算精度:
    • Ozaki方案II的第2-c部分和第3部分使用六倍字长算术来保证累加和归约的精度。

A4 实验结果

4.1 使用INT8 TCs进行FP64仿真（GPU）

• 精度分析 (图5):

  • 实验内容: 在GH200上比较DGEMM、OS II-fast-s和OS II-accu-s的精度，其中$s$从8到20变化。
  • 实验结果: OS II-accu-s由于对乘积上界的估计更精确（开销更大），其结果比OS II-fast-s更准确。对于$\phi=0.5$（模拟HPL场景），两种提出的方法都需要14或15个模数才能达到原生DGEMM的精度水平。OS II-accu-s在处理更大$\phi$值（更具挑战性的数据分布）时表现更佳。
  • 结论: 所提出的Ozaki方案II能够通过调整模数数量达到与原生FP64相当的精度。
    

• 性能分析 (表3, 表4):

  • 实验内容: 比较DGEMM, ozIMMU EF-8, OS II-fast-s和OS II-accu-s在不同矩阵尺寸下的吞吐量（TFLOPS）。
  • 实验结果:
    • 对于小尺寸问题，由于INT8 TCs未能充分发挥性能且额外开销较大，仿真方法的性能低于DGEMM。
    • 对于大尺寸问题（例如16384x16384），在GH200上，OS II-fast-14达到了80.2 TFLOPS，超过了DGEMM的60.9 TFLOPS（表3）。在RTX 4090上，OS II-fast-14达到了9.81 TFLOPS，远超DGEMM的0.62 TFLOPS（表4）。
  • 结论: 在大矩阵上，所提出的方法在达到DGEMM级别精度的同时，实现了比原生FP64更高的吞吐量。

• 时间开销分解 (图6, 图7):
  • 实验内容: 分析OS II-fast-s和OS II-accu-s在不同矩阵尺寸和模数下的时间构成。
  • 实验结果:
    • 小尺寸问题: 性能瓶颈是内核启动开销。
    • 中等尺寸问题: 瓶颈主要是从双精度到INT8的矩阵转换，以及矩阵乘积的累加和逆缩放开销。
    • 大尺寸问题: 瓶颈仍然是双精度到INT8的转换，而内核启动和累加的开销几乎可以忽略不计。
  • 结论: 性能优化的关键在于减少数据类型转换的开销。
    
    

4.2 使用FP64进行四倍字长算术仿真（CPU）

• 精度分析 (图8):

  • 实验内容: 在Core i9上比较Ozaki方案I和方案II在不同矩阵尺寸下的最大相对误差。
  • 实验结果: 当矩阵乘法次数较少时，方案I的精度更高。但随着乘法次数增加，方案II的精度改善更快，最终超越方案I。这是因为方案II的精度在对数尺度上与矩阵乘法次数呈线性关系。
  • 结论: 在需要更高精度时，方案II比方案I更有效率。
    

• 时间开销分解 (图9, 图10):

  • 实验内容: 分析方案II在Core i7和Core i9上不同矩阵尺寸下的各部分计算时间占比。
  • 实验结果: 对于小矩阵（n=1000），除矩阵乘法外的其他开销占主导。随着矩阵尺寸增大，矩阵乘法（Part 2-b）的时间占比显著提高，表明方法越来越依赖于GEMM的性能。
  • 结论: 该方法的设计目标——以GEMM为计算核心——在大规模问题上得以实现。
    
    

• 性能分析 (表5, 表6):

  • 实验内容: 比较方案I和方案II在Core i7和Core i9上的吞吐量（GFLOPS）。以精度相当的OS I-10和OS II-22为例。
  • 实验结果: 在小矩阵（n=1000）上方案II较慢，但在较大矩阵上性能更优。在Core i7上，当n=4000时，方案II实现了高达2.29倍的加速；在Core i9上，当n=8000时，实现了高达2.12倍的加速。
  • 结论: 对于中到大规模问题，Ozaki方案II在CPU上仿真高精度算术方面比传统方案I具有显著的性能优势。

A5 结论

本研究成功开发了一种名为Ozaki方案II的新型矩阵乘法仿真算法。该算法结合了中国剩余定理和基于GEMM的无误差计算思想，能够高效地利用现代硬件（特别是低精度矩阵引擎）来执行高精度计算。实验结果表明，该方法在多组分算术（如四倍精度）仿真方面也表现出色。

未来的工作将集中在以下几个方面：
1. 实现优化：进一步优化当前实现的性能，特别是解决数据类型转换等瓶颈问题。
2. 误差分析：对所提方法引入的舍入误差进行详细的理论分析。
3. 探索下一代架构：计划在具有更高INT8张量核心性能的GPU（如NVIDIA B200）上进行进一步实验，以挖掘该方法在未来硬件上的潜力。