Ozaki Scheme

背景

近几年 NVIDIA 的计算卡规格有一个现象：低精度和高精度算力虽然都在增长，但总体趋势上，低精度相对高精度的加速比在扩大，不过这个比例并非跨代单调递增。以相同的 INT8/FP32 口径比较，A100 上 INT8 Tensor Core（624 TOPS）相对 FP32 标量核（19.5 TFLOPS）是 32 倍，H100 约为 30 倍，到了 B200 则约为 56 倍。这也表明，用低精度模拟高精度计算会越来越有性价比，因为低精度算力总体上比高精度增长得更快。

值得注意的是，从 Blackwell 开始，官方规格中开始并列给出基于 Tensor Core emulation 的 FP64 吞吐：B200 的模拟峰值为 150 TFLOPS，而原生 FP64 为 37–40 TFLOPS；到 Rubin，官方 Table 3 已明确标注 emulation 峰值。这类方法使用 Tensor Core 对低精度计算做多次分解，从软件层面模拟高精度计算，也就是 Ozaki 类算法的技术路线。

精度 / 核心类型	A100 SXM 80GB (Ampere)	H100 SXM5 (Hopper)	H800 SXM5	B200 SXM (Blackwell)	Rubin R200
CUDA Core FP64	9.7	33.5	≈1 ᵇ	37–40	33
CUDA Core FP32	19.5	66.9	66.9(未砍)	80	130
Tensor Core FP64(原生硬件)	19.5	66.9	≈1 ᵇ	37–40(与CUDA Core同级,无独立加速电路)	未公布
Tensor Core FP64(Ozaki类算法模拟) ᵃ	-	66.9(原生硬件已达此值,无需模拟)	同H100	150(软件模拟)	200(软件模拟)
Tensor Core TF32	156 / 312	494.7	同H100	1,100	未公布
Tensor Core FP32(模拟真FP32精度) ᵃ	-	66.9(同CUDA Core,无独立矩阵加速)	同H100	227(软件模拟)	400(软件模拟)
Tensor Core FP16 / BF16	312 / 624	989.4	同H100	2,250	未公布
Tensor Core FP8	-	1978.9	同H100	4,500	≈16,000(非官方,dense)
Tensor Core INT8	624	1978.9	同H100	4,500	未公布
Tensor Core FP4 / INT4	1,248	-	-	9,000	训练35,000 / 推理50,000

单位:TFLOPS(浮点)/ TOPS(整数),均为单卡(单GPU封装)峰值理论值,不含互联带宽。

表格数据来源与口径说明

• ᵃ Ozaki 类算法模拟:数据来自NVIDIA官方2026年技术博客《Inside the NVIDIA Vera Rubin Platform》Table 3,原文明确标注”Peak performance using Tensor Core-based emulation algorithms”。即不是专用高精度电路直接输出,而是用低精度Tensor Core做多趟分解运算、软件层面模拟出高精度累加结果，和Ozaki scheme是同一类技术路线,现在已经被NVIDIA写进官方性能指标。Hopper一行没有这个标注,因为H100的FP64 Tensor Core是真正的专用硬件电路,原生吞吐本身就是66.9,不需要模拟。
• ᵇ H800 FP64:NVIDIA从未正式公布精确数字,业界(Lenovo/电子分销商等)估计在1 TFLOPS量级,不同渠道口径有出入,仅供参考。H800与H100的低精度Tensor Core(TF32/FP16/BF16/FP8/INT8)吞吐完全一致,只有FP64和NVLink带宽(900→400 GB/s)被砍。

---

基本算法（Ozaki Scheme）

来源

此处的算法实现主要参考 Ootomo、Ozaki 和 Yokota 的论文（arXiv 2023，IJHPCA 2024 刊出），“DGEMM on Integer Matrix Multiplication Unit”，和 Ozaki scheme 算法的代码实现 ozIMMU.

记号约定

基础算法沿原文使用 $s$ 表示切片数、$\alpha$ 表示每片有效位数；优化算法沿 Uchino 等人的论文使用 $k$ 和 $\beta$。$u_{\mathrm{acc}}$ 表示累加器的单位舍入误差。

核心思想

低精度模拟高精度的基本思想：把高精度数值（如 FP64）的尾数按固定长度切成若干段，每段对应不同的位权（即不同的 2 的幂次缩放）。把每段数值乘上对应位权再求和，即可无损还原原高精度数值。由于每段长度固定（如 8 bit），拆分出的每段本身就是一个 INT8 整数，两数相乘时也就退化为 INT8×INT8 的整数乘法。

矩阵拆分算法 SplitFP

Input: An $m\times k$ matrix $\mathbf{M}$ and the number of slices $s$

Output: Matrix slices ${\mathbf{M}}^{(1)},{\mathbf{M}}^{(2)},\cdots ,{\mathbf{M}}^{(s)}$

