Mixed-precision Interpolative Decomposition on GPUs

Low-rank approximation is a powerful tool for reducing storage and arithmetic complexity while preserving linearity. Traditional rank-revealing factorizations, such as the singular value decomposition (SVD) and pivoted QR, have been adapted to GPUs by increasing reliance on GEMM kernels, yet they remain throughput-limited. We present hyacin, a GPU implementation for constructing low-rank approximations via interpolative decomposition. hyacin combines column-pivoted Cholesky QR with extended integer quantization, which is an abstraction from the Ozaki-scheme error-free floating-point transformation, as well as mixed floating-point precisions. Integer-emulated GEMMs accelerate Gram matrix formation, while dynamically tuned quantization orders balance between accuracy and performance. On a single NVIDIA B200 GPU, our method computes low-rank approximations of large rectangular complex double-precision matrices over 14 × faster than cuSOLVER’s bi-diagonalization SVD and MAGMA’s column-pivoted QR, while preserving LAPACK-quality pivoting. Against cuSOLVER’s randomized SVD, hyacin remains consistently superior: 4.5 × and 1.8 × faster correspondingly with and without robustness options such as power iterations. Comparable speedups are observed across other GPU architectures, including the H100, A100, and RTX 4090. These results demonstrate that hyacin delivers LAPACK-level accuracy at unprecedented speed, outperforming both deterministic and randomized state-of-the-art GPU factorizations while maintaining strong numerical stability.