Ozaki-I Scheme: Exact Low-Precision GEMM

• Ozaki-I Scheme is an error-free transformation that decomposes high-precision GEMM into a structured sum of low-precision GEMMs using mantissa slicing.
• It employs a triangular schedule where non-overlapping slices are computed via optimized INT8 kernels and then accurately recombined for high-precision results.
• The approach balances precision and performance, with extensions for hardware-specific optimizations like Tensor Cores and FP8 pathways outperforming native high-precision GEMM.

Ozaki-I Scheme is an error-free transformation for matrix multiplication that reconstructs a target-precision product from a structured sum of lower-precision products. In the numerical linear algebra literature, it is used to compute a high-precision GEMM $C = AB$ by decomposing $A$ and $B$ into non-overlapping slices, evaluating a triangular set of slice-pair products with optimized low-precision kernels, and recombining them in higher precision (Lu et al., 24 Jun 2026). In modern implementations, the scheme is especially associated with integer or low-precision Tensor Core pathways, where exact INT8$\times$INT8$\to$INT32 accumulation or carefully bounded FP32 accumulation can be exploited to emulate FP32- or FP64-level results at substantially higher throughput than native high-precision arithmetic (Ootomo et al., 2023).

1. Core definition and nomenclature

The Ozaki scheme was introduced as an error-free transformation of dot products extended to matrix multiplication. Its central objective is to rewrite one high-precision GEMM into a finite sum of low-precision GEMMs whose individual products and blockwise sums are performed without rounding error, so that the final reconstruction achieves the accuracy of the original high precision (Mukunoki, 1 Aug 2025). In this setting, Ozaki-I denotes the slice-based formulation in which operands are decomposed into mantissa segments and only a triangular subset of slice products is needed.

A persistent terminological feature of the literature is that “Ozaki-I” is used in two closely related ways. Some papers use “Ozaki scheme I” for the conventional slice-based formulation in general, including floating-point realizations (Ozaki et al., 10 Apr 2025). Other papers use “Ozaki-I” more specifically for the integer-based variant adapted to integer matrix multiplication units, where fixed-point arithmetic with shared exponents per slice makes each low-precision integer GEMM exact and leaves only the final reconstruction to floating-point rounding (Ootomo et al., 2023). The common structure is unchanged: shared-place mantissa splitting, exact or rigorously controlled low-precision subproducts, and higher-precision accumulation.

This scheme differs from elementwise-place multiword methods such as double-double arithmetic. In the IMMU-oriented formulation, it is a shared-place mantissa-splitting algorithm: rows of $A$ and columns of $B$ are aligned by shared exponents, sliced into low-bit components, and then recombined after the low-precision GEMMs have been evaluated (Ootomo et al., 2023). That shared-place organization is what makes the method naturally GEMM-oriented and compatible with vendor-optimized BLAS or Tensor Core kernels.

2. Mathematical structure

A standard formulation decomposes the operands into $\beta$-bit mantissa slices with per-row and per-column power-of-two scales:

$$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$

where $A'_i$ and $A$0 are INT8 slices and $A$1 are powers of two chosen per row and column to align exponents (Lu et al., 24 Jun 2026). The corresponding product expands as

$A$2

The defining Ozaki-I observation is that all pairs with the same index sum $A$3 have the same positional weight. This yields the triangular accumulation

$A$4

followed by the reconstruction

$A$5

The result is a triangular schedule with $A$6 exact INT8 GEMMs and $A$7 accumulators, rather than a full $A$8-term materialization (Lu et al., 24 Jun 2026).

The same structure is often described in digit language. With base $A$9, one conceptually represents a scalar by a $B$0-bit segmentation

$B$1

while row- and column-wise power-of-two scaling ensures that each slice fits the chosen low-precision format (Lu et al., 24 Jun 2026). In the IMMU literature, this becomes block floating-point with shared exponents per row of $B$2 and per column of $B$3, and the slice width is chosen from the accumulator precision and inner-product length. A common formula is

$B$4

where $B$5 is the accumulator mantissa length and $B$6 is the input mantissa length (Ootomo et al., 2023).

Exactness hinges on three properties. First, the integer MMA path is exact: INT8$B$7INT8$B$8INT32 multiplication and accumulation introduce no rounding in the integer core (Lu et al., 24 Jun 2026). Second, the weights are powers of two, so the shift-reduce applies exact binary scaling. Third, correctness assumes that the slices fit their ranges and that the accumulator does not overflow; several papers note that explicit closed-form INT32 overflow bounds are not given for all implementations, although practical scaling rules are designed to stay within range (Lu et al., 24 Jun 2026).

3. Triangular schedule and algorithmic workflow

At the algorithmic level, Ozaki-I computes only the leading diagonals of the slice-slice product lattice. In a conventional presentation, if

$B$9

then the exact identity is

$\times$0

