Analysis of Floating-Point Matrix Multiplication Computed via Integer Arithmetic

• Reduced-precision formats can exceed binary64 throughput by over two orders of magnitude, motivating mixed-precision reformulations of scientific computations.Source: §1 · p. 2
• The Ozaki scheme converts rows of A and columns of B to block fixed-point form, splits them into exactly multipliable slices, and accumulates the partial products in floating-point arithmetic.Source: §1 · p. 2
• The paper’s error analysis shows that badly scaled rows of A or columns of B can make the scheme inaccurate, while unequal slice counts can reduce multiplications under an accuracy bound.Source: §1 · p. 3
• Experiments confirm that badly scaled matrices may require extremely many slices, whereas judiciously unequal slicing can preserve accuracy in selected cases.Source: §1 · p. 3
• GPU benchmarks compare high-performance Ozaki implementations, including the first reported results on an NVIDIA Blackwell GPU, and show reduced-product variants working under favorable scaling.Source: §1 · p. 3

• A floating-point number uses a sign s ∈ {0, 1}, exponent e ∈ Z, and significand m ∈ [1, 2); limiting significand precision and exponent range discretizes the reals.Source: §2 · p. 3
• The floating-point system F⟨emin, emax, p⟩ limits exponent range and significand precision, with binary64 using p = 53 and emax = 1023.Source: §2 · pp. 3–4
• Round-to-nearest and the standard model assume elementary floating-point operations incur relative error bounded by the unit roundoff u = 2^-p.Source: §2 · p. 4
• Signed integer format It represents values in [−2^t, 2^t − 1], and addition and multiplication are exact unless overflow occurs.Source: §2 · p. 4
• Integer MMA units compute AB+C for signed-integer matrices, using INT8 inputs and 32-bit signed accumulators on NVIDIA Hopper and Blackwell tensor cores.Source: §2 · p. 4

The Ozaki scheme

• The integer Ozaki scheme converts each row of A and column of B to block fixed-point form, then splits these representations into integer slices for MMA computation.Source: §3 · pp. 5–6
• Each slice product A(ℓ)B(h) is exact when computation and accumulation avoid overflow, after which the products are converted to floating point and accumulated.Source: §3 · p. 7
• Products with ℓ+h > s+1 are omitted because their scaling factor 2^−(ℓ+h)t makes their contributions small, reducing the number of matrix multiplications.Source: §3 · p. 7

Integer accumulation

• The accumulation format needs 2t + ⌈log2 k⌉ bits, and the slice width must respect both the MMA input width and the accumulator capacity.Source: §3 · pp. 7–8

Enhanced algorithm

• For hardware with fast integer MMA units, accumulating matrix products in integer arithmetic when ℓ+h is constant reduces floating-point sums per result element from s(s + 1)/2 to s − 1.This increases integer sums from k to k + s − 1 and may cause overflow if t is computed as in algorithm 3.1.Source: §3 · pp. 9–10
• Round-to-nearest can improve accuracy when the number of slices is fixed, but it costs more to implement and requires the first slice to contain t − 1 bits.Source: §3 · p. 10

4. Error analysis.

• The total error decomposes into truncation from slicing, rounding from floating-point accumulation, and their interaction.The analysis treats these sources separately before combining them into matrix-multiplication error bounds.Source: §4 · p. 10, §4.2 · p. 12

4.1. General results on truncation error.

• Allowing different slice counts sA and sB preserves exact-product representations only when the counts meet the minimum requirements, while smaller counts introduce bounded conversion errors.The cost depends directly on the slice counts, so reducing them can improve practicality at the expense of truncation accuracy.Source: §4.1 · p. 11
• Badly scaled rows of A or columns of B can produce large relative conversion errors, even when many slices are used.The bounds become large when entries vary widely in magnitude within a row or column.Source: §4.1 · p. 12

4.2. Computation of all products.

• Integer partial products are exact when the accumulator is wide enough, but converting and summing them in floating-point arithmetic introduces a separate accumulation error.The bounds use the number of summed matrices and the unit roundoff of the floating-point format.Source: §4.2 · pp. 12–13
• The error bound can be substantial for large scaling factors κA or κB, while increasing sA and sB improves accuracy but increases the χ(sA, sB) matrix multiplications.The slice counts therefore control a direct performance-accuracy trade-off.Source: §4.2 · p. 13, §4.5 · p. 17

4.3. Reduction in number of products.

• Reducing the number of products changes the error bound asymmetrically: when κA is large, the reduced-product bound can grow faster than the full-product bound.The analogous dependence on κB appears when the roles of A and B are switched.Source: §4.3 · pp. 13–15

4.4. Accumulation in floating-point arithmetic.

• The accumulation analysis divides final-accumulation sums into levels computed exactly internally, so rounding errors occur only when values from different levels are added together.The bound is tight as long as max{sA, sB} ≤255.Source: §4.4 · pp. 15, 17

4.5. Discussion.

• The proposed practical choice evaluates small integer pairs (sA, sB) and selects the pair minimizing χ(sA, sB) among those satisfying the target error bound.The quantities κA and κB can be computed while forming the scale vectors, making this selection inexpensive.Source: §4.5 · p. 17

5.1. Behaviour on badly scaled matrices.

• More severe vector scaling increases the slices needed for binary64 accuracy: about 7 slices suffice at φ = 0, whereas over 20 are needed at φ = 100.The two tested Ozaki variants perform similarly on this example.Source: §5.1 · p. 18
• For φ = 13, even the 10-slice variant of the Ozaki scheme produces large errors, while using more slices can improve accuracy but increase runtime and make it impractical.More slices reduce truncation error in the inputs; impracticality depends on the relative performance of the available integer and binary64 arithmetic.Source: §5.1 · p. 18

5.2. Solving Ax = b with block LU factorization.

• In block LU factorization, all tested matrices meet the HPL threshold with eight slices for both factors, while some matrices meet it with asymmetric slice counts.For the Anymatrix core/blockhouse case, the error is below the HPL threshold when both matrices use at least seven slices.Source: §5.2 · p. 21, Figure 5.4 · p. 21

5.3. Benchmarking on GPUs.

• On GPUs, three slices deliver the best cuIMMA performance but yield roughly 2.4× and 7.6× speedups with forward error of order 10^-10 versus cublas’s 10^-19; seven slices restore comparable accuracy.On B200, cuIMMA can outperform binary64 cublas while maintaining similar accuracy because the INT8-to-binary64 tensor-core peak-performance ratio is 112.5; higher-level MAGMA algorithms generally require eight slices.Source: §5.3 · pp. 22–24
• For badly scaled GPU inputs, the safe random-matrix choice s = 8 is insufficient, and κD values of 10^30 or 10^40 remain too inaccurate even with 18 slices.Lower κD values can reach binary64 accuracy only by increasing slices, severely reducing performance.Source: §5.3 · p. 25

6. Conclusion.

• Badly scaled input matrices can make the emulation highly inaccurate even when many slices are used, preventing it from fully replacing binary64 arithmetic.The current scheme also does not correctly handle negative zeros, infinities, or NaNs.Source: §6 · p. 26