---
title: Ozaki Scheme I for High-Precision GEMM
url: https://www.emergentmind.com/topics/ozaki-scheme-i
type: topic
---

# Ozaki Scheme I for High-Precision GEMM

Ozaki Scheme I is a decomposition-based algorithm for computing high-precision matrix products using only multiple calls to highly optimized, lower-precision matrix multiplication routines. By slicing each operand into lower-precision “digits” or “slices,” the scheme expresses the exact or nearly exact high-precision result as a sum of low-precision cross-products, allowing the aggregation of modern hardware’s drastic performance advantage for small-precision GEMM units to solve high-precision scientific problems efficiently.

## 1. Principle and Mathematical Formulation

Given matrices $A \in \mathbb{R}^{m \times k}$ and $B \in \mathbb{R}^{k \times n}$, both stored in a high (long) precision (mantissa $L$ bits), the core idea is to decompose each entry into sums of $S$-bit slices:
\[
A = \sum_{i=1}^D A^{(i)}, \qquad B = \sum_{j=1}^D B^{(j)},
\]
with each $A^{(i)}, B^{(j)}$ exactly representable in a short-precision format (mantissa $S \ll L$).

The matrix product is then evaluated as
\[
C = A B = \sum_{i=1}^D \sum_{j=1}^D A^{(i)} B^{(j)} = \sum_{i,j} {\rm GEMM}_S(A^{(i)}, B^{(j)}),
\]
where each ${\rm GEMM}_S$ is performed in $S$-bit precision, typically with highly efficient vendor BLAS or hardware tensor core routines. The accumulation into $C$ is performed in high-precision arithmetic (exact in $L$ bits or higher).

The parameter $D$ is the number of slices per operand. For full-precision results, it must satisfy
\[
D \ge \left\lceil \frac{L + \log_2 k}{S} \right\rceil,
\]
ensuring that products of slices do not overlap and each product is computed without rounding error in $S$-bit arithmetic [2301.09960].

## 2. Decomposition Algorithm (Slicing)

Any high-precision scalar or matrix entry $x$ is recursively split:
- $R^{(0)} \gets x$
- For $i=1, \ldots, D-1$:
  - $x^{(i)} \gets$ round or truncate $R^{(i-1)}$ to $S$ bits (short-precision)
  - $R^{(i)} \gets R^{(i-1)} - x^{(i)}$ (computed in high precision)
- $x^{(D)} \gets R^{(D-1)}$

This process is repeated for all entries. The blockwise (matrix) version applies the same logic row- or column-wise.

For hardware integer units (e.g., INT8 tensor cores), slices correspond to blocks of bits in the mantissa, often with a scaling factor per row/column to ensure exponent alignment. The decomposition can also leverage two-level blocking for memory efficiency [2508.00441, 2306.11975, 2606.25453].

## 3. Computational Workflow and Pseudocode

The generic Ozaki Scheme I pseudocode is:

```python
def OzakiGEMM(A, B, D):
    # Step 1: Decompose matrices into D slices
    A_slices = Slice(A, D)  # Each in S bits
    B_slices = Slice(B, D)
    # Step 2: Multiply and accumulate
    C = zeros(shape=(A.shape[0], B.shape[1]), dtype=high_precision_type)
    for i in range(D):
        for j in range(D):
            T = GEMM_S(A_slices[i], B_slices[j])  # Low-precision GEMM
            C += cast_to_high_precision(T)
    return C
```
[2301.09960, 2511.13778].

In fused-kernel (GPU) implementations, all cross-products and their accumulation are performed in a single kernel, eliminating off-chip round-trips for partial sums and maximizing arithmetic intensity [2606.25453].

## 4. Error Analysis and Parameter Choice

The approximation error comprises:
- **Truncation error** from residual slices, bounded by $2^{-S D}\|A\|\|B\|$.
- **Low-precision GEMM error** from each slice product; if the cross-product is exact (i.e., no rounding because of sufficiently small inputs and large enough accumulator/range), this is zero.
- **Accumulation error** from summing $D^2$ terms in high precision.

The forward error satisfies:
\[
\left\| C_{\text{Ozaki}} - AB \right\| \le O(2^{-SD}) \|A\|\|B\|
\]
Choosing $D$ as above ensures the error is within the target precision's machine epsilon.

For hardware with limited accumulator precision (e.g., INT32 in INT8 x INT8 → INT32 TCs), the slice width is constrained by accumulator size, and the number of slices $D$ must respect overflow bounds [2606.25453, 2306.11975]. Exponent-span-based estimators further refine $D$ for diverse input distributions [2511.13778].

## 5. Hardware Mapping, Extensions, and Implementational Variants