and the subproducts $\times$1 with $\times$2 are the error-free part computed by BLAS kernels (Ozaki et al., 10 Apr 2025). This is the matrix-level counterpart of the triangular accumulation used in current Tensor Core implementations.

In the integer formulation, the splitting routine first computes shared exponents by row-wise and column-wise maxima, then extracts non-overlapping mantissa blocks into signed INT8 slices while preserving sign. The compute phase iterates over admissible slice pairs, executes an integer GEMM for each pair, rescales the INT32 result by a power-of-two weight and the outer product of the row and column exponents, and accumulates in FP64 (Ootomo et al., 2023). The dominant arithmetic count is the triangular number $\times$3, which explains both the controllable precision increase and the quadratic growth in the number of low-precision GEMMs.

A modern fused realization preserves the same mathematics but changes the dataflow. EmuGEMM interleaves the $\times$4 slices of each operand along the contraction dimension $\times$5 at $\times$6 granularity, so a single TMA load fetches all slices for a $\times$7-step. For each $\times$8-step, the kernel issues all pairings $\times$9 for $\to$0, keeps the $\to$1 INT32 accumulators on chip, and performs the weighted shift-reduce epilogue before storing to global memory (Lu et al., 24 Jun 2026). This eliminates the repeated global-memory materialization of intermediate INT32 results that dominated earlier realizations.

4. Hardware realizations and performance evolution

The recent literature on Ozaki-I is largely the literature on how to map the triangular schedule onto increasingly specialized low-precision hardware. The main progression runs from IMMU-based slice GEMMs, through reduced-accumulation variants, to fused persistent kernels and automatic runtime selection (Ootomo et al., 2023, Uchino et al., 2024, Lu et al., 24 Jun 2026, Schwarz et al., 16 Nov 2025).

Implementation	Mechanism	Reported result
Ootomo et al.	INT8 Tensor Cores with FP64 reconstruction	up to $\to$2 over cuBLAS DGEMM on consumer GPUs
ozIMMU_EF / ozIMMU_H	group-wise exact INT32 accumulation; reduced FP64 additions	$\to$3–$\to$4 over prior ozIMMU at equal $\to$5
EmuGEMM-I	fused integer Tensor Core kernels; on-chip triangular accumulation	up to $\to$6 Top/s on Hopper and $\to$7 Top/s on Blackwell
ADP with ESC	unsigned slicing; runtime slice selection and fallback	less than $\to$8 run time overhead; up to $\to$9 and $A$0 over native FP64 GEMM

The first generation of IMMU realizations established that integer Tensor Cores can be used for FP64-accurate DGEMM when the mantissa space is covered by enough slices and the final accumulation is done in FP64. On NVIDIA consumer GPUs, Ozaki-I outperformed cuBLAS DGEMM by up to $A$1, and the same machinery accelerated a quantum circuit simulation by up to $A$2 while maintaining FP64 accuracy (Ootomo et al., 2023).

The next major optimization targeted the reduction phase. Since FP64 accumulation could consume $A$3–$A$4 of runtime when INT8 matmuls were much faster than FP64 arithmetic, group-wise exact INT32 accumulation was introduced: several products along one diagonal are summed in INT32 before a single FP64 conversion and rescaling. This reduces the number of FP64 accumulations from $A$5 to a smaller count $A$6, with observed $A$7–$A$8 speedups over prior ozIMMU at equal $A$9 (Uchino et al., 2024).

EmuGEMM represents a more radical kernel-level fusion. Using Scheme I, it sustains up to $B$0 Top/s on NVIDIA GH200 Hopper ($B$1 of INT8 peak) and $B$2 Top/s on NVIDIA B200 Blackwell ($B$3). For large matrices, it surpasses cuBLAS TF32 throughput by up to $B$4 on Hopper and $B$5 on Blackwell at comparable accuracy (Lu et al., 24 Jun 2026). The critical change is not the algebraic formulation but the elimination of redundant memory round-trips: operand slices are loaded once per $B$6-step, and intermediate INT32 tiles never leave on-chip storage.

A separate line of work adapts the scheme to hardware that emphasizes FP8 rather than INT8. On Blackwell-class GPUs, FP8-based Ozaki with FP32 accumulation is presented as the better match when FP8 Tensor Cores have much higher throughput and software support is robust, while integer Ozaki-I remains advantageous on architectures emphasizing INT8/INT32 matrix engines (Mukunoki, 1 Aug 2025). Runtime systems such as ADP add automatic dynamic precision, Exponent Span Capacity, exception handling, and seamless fallback to native FP64, with all guardrails kept GPU-resident (Schwarz et al., 16 Nov 2025).

5. Accuracy, parameter selection, and limitations