$$\begin{array}{rl}1.& \alpha =\lfloor \frac{-{\log }_2{u}_{\mathrm{acc}}-{\log }_{2}k}{2}\rfloor \\ 2.& \mathbf{R}:=\mathbf{M}\\ 3.& \text{for }p=1,\dots ,s-1\text{ do}\\ 4.& \text{Compute }{\mathbf{m}}^{(p)}\text{ where }{m}_i^{(p)}\leftarrow \underset{j}{\max }|R_{i,j}|\\ 5.& \text{Compute }{\sigma }^{(p)}\text{ where }{\sigma }_i^{(p)}\leftarrow \{\begin{array}{ll}0.75\times 2^{\lceil {\log }_2{m}_i^{(p)}\rceil }\times 2^{54-\alpha },& \text{if }{m}_i^{(p)}\ne 0\\ 0,& \text{else}\end{array}\\ 6.& \text{Compute }{\mathbf{M}}^{(p)}\text{ where }{M}_{i,j}^{(p)}\leftarrow \left(R_{i,j}+{\sigma }_i^{(p)}\right)-{\sigma }_i^{(p)}\\ 7.& \mathbf{R}=\mathbf{R}-{\mathbf{M}}^{(p)}\\ 8.& \text{end for}\\ 9.& {\mathbf{M}}^{(s)}=\mathbf{R}\end{array}$$

算法解释：

1. 计算得到每个切片的 bit 数。这里的 $u_{\mathrm{acc}}$ 是累加器的单位精度，对于 FP32 来说就是 $u_{\mathrm{acc}}=2^{-24}$，所以 $-{\log }_2{u}_{\mathrm{acc}}=24$，这意味着累加器大约有 24 bit 的“预算”空间。现在来看这 24 bit 预算会被谁占用：

   1. 乘法本身要花 $\alpha$ bit。每个切片是 $\alpha$ bit 宽的 INT 整数，量级大约是 $2^{\alpha }$。两个切片相乘 $A^{(i)}\cdot B^{(j)}$，量级就变成 $2^{\alpha }\times 2^{\alpha }=2^{2\alpha }$，也就是需要 $2\alpha$ bit 才能精确装下这个乘积。

   2. 对 k 个乘积求和还需要再花 ${\log }_{2}k$ bit。一次 GEMM 内部会把内积中 reduction 维度上的 k 个乘积结果累加成一个数。k 个量级为 $2^{2\alpha }$ 的数加在一起，最坏情况下总和能达到 $k\times 2^{2\alpha }=2^{2\alpha +{\log }_{2}k}$，所以还需要再预留 ${\log }_{2}k$ bit 的进位空间，k 越大，需要的余量就越大。 所以要想整个累加过程绝对不溢出（即 error-free），就需要“总开销 $\le$ 总预算”：

      $$2\alpha +{\log }_{2}k\le -{\log }_2{u}_{\mathrm{acc}}$$

      解出 $\alpha$：

      $$\alpha \le \frac{-{\log }_2{u}_{\mathrm{acc}}-{\log }_{2}k}{2}$$

      因为 $\alpha$ 必须是整数，而这个不等式给出的是上界，所以只能向下取整，取满足不等式的最大整数：

      $$\alpha =\lfloor \frac{-{\log }_2{u}_{\mathrm{acc}}-{\log }_{2}k}{2}\rfloor$$

2. 将残差矩阵 $R$ 初始化为输入矩阵 $\mathbf{M}$，$R$ 相当于是一个工作副本，之后每轮循环都会从里面取走一层 $\alpha$ bit，剩下的继续留在 $R$ 里面留给下一轮使用。

3. 操作 $s-1$ 次。

4. 取残差矩阵的每一行的最大绝对值。因为之后是矩阵的行和另一个矩阵的列做点积，所以矩阵的每一行共用一套指数参考，而不同行之间可以有各自的参考，这是一个比整个矩阵共用一个指数更精细的方案，能减少大范围输入造成的浪费。

5. 构建大指数 $\sigma$。

$${\sigma }_i^{(p)}\leftarrow 0.75\times 2^{\lceil {\log }_2{m}_i^{(p)}\rceil }\times 2^{54-\alpha }$$

关于 0.75 与 ±1 位余量的细节

$\sigma$ 的位置由 ceiling 后的行最大值指数决定，这个向上取整本身带来约 1 位的保守余量，因此每轮实际挤出来的是约 $\alpha$ 位（±1，取决于舍入是否发生进位）。实际上这个算法在原论文里面也被描述为 “each element substantially keeps slightly more bits than α when using RN for the rounding”，只能说大致 $\alpha$ 位，用 RN 舍入时略多一点。所以不能达到恰好 $\alpha$ 位。但是由于 $\alpha$ 的取值本身自带了一层保守余量，已经能保证 INT32 累加器不会溢出，所以这里 ±1 的误差已经被覆盖了。

0.75 的作用是让 $\sigma$ 落在自己所在二进制区间的正中间，从而保证 $R_{i,j}+\sigma$（无论 $R_{i,j}$ 正负）都不会意外跨越区间边界导致舍入网格突变；这跟 ceiling 带来的 ±1 位模糊是两个独立的问题，0.75 并不解决后者。

6. $M_{i,j}^{(p)}\leftarrow \left(R_{i,j}+{\sigma }_i^{(p)}\right)-{\sigma }_i^{(p)}$ 用了一个加上再减掉的方法，利用的是硬件的舍入机制。因为 $\sigma$ 的量级比 $R_{i,j}$ 大得多（大出约 $54-\alpha$ 个 bit 的指数差），$R_{i,j}+{\sigma }_i^{(p)}$ 这一步浮点加法会被迫按 $\sigma$ 的精度网格舍入，也就是说 $R_{i,j}$ 低于这个网格的 bit 全部被舍掉了；再减掉 ${\sigma }_i^{(p)}$，剩下的正好是 $R_{i,j}$ 这一行相对当前最高位往下数、宽度约为 $\alpha$ bit 的那一段，也就纯靠浮点加法自身的舍入规则完成“切片”。
7. 将刚刚抠出来的那一层从残差里面减掉，剩下更低位的部分留给下一轮继续切。
8. 循环结束后，最后一片直接等于剩下的残差。

FP 矩阵乘法算法

原来的 $A,B$ 被拆分成 $s$ 个 $A$ 和 $t$ 个 $B$，原来的矩阵乘法被拆分成了 $s$ 个矩阵和 $t$ 个矩阵的乘加。

论文中的这张图展示了两种切分向量的方法：左图是逐元素切分，$a_1,a_2,...,a_k$ 是待拆分向量里的元素，各自保留自己的浮点指数 $2^{e_i}$。对于不同元素而言，这些切片对应的实际数值位权是不同的，因为每个元素的 $e_i$ 不一样。这种切法的问题就是：当我想用整数矩阵乘法去算 ${\sum }_k{a}_{\cdot k}^{(i)}{b}_{k\cdot }^{(j)}$ 时，两个向量中不同 $k$ 的切片对应的位权不同，没法用一个统一的 scale factor 去恢复真实值，也就没法用若干次干净的、无损的 Tensor Core GEMM 来计算，也不能靠增加切片数量逼近目标精度。

右图就是 Ozaki scheme 采用的方式：这里的切分点不再基于每个元素自己的尾数起始位，而是基于一个所有元素共享的、按真实数值位权对齐的坐标轴。也就是说，只要切片编号相同，其位权在整个向量里就是统一的，于是 $A^{(i)}\cdot B^{(j)}$ 可以直接用一次整数矩阵乘法完成，再按照 $(i,j)$ 对应的统一 scale 做加权求和，从而实现无误差变换。

Input: Matrices $\mathbf{A},\mathbf{B}$ and the number of slices $s$

Output: $\mathbf{C}\leftarrow \mathbf{A}\cdot \mathbf{B}$

$$\begin{array}{rl}1.& {\mathbf{A}}^{(1)},{\mathbf{A}}^{(2)},\cdots ,{\mathbf{A}}^{(s)}\leftarrow \mathrm{SplitFP}(\mathbf{A},s)\\ 2.& {\mathbf{B}}^{(1)},{\mathbf{B}}^{(2)},\cdots ,{\mathbf{B}}^{(s)}\leftarrow \mathrm{SplitFP}(\mathbf{B},s)\\ 3.& \mathbf{C}=0\\ 4.& \text{for }i=1,\dots ,s\text{ do}\\ 5.& \text{for }j=1,\dots ,s-i+1\text{ do}\\ 6.& {\mathbf{C}}_{\mathrm{tmp}}\leftarrow {\mathbf{A}}^{(i)}\cdot {\mathbf{B}}^{(j)}//\text{Low-precision. No rounding error}\\ 7.& \mathbf{C}\leftarrow \mathbf{C}+{\mathbf{C}}_{\mathrm{tmp}}//\text{High-precision accumulation}\\ 8.& \text{end for}\\ 9.& \text{end for}\end{array}$$

算法解释：

1. 先将 $A,B$ 进行切片，分别得到 s 个矩阵切片（为了方便起见，这里将 A 和 B 的切片数量设为相同）
2. 第 $i$ 片 ${\mathbf{A}}^{(i)}$ 的量级约为最高片的 $2^{-\alpha (i-1)}$ 倍。所以乘积项 ${\mathbf{A}}^{(i)}{\mathbf{B}}^{(j)}$ 相对于最大项 ${\mathbf{A}}^{(1)}{\mathbf{B}}^{(1)}$ 的量级约为 $2^{-\alpha (i+j-2)}$。
3. 通过 $i=1,\dots ,s$ 和 $j=1,\dots ,s-i+1$ 控制 $i+j\le s+1$。因为当 $i+j>s+1$ 的时候，这些项已经低于目标精度能表示的最低位了，所以论文说它们 “only have a negligible effect on the resulting accuracy”。这里其实存在截断误差，但是可以通过增加切片数 $s$ 将其压到任意小，直至小于目标精度（如 FP64）的舍入误差水平。论文中 FP64 通常需要 $s=8\sim 13$。
4. 针对每个 $i,j$ 做 GEMM。注意，这里每个切片由于尾数较短，可以视为较低精度的浮点数（如FP16）进行计算，浮点数自带指数，位权在浮点乘法里被硬件自动处理了，所以不需要显式的指数缩放。
5. 由于这里指数缩放都由硬件管理，所以 GEMM 和累加过程中都不需要显式的指数缩放。

相关参数：

• $s$ 是拆分后的矩阵数量，一般来说是静态地按照目标精度计算得到的：切片数量是精度、计算量、显存三者的 trade-off，具体来说，整个切片组能保留的精度等于：$\text{mantissa space}=s\times \mathrm{BPS}$（BPS 就是每片的有效 bit 数）。但这里的目标 mantissa space 不是简单地等于原始格式的尾数长度（比如 FP64 的 53 bit），由于 shared-place splitting 是按照行内最大值统一定坐标轴的，所以实际需要的 mantissa space = 目标精度位数 + 行内指数分布宽度带来的额外开销。

矩阵拆分算法 SplitINT

考虑到既然向量采用 shared-place splitting 进行切分，那向量中每个元素其实共用一个指数位，自然也就想到能将尾数切片拆出来作为 INT 参与矩阵乘法，从而转化为 IMMU (Integer matrix multiplication unit) 上的计算。

Input: An $m\times k$ matrix $\mathbf{M}$ and the number of slices $s$

Output: Matrix slices ${\mathbf{M}}^{(1)},{\mathbf{M}}^{(2)},\cdots ,{\mathbf{M}}^{(s)}$ and an $m$-dimensional exponent-offset vector $\mathbf{e}$

$$\begin{array}{rl}1.& \alpha =\lfloor \frac{-{\log }_2{u}_{\mathrm{acc}}-{\log }_{2}k}{2}\rfloor \\ 2.& \text{Compute }\mathbf{e}\text{ where }{e}_i\leftarrow \underset{j}{\max }{2}^{\lfloor {\log }_2|R_{i,j}|\rfloor }\\ 3.& \text{for }p=1,\dots ,s\text{ do}\\ 4.& \text{Compute }{\mathbf{M}}^{(p)}\text{ where}\\ & \text{the sign of }{M}_{i,j}^{(p)}\text{ is the same as }{M}_{i,j}\\ & \text{the mantissa of }{M}_{i,j}^{(p)}\text{ is part of }{M}_{i,j}\text{ mantissa,}\\ & \text{which corresponds to the part from }{2}^{(p-1)\alpha }\text{ bit to }{2}^{p\alpha }-1\text{ bit}\\ & \text{in the mantissa space with }{e}_i\text{ as the first bit, as shown at the right of Fig. 2.}\\ 5.& \text{end for}\end{array}$$

算法解释：

1. 基本思想和 SplitFP 一致，先算出每个切片的 bit 数，然后算出每一行的最大指数，作为共享指数。
2. 进行切片循环，相当于是将一个长尾数按 $\alpha$ 位一段切成 $s$ 段，高位段在前、低位段在后。

代码实现：

• 实际 ozIMMU 实现这个算法的时候，是这样操作的：

  1. 先通过判断浮点数正负确定符号

  2. 取出浮点数的尾数部分，并放到一个 128-bit 整数的最顶端，方便从最高位开始切

  3. 将当前浮点数的指数位和最大指数的指数位相减，再右移尾数长度的位数，得到纯指数差 $offset=E_{max}-E_a$。再将当前浮点数的尾数右移 $offset$ 位，即按照 max_exp 对齐。元素越小，右移越多，有效位就会落到共同标度的更低段里面。

  4. 将之前得到的 128-bit 数右移 $(128-b)$ 位，即取最高 $b$ 位，因为 significant 是无符号，切出来的是 $[0,2^b-1]$ 的非负整数。再将 128-bit 数左移 $b$ 位，窗口下移，下一轮取下一段。

  5. 最终结果可以表示为

     $$a\approx \mathrm{sign}\cdot 2^{E_{\max }}\cdot \left({int8}_0\cdot 2^{-b}+{int8}_1\cdot 2^{-2b}+\cdots \right)$$
  ozIMMU cut_int8_core 函数源码

  template <class INPUT_T, class MANTISSA_T>
  __device__ void cut_int8_core(std::int8_t *const out_ptr, const std::size_t inc,
                                const INPUT_T a, const INPUT_T max_exp,
                                const unsigned num_split,
                                const unsigned mantissa_length) {
    const std::uint8_t sign_flag = a > 0;
    // When the input number is not normalized, don't set the implicit one bit.
    const std::uint64_t implict_one_bit =
        cutf::experimental::fp::mask_exponent(a) ? 1lu : 0lu;
    const auto mantissa =
        static_cast<MANTISSA_T>(
            cutf::experimental::fp::mask_mantissa(a) |
            (implict_one_bit
             << cutf::experimental::fp::get_mantissa_size<INPUT_T>()))
        << ((sizeof(MANTISSA_T) - sizeof(INPUT_T)) * 8 +
            cutf::experimental::fp::get_exponent_size<INPUT_T>());
    const auto mantissa_shift_offset =
        (cutf::experimental::fp::reinterpret_as_uint(max_exp) -
         cutf::experimental::fp::mask_exponent(a)) >>
        cutf::experimental::fp::get_mantissa_size<INPUT_T>();
   
    auto shifted_mantissa = mantissa >> mantissa_shift_offset;
    for (unsigned s = 0; s < num_split; s++) {
      const std::int8_t int8 =
          static_cast<std::int8_t>(shifted_mantissa >>
                                   (sizeof(MANTISSA_T) * 8 - mantissa_length)) *
          (sign_flag ? 1 : -1);
      shifted_mantissa <<= mantissa_length;
   
      out_ptr[s * inc] = int8;
    }
  }
   
  __device__ cuDoubleComplex x2(const cuDoubleComplex a) {
    return make_cuDoubleComplex(a.x * 2, a.y * 2);
  }

• 这张图可以直观地解释 SplitFP 算法和 SplitINT 算法的区别

  注意：在图中的右边部分，段编号是从低位往高位段数的（$a^{(1)}$ 最低，$a^{(s)}$ 最高），量级逐段升高，所以缩放指数逐段 $+\alpha$。但是实际代码中是从高位段往低位段切的。

INT 矩阵乘法算法

Input: Matrices $\mathbf{A},\mathbf{B}$ and the number of slices $s$

Output: $\mathbf{C}\leftarrow \mathbf{A}\cdot \mathbf{B}$

$$\begin{array}{rl}1.& \left({\mathbf{A}}^{(1)},{\mathbf{A}}^{(2)},\cdots ,{\mathbf{A}}^{(s)}\right),{\mathbf{e}}_A\leftarrow \mathrm{SplitInt}(\mathbf{A},s)\\ 2.& \left({\mathbf{B}}^{(1)},{\mathbf{B}}^{(2)},\cdots ,{\mathbf{B}}^{(s)}\right),{\mathbf{e}}_B\leftarrow \mathrm{SplitInt}(\mathbf{B},s)\\ 3.& \mathbf{C}=0\\ 4.& \text{for }i=1,\dots ,s\text{ do}\\ 5.& \text{for }j=1,\dots ,s-i+1\text{ do}\\ 6.& {\mathbf{C}}_{\mathrm{tmp}}\leftarrow {\mathbf{A}}^{(i)}\cdot {\mathbf{B}}^{(j)}//\text{Integer}\\ 7.& \mathbf{C}\leftarrow \mathbf{C}+{\mathbf{C}}_{\mathrm{tmp}}\odot \left(2^{-(i+j)\alpha }\cdot {\mathbf{e}}_A\cdot {\mathbf{e}}_B^{\mathrm{T}}\right)\\ 8.& \text{end for}\\ 9.& \text{end for}\end{array}$$

算法解释：

1. 与 FP 矩阵乘法类似，先将 $A,B$ 分别拆分。

2. 同样采用三角形筛选，控制 $i+j\le s+1$。

3. 将 $A^{(i)},B^{(j)}$ 扔给 INT8 Tensor Core，累加到 INT32，这个过程无舍入误差。

4. 单元素展开：

   $$a\approx e_{A,\mathrm{row}}\cdot \sum _{p=1}^s{I}_p\cdot 2^{-p\alpha },b\approx e_{B,\mathrm{col}}\cdot \sum _{q=1}^s{J}_q\cdot 2^{-q\alpha }$$

   两式相乘：

   $$a\cdot b\approx e_{A,\mathrm{row}}\cdot e_{B,\mathrm{col}}\cdot \sum _{p,q}{I}_p{J}_q\cdot 2^{-(p+q)\alpha }$$

   对应到矩阵形式，$i,j$ 就是 $p,q$，${\mathbf{e}}_A\cdot {\mathbf{e}}_B^{\mathrm{T}}$ 做外积将行指数和列指数分别广播到每个元素，$\odot$ 是逐元素乘。

为什么选择 INT8-INT32

结论： 在论文关注的 $[2^{11},2^{20}]$ 矩阵规模范围内，INT8-INT32 是最适合 Ozaki scheme 的 IMMU 组合：INT32 累加器允许单片使用更多有效位，从而减少切片数、显存占用和 GEMM 次数。

推导链条是：

$$INT32\text{ 累加器更宽}\Rightarrow \mathrm{BPS}\text{ 更大}\Rightarrow s\text{ 更少}\Rightarrow \{\begin{array}{l}\text{内存占用}\propto s\text{，更低}\\ \mathrm{GEMM}\text{ 次数}\propto s^2\text{，更少}\end{array}$$

---

基础实现的 Profiling

结论： INT8x 是否更快取决于 GPU 的原生 FP64 能力；它在 FP64 较弱的消费级或受限 GPU 上更有优势，在 A100 上则会被多次整数 GEMM 的开销抵消。

论文中，作者分别在计算卡和消费卡上测试了吞吐量和功耗。结果表明，INT8 的优势取决于 GPU 的 FP64 算力有多强：

• 在 A100 上（FP64 强）：DGEMM 能跑到 FP64 Tensor Core 理论峰值的 90% 以上，而 INT8x 反而慢 3–5 倍。原因是 INT8x 内部需要做 45–91（$s=9\sim 13$）次整数 GEMM，而 INT8 Tensor Core 的峰值只有 FP64 Tensor Core 的约 32 倍，所以净结果是变慢。
• 在其他消费级 GPU（FP64 被阉割）上：这些卡的 FP64 理论峰值很低，DGEMM 本身就跑得很慢，此时 INT8x 反而更快。

这也展示了低精度模拟高精度的核心 motivation：消费级以及特供版 GPU 的 FP64 单元被厂商刻意砍得很弱，但 INT8 Tensor Core 算力依然很强，于是低精度模拟高精度在这些卡上有很高性价比。

结论： 基础实现的主要瓶颈是 INT8 GEMM 数量和 memory-bound 的 FP64 累加；后续优化应优先减少 FP64 累加次数。

同时，作者还用 ncu 做了时间 profiling，发现 INT8 GEMM 和 FP64 累加是两大时间消耗。用 ncu 深挖 FP64 累加 kernel 发现：$m\le 2^{11}$ 时 kernel 没有完全利用 GPU 资源，表明性能欠佳；$m\ge 2^{12}$ 时 kernel 用到约 90% 显存带宽，说明 FP64 累加是 memory-bound（累加过程中每个元素需要读写 3 次，但是只做一次缩放乘 + 一次加）。所以想要进一步提速，需要修改算法、减少 FP64 累加次数。

---

优化算法（Uchino et al. 2024）

来源

此处的优化主要参考 Uchino、Ozaki 和 Imamura 的论文（arXiv 2024，IJHPCA 2025 刊出），Performance Enhancement of the Ozaki Scheme on Integer Matrix Multiplication Unit。

基本算法存在的优化空间

1. 切片可能存在较多前导零，浪费有效位
2. 仍然存在截断误差，而这个截断误差可以再减小

方向 1：bitmask → round-to-nearest

• 基本算法中是每次截取同样长度的一段尾数，从最高位开始每 $\alpha$ 位切一刀。而 round to nearest 是每次对残差做 round-to-nearest，只保留 $\alpha$ 位有效位的结果，然后再对残差做同样的事情。也就是每一片都取残差的最近可表示值（而不是截断），而且允许出现负的分片（这在 INT8 是有符号的前提下没有问题）。
• 通俗一点说就是：每次尝试用有限的位数尽可能接近残差值。

Magic constant

具体算法中会使用一个 “magic constant” 实现舍入。

• 先介绍 ulp（unit in the last place）这个概念：IEEE-754 里面，浮点数不是均匀分布的，在区间 $[2^e,2^{e+1})$ 内，所有可表示的浮点数间距是恒定的，这个间距就叫 ulp：

  $$\mathrm{ulp}([2^e,2^{e+1}))=2^{e-(p-1)}$$

  其中 $p$ 是有效位数（FP64 是 $p$ = 53）。指数 $e$ 每增加 1，网格就加粗一倍。所以网格间距是被结果落在哪个指数区间决定的。 现在看一个特殊区间。取 $e=p-1=52$，即区间 $[2^{52},2^{53})$。代入得：

  $$\mathrm{ulp}=2^{52-52}=2^0=1$$

  这个区间里面，浮点数间距恰好是 1，也就是说这段里能表示的全是整数，且相邻整数都能表示。只要能把一个数搬进这个区间，浮点加法器就会自动把它舍入到最近的整数。

• 这时候，我们取 magic constant $\sigma =c\cdot 2^{52}$。考虑 $a+\sigma$。如果 $|a|$ 不太大，$a+\sigma$ 落在 $[2^{52},2^{53})$ 这个 $\mathrm{ulp}=1$ 的区间里面，浮点加法器会完成最近整数舍入：

  $$\mathrm{fl}(a+\sigma )=\sigma +\mathrm{round}(a)$$

  减回 $\sigma$ 后即可得到舍入结果：

  $$\mathrm{fl}(a+\sigma )-\sigma =\mathrm{round}(a)$$

为什么取 c = 1.5：两种 σ 取值的落点分析

• 若 $\sigma =2^{52}$（即 $c=1$）：当 $a>0$，$a+\sigma \in [2^{52},2^{53})$，$\mathrm{ulp}=1$。但当 $a<0$，$a+\sigma <2^{52}$，会进入 $[2^{51},2^{52})$，此时 $\mathrm{ulp}=0.5$，负数会被舍入到最近的 $0.5$ 倍数，而不是整数。
• 若 $\sigma =1.5\cdot 2^{52}$（即 $c=1.5$）：$\sigma$ 落在区间正中间。只要 $|a|<2^{51}$，$a+\sigma$ 就始终位于 $[2^{52},2^{53})$，正负数都会落在 $\mathrm{ulp}=1$ 的整数网格上。

因此，具体算法取 $c=1.5$，给 $a$ 留出上下各 $2^{51}$ 的安全余量。对应的使用条件是：

$$|a|<2^{51}\Rightarrow \text{magic constant 给出正确的 round-to-nearest-even}$$

• 推广到任意网格间距 $\delta =2^t$，直接令：

  $$\sigma =1.5\cdot 2^{52}\cdot \delta =1.5\cdot 2^{52+t}$$

  此时 $\sigma$ 所在区间的 ulp 恰好等于 $\delta$，于是 $(a+\sigma )-\sigma$ 给出 $a$ 在 $\delta$ 网格上的最近点。这正是论文中的写法：$\sigma =0.75\cdot 2^{53}\cdot \mu =1.5\cdot 2^{52}\cdot \mu$。使用条件相应缩放为 $|a|<2^{51}\cdot \delta$。

算法分析

• 我们先确定切片的长度 $\beta$，有两个限制：

  1. 我们的目标是将浮点数切分为 INT8，INT8 是有符号的，符号占一位，单片最多装 7 位有效数字。

  2. $\lfloor \frac{31-{\log }_{2}n}{2}\rfloor$：同基本算法中提到的限制一样，INT8 GEMM 累加器是 INT32，一次内积累加 $n$ 项 $2^{2\beta }$ 量级的乘积，不能撑爆 $2^{31}$。 于是有：

     $$\beta =\min \left(7,\lfloor \frac{31-{\log }_{2}n}{2}\rfloor \right)$$

• 在确定完切片长度之后，我们需要知道应该将浮点数放置在多大间距的网格里面，也就是“保留到哪一位”：

  • 设残差行最大值指数为 $E=\lceil {\log }_2{\max }_j|a_{ij}|\rceil$，则残差的量级上界是 $2^E$（残差元素满足 $|a|\lesssim 2^E$）。

  • 我们假设取得的间距为 $\mu$，则通过 $(a+\sigma )-\sigma$ 提取出来的切片必然是 $\mu$ 的整数倍（它就是从 $\mu$ 间距网格上取的点），而切片除以 $\mu$ 就是想要的 INT8。

  • 我们希望切片除以 $\mu$ 后是个 $\beta$ 位带符号整数。由于 RN 可能向上舍入，因此有：

    $$\left|\frac{\mathrm{slice}}{\mu }\right|\lesssim \frac{2^E}{\mu }\stackrel{!}{\le }{2}^{\beta -1}$$

    解出：

    $$\mu \ge 2^{E-(\beta -1)}=2^{E-\beta +1}$$

    取等号（用满 $\beta$ 位）就得到：

    $$\mu =2^{E-\beta +1}=2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }$$

• 接着我们构建 magic constant，并实现切分 + RN：

  $$\sigma =0.75\cdot 2^{53}\cdot \mu ,(A_s{)}_{ij}=(a+\sigma )-\sigma$$

• 现在 $A_s$ 已经在网格上了，再除以 $\mu$ 、加上符号位就得到所需 INT8。

RN 切分算法

Input / Output 与伪代码

Input: $\mathbf{A}\in {\mathbb{F}}_{64}^{m\times n},k\in \mathbb{N}$

Output: ${\mathbf{A}}_s^{\prime }\in {\mathbb{I}}_8^{m\times n},{\mathbf{\mu }}^{(s)}\in {\mathbb{F}}_{64}^m,s=1,\dots ,k$

$$\begin{array}{rl}1.& \beta \leftarrow \min \left(7,\lfloor \frac{31-{\log }_{2}n}{2}\rfloor \right)\\ 2.& \text{for }s=1,\dots ,k\text{ do}\\ 3.& {\mu }_i^{(s)}\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }\text{for }i=1,\dots ,m\\ 4.& {\sigma }_i\leftarrow 0.75\cdot 2^{53}\cdot {\mu }_i^{(s)}\text{for }i=1,\dots ,m\\ 5.& (A_s{)}_{ij}\leftarrow (a_{ij}+{\sigma }_i)-{\sigma }_i\text{for }i=1,\dots ,m\{\text{Extract in round-to-nearest-even mode}\}\\ 6.& (A_s^{\prime }{)}_{ij}\leftarrow INT8\left(({\mu }_i^{(s)}{)}^{-1}\cdot (A_s{)}_{ij}\right)\forall i,j\{\text{Convert to INT8}\}\\ 7.& a_{ij}\leftarrow a_{ij}-(A_s{)}_{ij}\forall i,j\\ 8.& \text{end for}\end{array}$$

