Accelerating Multiple-Precision LU Decomposition with Ozaki Scheme II
Abstract
Solving ill-conditioned linear systems needs LU decomposition in precisions beyond binary64. The standard approaches -- GMP/MPFR, or multi-component arithmetic such as double-double -- leave every scalar multiply-add inside the O(n^3) update in multiple precision, and so cannot exploit the low-precision matrix engines that dominate current hardware. We build a blocked, partially pivoted LU decomposition on Ozaki scheme II, which replaces the multiple-precision GEMM by exact integer modular products followed by an explicit CRT reconstruction, and implement it for multi-component (DD/TD/QD) and arbitrary precision, on CPUs and GPUs. Two contributions make this practical: a direct conversion between the non-overlapping expansion format and the internal fixed-point representation, which removes the MPFR round trip and is worth a factor of 2.4--4.7; and new FP16, FP8 and binary64 GPU back-ends, the FP8 one using a balanced base-17 two-digit encoding that keeps the full 362.8-bit CRT capacity of INT8. On an Arm/GB10 and an x86/H100 the fastest back-end changes with the machine: INT8 wins on GB10, whereas on H100 binary64 is fastest for QD, beating both the native implementation and INT8. The essential point of Ozaki scheme II is thus not to use a low-precision engine but to choose the format maximising bits-per-modulus times engine throughput. For the Lotkin matrix at p ~ 1.2 log2 cond(A) we reach relative errors of 10^-645 at n=2048, up to 2.84x faster than a fully OpenMP-parallel multiple-precision LU. We also show by measurement that the O(p) advantage of scheme II over scheme I applies only to the GEMM term: modular reduction and CRT grow as O(p^2) and dominate the runtime at the matrix sizes considered here. The implementation is released as open source.
Related
- Ozaki Scheme II Is Fast on CPUs Too: Multiple-Precision Matrix Multiplication on Intel AMX-INT8 and Arm SVE2-i8mm
- Emulating DGEMM and ZGEMM with INT8 Matrix Engines on Intel Arc GPUs: Ozaki Scheme II on the Xe2 Architecture
- Accelerating phase-field simulations of grain growth via tensor decomposition
- Linear Robin-type Domain Decomposition scheme for Mixed Formulation of Richards' Equation
- Mixed Precision Compression of Tucker Decomposition
BibTeX
@misc{kouya2026accelerating,
title = {{Accelerating Multiple-Precision LU Decomposition with Ozaki Scheme II}},
author = {Tomonori Kouya},
howpublished = {arXiv},
year = {2026},
month = sep,
doi = {10.48550/arxiv.2609.29554},
eprint = {2609.29554},
archivePrefix = {arXiv},
url = {https://arxiv.org/abs/2609.29554},
note = {CC BY 4.0},
}
RIS
TY - UNPB
TI - Accelerating Multiple-Precision LU Decomposition with Ozaki Scheme II
AU - Tomonori Kouya
PY - 2026
DA - 2026/09/01
PB - arXiv
DO - 10.48550/arxiv.2609.29554
UR - https://arxiv.org/abs/2609.29554
AB - Solving ill-conditioned linear systems needs LU decomposition in precisions beyond binary64. The standard approaches -- GMP/MPFR, or multi-component arithmetic such as double-double -- leave every scalar multiply-add inside the O(n^3) update in multiple precision, and so cannot exploit the low-precision matrix engines that dominate current hardware. We build a blocked, partially pivoted LU decomposition on Ozaki scheme II, which replaces the multiple-precision GEMM by exact integer modular products followed by an explicit CRT reconstruction, and implement it for multi-component (DD/TD/QD) and arbitrary precision, on CPUs and GPUs. Two contributions make this practical: a direct conversion between the non-overlapping expansion format and the internal fixed-point representation, which removes the MPFR round trip and is worth a factor of 2.4--4.7; and new FP16, FP8 and binary64 GPU back-ends, the FP8 one using a balanced base-17 two-digit encoding that keeps the full 362.8-bit CRT capacity of INT8. On an Arm/GB10 and an x86/H100 the fastest back-end changes with the machine: INT8 wins on GB10, whereas on H100 binary64 is fastest for QD, beating both the native implementation and INT8. The essential point of Ozaki scheme II is thus not to use a low-precision engine but to choose the format maximising bits-per-modulus times engine throughput. For the Lotkin matrix at p ~ 1.2 log2 cond(A) we reach relative errors of 10^-645 at n=2048, up to 2.84x faster than a fully OpenMP-parallel multiple-precision LU. We also show by measurement that the O(p) advantage of scheme II over scheme I applies only to the GEMM term: modular reduction and CRT grow as O(p^2) and dominate the runtime at the matrix sizes considered here. The implementation is released as open source.
ER -