Ozaki-I accuracy is governed primarily by slice width, slice count, and accumulation strategy. In the integer scheme, each additional slice adds approximately $B$7 bits of mantissa precision, and the use of a signed leading slice with unsigned remaining slices raises the usable $B$8 to $B$9 rather than $\beta$0 (Lu et al., 24 Jun 2026). In practice, EmuGEMM reports that $\beta$1 is suitable for TF32-like accuracy, $\beta$2–$\beta$3 for FP32-like accuracy, and that approaching FP64 requires larger $\beta$4 or Scheme II (Lu et al., 24 Jun 2026).

A more explicit IMMU-oriented parameterization starts from

$\beta$5

which ensures that the worst-case inner product in an integer GEMM cannot overflow the INT32 accumulator, and then sets $\beta$6 (Ootomo et al., 2023). This immediately couples achievable precision to the contraction length $\beta$7: larger $\beta$8 reduces $\beta$9, so more slices are needed for the same target accuracy.

Error is not dominated by the low-precision matmuls themselves when the exactness conditions hold. In EmuGEMM, errors stem only from finite-slice mantissa truncation and the final floating-point shift-reduce, and the paper explicitly states that it does not present a new closed-form error bound, instead referring to Ozaki et al. (2012) for a priori estimates (Lu et al., 24 Jun 2026). The 2024 IMMU study makes the truncation-accumulation tradeoff explicit:

$$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$0

while the accumulation error is bounded either by $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$1 in the baseline formulation or by $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$2 with group-wise accumulation (Uchino et al., 2024). These bounds formalize a central practical fact: adding slices reduces truncation error, but the resulting growth in the number of partial products can increase accumulation cost and memory traffic.

Several limitations recur across implementations. Wide exponent distributions increase the mantissa space length and require larger slice counts; insufficient slices reduce accuracy (Ootomo et al., 2023). Extremely large dynamic ranges may require more slices or a different scheme (Lu et al., 24 Jun 2026). Small matrices can be unfavorable because preprocessing and reconstruction overhead dominate (Lu et al., 24 Jun 2026). Memory pressure is also intrinsic to the unfused formulation, since storing many slices and partial products can be expensive; this is one reason fused kernels and on-chip accumulation materially change performance (Ozaki et al., 10 Apr 2025).

6. Extensions, contrasts, and related uses

Ozaki-I is the foundation for a broader family of emulation strategies, but it is not the only one. Scheme II replaces the triangular slice lattice with a modular method based on symmetric modulo arithmetic and the Chinese Remainder Theorem. In GEMM-oriented comparisons, Scheme II offers linear control over the number of matrix multiplications, whereas Scheme I requires $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$3 GEMMs; consequently, Scheme II is presented as preferable beyond FP32 toward FP64, while Scheme I is best up to FP32-like accuracy and uses lower workspace at matched $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$4 (Ozaki et al., 10 Apr 2025, Lu et al., 24 Jun 2026).

The same slicing idea has also been transported beyond GEMM. In FFT computation, Ozaki-I is applied to the cyclic convolution inside Bluestein’s algorithm: floating-point inputs are split into non-overlapping components, scaled to integers, convolved exactly by NTTs over 32-bit prime fields, and reconstructed by CRT. The resulting method computes target-precision FFTs using lower-precision backends and shows reduced relative error compared with FFTW- and Triple-Single-based alternatives, albeit at much higher runtime (Kawakami et al., 31 Mar 2026).

In eigensolvers, the Ozaki scheme has been adapted to forward error-oriented iterative refinement for eigenvectors of a real symmetric matrix. There, the accurate product $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$5 is formed by splitting only $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$6 deeply and keeping only the leading slice of the current eigenvector iterate, which reduces the work per iteration from $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$7 to $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$8 while preserving the role of exact block products in the refinement analysis (Terao et al., 22 Feb 2026). In arbitrary-precision complex GEMM, the same slice-and-recompose principle is combined with 3M and 4M complex multiplication, with the 3M-based Ozaki complex GEMM reported as $$A \approx \operatorname{diag}(\mu)\sum_{i=0}^{p-1} 2^{-\beta i} A'_i,\qquad B \approx \sum_{j=0}^{p-1} 2^{-\beta j} B'_j \operatorname{diag}(\nu),$$9–$A'_i$0 faster than its 4M-based counterpart at the same division count (Kouya, 2023).

A final terminological note is necessary because the label “Ozaki-I” is not unique across mathematics. In geometric function theory, “Ozaki-I” may refer to the class

$A'_i$1

an unrelated usage concerning close-to-convex functions rather than numerical linear algebra (Kumar et al., 2023). In computational mathematics and high-performance GEMM, however, Ozaki-I Scheme denotes the slice-based emulation framework centered on non-overlapping decomposition, triangular low-precision subproducts, and high-precision reconstruction.

Within that numerical meaning, Ozaki-I remains a canonical example of how algorithmic error-free transformation can be aligned with hardware specialization. Its enduring features are shared-place mantissa splitting, exact or carry-free low-precision inner kernels, and a reconstruction stage whose cost and stability can be traded against slice count, memory movement, and target precision.