方向 2：Group-wise error-free accumulation

• 回顾之前的 profile 小节，FP64 accumulation 这一步占了总时间的 40%–50%，主要是因为对每个 $(s,g-s)$ 都要做：

  $$\mathbf{C}\leftarrow \mathbf{C}+\mathrm{diag}({\mathbf{\mu }}^{\prime \prime })2^{-\beta s+1-\beta (g-s)+1}FP64({\mathbf{C}}^{\prime \prime })\mathrm{diag}({\mathbf{\nu }}^{\prime \prime })$$

  这里面每一项都要：（1）把 INT32 结果 ${\mathbf{C}}^{\prime \prime }$ 转成 FP64（慢）；（2）乘两个对角缩放 $\mathrm{diag}({\mathbf{\mu }}^{\prime \prime })$、$\mathrm{diag}({\mathbf{\nu }}^{\prime \prime })$ 和 $2$ 的幂；（3）将 FP64 结果加到 $\mathbf{C}$ 上。而这样的项一共有 $\frac{k(k+1)}{2}$ 个。每一项都触发一次 INT32 $\to$ FP64 转换、缩放 + FP64 加法，这就是 $40\text{–}50\%$ 的来源。

• 但是可以观察到上面这些项里面，有很多项共享同一个幂缩放因子，因此可以先在 INT32 累加器里面加起来，再一次性转 FP64，而不是每项转一次。于是考虑分组：

  $$G_g:=\{(i,j)\mid i+j=g\}.$$

  也就是所有下标和等于 $g$ 的分片对归为一组。对同一组 $G_g$ 内的项，它们的缩放因子 $2^{-\beta i+1}\cdot 2^{-\beta j+1}=2^{-\beta (i+j)+2}=2^{-\beta g+2}$ 对全组来说是共用的。于是原累加和公式可以改写成：

  $$\sum _{s=1}^{g-1}\mathrm{diag}({\mathbf{\mu }}^{\prime \prime })2^{-\beta i_s+1}{2}^{-\beta j_s+1}{\mathbf{A}}_{i_s}^{\prime \prime }{\mathbf{B}}_{j_s}^{\prime \prime }\mathrm{diag}({\mathbf{\nu }}^{\prime \prime })=2^{-\beta g+2}\mathrm{diag}({\mathbf{\mu }}^{\prime \prime })\left(\sum _{s=1}^{g-1}{\mathbf{A}}_{i_s}^{\prime \prime }{\mathbf{B}}_{j_s}^{\prime \prime }\right)\mathrm{diag}({\mathbf{\nu }}^{\prime \prime }).$$

  此时括号里面的求和就是纯 INT32 累加，不用每次转 FP64。整组算完，才做一次 INT32 → FP64 转换、一次缩放和一次 FP64 加法。