| Scenario                  | Slicing Format   | Kernel Type            | Accumulation | Notes                              |
|---------------------------|------------------|------------------------|--------------|-------------------------------------|
| FP64 → FP16/FP8 gemms     | Float (e.g. 16b) | cuBLAS/TC FP16/FP8     | FP64         | Biased for hardware throughput      |
| FP64 → INT8 TCs           | Integer (8b)     | IMMA (INT8×INT8→INT32) | FP64         | Exponent alignment for range        |
| MPFR arbitrary precision  | Float (e.g. 64b) | High-precision SW BLAS | MPFR         | Arbitrarily many slices             |
| Fused GPU kernel          | INT8/FP8 (8b)    | Persistent/block kernel | FP32/FP64    | On-chip accumulation, optimized for minimal memory traffic [2606.25453] |

*Slice decomposition* employs scaling strategies for range, and unsigned-integer slice encoding further reduces overhead [2511.13778]. On modern GPU hardware, the triangular product schedule and in-register accumulation enable high utilization, with sustained performance up to 80%–90% of INT8 peak [2606.25453]. FPGA, CPU (MKL), and custom quantum circuit frameworks employ similar schemes [2306.11975].

## 6. Complexity and Performance

Let $n$ denote matrix dimension, $L$ the high-precision mantissa bits, $S$ the slice precision, and $\alpha$ the speedup factor of $S$-bit GEMM over $L$-bit GEMM. Then:

- Classical high-precision GEMM: $O(n^3)$ $L$-bit operations
- Strassen: $O(n^{\log_2 7}) \approx O(n^{2.807})$
- Ozaki I: $D^2$ calls to $O(n^3)$ $S$-bit GEMMs, plus $O(D n^2)$ slicing/accumulation

The total leading cost is $O((L/S)^2 n^3) / \alpha$ for square matrices. For moderate $L/S$ (e.g., $3~\text{to}~8$), and $n$ in the low thousands, Ozaki I frequently outperforms high-precision Strassen and direct software GEMMs, sometimes by factors exceeding $5\times$ [2301.09960, 2307.06072].

On GPUs with fused persistent-kernel implementations, Scheme I achieves throughput of up to $1,639$ Top/s (Hopper, $83\%$ INT8 peak) and $3,654$ Top/s (Blackwell, $81\%$) [2606.25453], with up to $4.1\times$ (FP8TC) acceleration over hardware FP64 [2508.00441].

## 7. Applicability, Extensions, and Limitations

Ozaki Scheme I generalizes to:
- **Complex matrices**: Admits the 3M/4M strategies, combining slicing with minimal real GEMMs [2307.06072].
- **LU decomposition and other BLAS**: Enables efficient batched blocks for high-precision panel updates [2307.06072].
- **FFT emulation**: Underpins high-precision FFT by splitting input vectors/convolutions and summing low-precision transforms, including integer-based NTT/CRT for error-free accumulation [2603.29129].

Limitations:
- For very high bit-widths ($L/S \gg 10$) or extremely large $n$, $D^2$ scaling becomes prohibitive, so Strassen or direct approaches dominate.
- Fused-kernel memory footprint grows with $D$; on-chip buffer budgets may force kernel tiling for large $D$ on current GPU architectures [2606.25453].
- For small input sizes, decomposition and pre/post-processing overheads can exceed the GEMM cost.
- The method assumes reliable higher-precision accumulation on the host hardware; if not available, accuracy is not guaranteed.

Practical variants integrate automatic slice tuning (via exponent-span estimators), unsigned/signed residue tracking, and host-transparent fallback policies to ensure correctness across diverse input matrices [2511.13778].

## References

- [2301.09960] Acceleration of Multiple Precision Matrix Multiplication using Ozaki scheme
- [2307.06072] Acceleration of complex matrix multiplication using arbitrary precision floating-point arithmetic
- [2508.00441] DGEMM without FP64 Arithmetic -- using FP64 Emulation and FP8 Tensor Cores with Ozaki Scheme
- [2306.11975] DGEMM on Integer Matrix Multiplication Unit
- [2606.25453] EmuGEMM: Fused Tensor Core Kernels for Precision Emulation in Matrix Multiplication
- [2511.13778] Guaranteed DGEMM Accuracy While Using Reduced Precision Tensor Cores Through Extensions of the Ozaki Scheme
- [2603.29129] Computing FFTs at Target Precision Using Lower-Precision FFTs
- [2602.19090] Forward Error-Oriented Iterative Refinement for Eigenvectors of a Real Symmetric Matrix

Source: https://www.emergentmind.com/topics/ozaki-scheme-i