Numerically Stable Cholesky-QR on GPU via Mixed-Precision Randomized Preconditioning

Abstract

Cholesky-QR is among the fastest algorithms for computing the thin QR factorization of tall-and-skinny matrices on GPUs, relying entirely on BLAS-3 operations. However, it is numerically unstable: forming the Gram matrix squares the condition number, causing breakdown when $\hat{I}{º}_2(\mathbf{A})\gtrsim 10^8$. We present MRCQR (Mixed-Precision Randomized Cholesky-QR), a stable GPU algorithm that addresses this limitation. MRCQR uses a subsampled randomized trigonometric transform to construct a preconditioner ${\mathbf{R}}_s$ that reduces $\hat{I}{º}_2(\mathbf{A}{\mathbf{R}}_s^{-1})$ to near unity with high probability, then applies Cholesky-QR in double precision to the preconditioned matrix. The key insight -- supported by perturbation analysis -- is that the preconditioner requires far less accuracy than the final result: single (FP32) precision suffices when $\hat{I}{º}_2(\mathbf{A})\lesssim 10^8$, and half (FP16) when $\hat{I}{º}_2(\mathbf{A})\lesssim 10^4$. MRCQR produces an explicit orthogonal factor $\hat{\mathbf{Q}}$ satisfying $\Vert \mathbf{I}-{\hat{\mathbf{Q}}}^{\top }\hat{\mathbf{Q}}{\Vert }_2=\mathcal{O}(\mathbf{u})$ ($\mathbf{u}\approx 10^{-16}$, double-precision unit roundoff) for condition numbers up to $10^{16}$, far beyond the $10^8$ limit of CholQR2. Experiments on an NVIDIA H100 GPU show that MRCQR (FP16) outperforms rand-cholQR by $1.4$--$1.8\times$ across all tested column counts and is $1.8$--$13.5\times$ faster than cuSOLVER geqrf, while the FP16 sketch (used when $\hat{I}{º}_2(\mathbf{A})\lesssim 10^4$) is $2\times$ cheaper than FP64 at no accuracy cost.