• 但是由于 INT32 累加器只有 31 位数值，加太多项会溢出。所以容易求得，能塞进一个累加器的项数有上限 $r$：

  $$r:=\max \left(1,2^{31-2\beta -\lceil {\log }_{2}n\rceil }\right)$$

  组内前 $r$ 项在 INT32 累加器里能精确无误差地加；若一组超过 $r$ 项，就分批：每攒够 $r$ 项一起转 FP64。

Group-wise 累加算法

Input / Output 与伪代码

Input: ${\mathbf{A}}_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n},{\mathbf{\mu }}^{\prime \prime (s)}\in {\mathbb{F}}_{64}^m,{\mathbf{B}}_s^{\prime \prime }\in {\mathbb{I}}_8^{n\times p},{\mathbf{\nu }}^{\prime \prime (s)}\in {\mathbb{F}}_{64}^p,s=1,\dots ,k$

Output: $\mathbf{C}\in {\mathbb{F}}_{64}^{m\times p}$

$$\begin{array}{rl}1.& \beta \leftarrow \min \left(7,\lfloor \frac{31-{\log }_{2}n}{2}\rfloor \right)\\ 2.& r\leftarrow \max \left(1,2^{31-2\beta -\lceil {\log }_{2}n\rceil }\right)\{\text{\#addends for error-free accum.}\}\\ 3.& \mathbf{C}\leftarrow {\mathbf{O}}_{m,p}\\ 4.& \text{for }g=2,\dots ,k+1\text{ do}\\ 5.& q\leftarrow 0\\ 6.& {\mathbf{C}}^{\prime \prime }\leftarrow {\mathbf{O}}_{m,p}\{{\mathbf{C}}^{\prime \prime }\in {\mathbb{I}}_{32}^{m\times p}\}\\ 7.& \text{for }s=1,\dots ,g-1\text{ do}\\ 8.& q\leftarrow q+1\\ 9.& {\mathbf{C}}^{\prime \prime }\leftarrow {\mathbf{C}}^{\prime \prime }+{\mathbf{A}}_s^{\prime \prime }{\mathbf{B}}_{g-s}^{\prime \prime }\{\text{Compute using GEMM on INT8 Tensor Core}\}\\ 10.& \text{if }q==r\Vert s==g-1\text{ then}\\ 11.& \mathbf{C}\leftarrow \mathbf{C}+2^{-\beta g+2}\mathrm{diag}({\mathbf{\mu }}^{\prime \prime })FP64({\mathbf{C}}^{\prime \prime })\mathrm{diag}({\mathbf{\nu }}^{\prime \prime })\{\text{Compute in FP64}\}\\ 12.& q\leftarrow 0\\ 13.& {\mathbf{C}}^{\prime \prime }\leftarrow \mathbf{O}\\ 14.& \text{end if}\\ 15.& \text{end for}\\ 16.& \text{end for}\end{array}$$

方向 3：Hybrid

• 很自然想到，可以把方向 1 与方向 2 结合起来。但 group-wise 累加成立的前提是 $\mathrm{diag}({\mathbf{\mu }}^{\prime \prime })$ 和 $\mathrm{diag}({\mathbf{\nu }}^{\prime \prime })$ 在循环内恒定，才能作为公因子提到组外。bitmask split 天然满足，RN split 则不满足：

  $${\mu }_i^{(s)}\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }$$

  这里的 ${\mu }^{(s)}$ 取自当前残差，会逐片变化，因此不能直接使用分组加速。

• 解决思路：将“每片衰减”从 $\mu$ 转移到 $\sigma$：

  RN 原来的每片量级衰减通过重算 ${\mu }^{(s)}$ 实现。可以让 ${\mu }^{\prime \prime }$ 固定（在循环外根据原始 $A$ 计算一次），把每片缩小量级的工作转移到 magic constant $\sigma$，令其带有随 $s$ 变化的二次幂因子：

  $$\begin{array}{rl}{\mu }_i^{\prime \prime }& \leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }\text{for }i=1,\dots ,m\\ {\sigma }_i& \leftarrow 0.75\cdot 2^{53-\beta (s-1)}\cdot {\mu }_i^{\prime \prime }\end{array}$$

Hybrid 切分算法

Input / Output 与伪代码

Input: $\mathbf{A}\in {\mathbb{F}}_{64}^{m\times n},k\in \mathbb{N}$

Output: ${\mathbf{A}}_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n},{\mathbf{\mu }}^{\prime \prime }\in {\mathbb{F}}_{64}^m,s=1,\dots ,k$

$$\begin{array}{rl}1.& \beta \leftarrow \min \left(7,\lfloor \frac{31-{\log }_{2}n}{2}\rfloor \right)\\ 2.& {\mu }_i^{\prime \prime }\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }\text{for }i=1,\dots ,m\\ 3.& \text{for }s=1,\dots ,k\text{ do}\\ 4.& {\sigma }_i\leftarrow 0.75\cdot 2^{53-\beta (s-1)}\cdot {\mu }_i^{\prime \prime }\text{for }i=1,\dots ,m\\ 5.& (A_s{)}_{ij}\leftarrow (a_{ij}+{\sigma }_i)-{\sigma }_i\forall i,j\{\text{Extract in round-to-nearest-even mode}\}\\ 6.& (A_s^{\prime \prime }{)}_{ij}\leftarrow INT8\left(({\mu }_i^{\prime \prime (s)}{)}^{-1}\cdot 2^{-\beta (s-1)}\cdot (A_s{)}_{ij}\right)\forall i,j\{\text{Convert to INT8}\}\\ 7.& a_{ij}\leftarrow a_{ij}-(A_s{)}_{ij}\forall i,j\\ 8.& \text{end for}\end{array}$$

---

优化效果与性能对比

结论： EF（group-wise error free）和 Hybrid 版本的 ozIMMU 在 FP64 accumulation 上花费的时间明显少于 RN 版本。

结论： 切片数增加会降低吞吐；EF 和 Hybrid 通过减少 FP64 开销，在多数测试规模上相对基础 ozIMMU 获得了稳定加速。

• 这张图左半是吞吐绝对值，右半是相对 ozIMMU 的加速比。
• 左半张：
  • 每条曲线是一个方法在某个矩阵规模 $n$ 下的吞吐。核心趋势：切片数越多，吞吐越低。因为切片越多，要算的 INT8 小 GEMM 数量按 $k(k+1)/2$ 增长，总计算量上去了，吞吐也就掉了。
  • 图中还有一条水平的 FP64 基准线。只有还没跌破基准线的切片数是值得的优化。
  • 同色系里 $n$ 越大的曲线整体越高：大矩阵下 INT8 Tensor Core 利用率更充分，低精度模拟的优势才发挥得出来。
• 右半张：
  • 这是四个改进变体相对基础版 ozIMMU 的加速比。EF 和 H 的曲线基本都在 1.2 以上，说明 group-wise 累加确实减少了 FP64 开销、带来净加速。

结论： 在相近吞吐下，Hybrid 和 EF 通常能取得更低误差，把精度—性能折中曲线推向更理想的左上方。

• 这张图中横轴是最大相对误差（越往左越准），纵轴是吞吐 TFLOPS（越往上越快），每个点上标的数字是该方法使用的切片数。固定 $n=4096$。
• 理想点在左上角（又快又准）。这张图直观地画出了精度-性能的 trade-off 曲线。
• 在同样吞吐下，H 和 EF 的点比 ozIMMU 更靠左（更准）。这说明改进方法能用更少或相当的代价拿到更高精度。

参考

• Ootomo, Ozaki, and Yokota, DGEMM on Integer Matrix Multiplication Unit.
• ozIMMU reference implementation.
• Uchino, Ozaki, and Imamura, Performance Enhancement of the Ozaki Scheme on Integer Matrix Multiplication Unit.
• RIKEN accelerator for ozIMMU.