Source-linked AI summary
A GPU-Accelerated Blocked Adaptive Randomized Range Finder Based on an Implicit Householder QR Decomposition
Carolin Penke, Andreas Herten
TL;DR
Low-rank applications often need only an adaptive orthonormal basis rather than an expensive full SVD. This paper introduces a blocked adaptive randomized range finder using implicit Householder QR and hybrid CPU–GPU execution, preserving numerical quality while reducing the largest tested runtime from 9.91 s on the CPU to 0.407 s on an NVIDIA GH200.
Problem
Computing a dominant low-dimensional subspace basis can require substantially less work than a full SVD, while efficient basis construction remains important for reduced-space methods such as GaLore.
Method
The method uses blocked Householder reflectors to represent the adaptive randomized basis implicitly, avoid explicit reorthogonalization, and overlap CPU panel factorization with GPU matrix updates.
Results
0.407 s versus 9.91 s on the CPU for the largest tested case, while the Householder variant preserves orthogonality and low residuals in extreme regimes where reorthogonalized BCGS deteriorates.
Takeaways & Limitations
The formulation provides a numerically robust and efficient GPU-oriented approach to adaptive low-rank matrix approximation.
Takeaways & Limitations
Further operation and data-movement analysis is needed to determine when proposed variants, including fused updates or compact WY representations, are beneficial.
Abstract
from arXiv · showhide
Low-rank methods can reduce the memory and computational requirements of deep neural network training in approaches such as GaLore. Randomized range finders offer an attractive alternative to singular value decompositions, particularly when the required rank is determined adaptively from a prescribed approximation tolerance. We introduce a blocked adaptive randomized range finder based on an implicit Householder QR decomposition and an optimized hybrid CPU--GPU implementation. The proposed method avoids explicit reorthogonalization. Numerical experiments show that it preserves orthogonality and approximation accuracy in regimes where block Gram--Schmidt without reorthogonalization becomes unstable. The blocked formulation exposes matrix--matrix operations and enables overlap of CPU panel factorization with GPU updates. On an NVIDIA GH200, the overlapped implementation reduces the runtime for the largest tested matrix from 9.91 seconds on the CPU to 0.407 seconds. The method provides a stable and efficient building block for low-rank approximation on heterogeneous systems with applications in computational science and engineering.
I. INTRODUCTION
The paper motivates adaptive low-rank subspace approximation for computational science, engineering, and machine learning, where full SVDs are expensive and only a dominant basis may be needed. It introduces a blocked Householder-based randomized range finder designed for numerical stability and efficient hybrid CPU–GPU execution.
- Dominant low-dimensional subspaces support applications in computational science, engineering, and machine learning.
- Low-rank training methods such as GaLore reduce memory consumption by projecting gradients and optimizer states into reduced subspaces.
- Full SVD computation can perform substantially more work than necessary when only a low-dimensional dominant subspace basis is required.
- Randomized range finders construct dominant-subspace bases either for a prescribed rank or adaptively until a specified approximation tolerance is reached.
- The proposed blocked Householder formulation avoids explicit reorthogonalization, exposes matrix–matrix operations, and overlaps CPU panel factorization with GPU updates.
II. BACKGROUND AND RELATED WORK
Prior adaptive randomized range finders use sampled blocks, residual-based stopping criteria, and Gram–Schmidt orthogonalization, while blocked Householder QR offers a stable alternative suited to level-3 operations and accelerator execution.
- Adaptive randomized range finders increase rank until a prescribed tolerance σ is met, rather than fixing the rank in advance.
- Residual norms provide stopping criteria, including an in-place residual norm and a Frobenius-norm criterion based on generated panels.
- Blocked adaptive formulations process sampled panels and commonly rely on Gram–Schmidt with explicit reorthogonalization to maintain basis orthogonality.
- Blocked Gram–Schmidt improves computational efficiency through BLAS level-3 operations but does not provide unconditional numerical stability.
- Householder QR avoids explicit reorthogonalization, offers strong computed-factor error bounds, and supports blocked level-3 operations through compact WY representations.
III. A RECURSIVE QB FACTORIZATION BASED ON HOUSEHOLDER QR
This section derives a recursive QB factorization that represents the randomized range basis implicitly with blocked Householder reflectors. The formulation avoids Gram–Schmidt orthogonalization while exposing matrix–matrix operations for efficient updates.
- The method implicitly performs a classic blocked QR factorization to construct a randomized QB decomposition.
- The adaptive formulation computes reflector representations successively because the target rank r is not known in advance.
- Applying the compact WY representation reduces operations to a series of matrix–matrix products.
- Householder vectors and scaling or triangular block factors are stored without explicitly forming the product of orthogonal reflectors.
- Unlike standard QR, the updated matrix width remains unchanged even as processed panel heights decrease.
- The recursive update uses the lower updated matrix block to determine when a low-rank approximation A ≈ Q_1B_1 is sufficiently accurate.
IV. THE HOUSEHOLDER BLOCKED ADAPTIVE RANDOMIZED RANGE FINDER
This section turns the recursive Householder construction into an adaptive blocked range finder. It samples Gaussian panels, builds Householder block columns and projected matrix rows, and stops when the residual criterion reaches the prescribed tolerance.
- The algorithm partitions A into blocks and stores Householder vectors in V while accumulating block rows of B in A’s memory locations.
- Algorithm 2 initializes E with the squared Frobenius norm of A, copies A into B, and repeatedly samples Gaussian matrices Ω.
- Each iteration factors a sampled panel with QR, producing a block of Householder vectors V and a triangular factor T.
- The procedure returns rank r, Householder vectors V, projected matrix B, and block factors T_0,...,T_{i−1} satisfying a Frobenius-norm error tolerance.
- The loop is guaranteed to terminate after at most min(j + 1, k + 1) steps, while early stopping can occur when E ≤ σ^2.
V. HYBRID CPU–GPU IMPLEMENTATION
The hybrid implementation assigns sequential panel factorizations to the CPU and matrix updates to the GPU. An overlapped variant independently updates the sampled panel, shortening the critical path at the cost of extra operations.
- The implementation aims to keep the GPU saturated with high-arithmetic-intensity BLAS-3 operations.
- The sequential implementation waits for the full update of B before sampling the next panel, whereas the overlapped variant updates the sampled panel independently.
- CPU panel factorization is overlapped with GPU matrix updates to expose parallelism in the hybrid implementation.
- The overlap strategy shortens the critical path but increases the operation count through redundant sampled-column updates.
- The method assumes input matrices reside on the GPU and uses MAGMA as the implementation basis, with GPU-suitable operations implemented primarily as matrix–matrix multiplications.
VI. NUMERICAL STABILITY AND PERFORMANCE
The Householder range finder preserves approximation quality and orthogonality more reliably than non-reorthogonalized BCGS, while its blocked CPU–GPU implementation accelerates large problems through overlap and integrated operations.
- Numerical stability: Without reorthogonalization, BCGS loses orthogonality immediately, and its residual stopping quantity remains misleadingly low as B absorbs basis defects.Reorthogonalization remedies these issues, but the Householder method avoids explicit reorthogonalization.
- Numerical stability: The Householder variant preserves orthogonality and low residuals near the single-precision noise floor, including when the prescribed singular values fall below machine precision.Reorthogonalized BCGS accuracy deteriorates near the end, whereas the Householder method does not show this problem.
- Numerical stability: The Frobenius-norm stopping quantity becomes inaccurate at k = 384 because cancellation makes it negative and forces it to be capped at zero.A spectral-norm stopping criterion for the lower part of B is identified as future work.
- Runtime performance: 0.407 s is the largest-case runtime for the two-queue implementation, compared with 9.91 s on the CPU and 0.728 s with one GPU queue.The two-queue variant becomes increasingly advantageous as matrix size grows because overlapping panel processing with matrix updates exposes additional parallelism.
- Runtime performance: 56× is the reported maximum speedup over RSVDpack for m = n = 16384.RSVDpack runtimes include host–device and device–host transfers, while the proposed implementation assumes input and output remain in GPU memory.
- Runtime performance: Even with the same number of basis vectors, the proposed method outperforms fixed-rank variants for large matrices by interleaving sampling, QR factorization, and B = Q^T A computation.The integrated blocked iteration improves GPU utilization relative to separate sequential stages.
- Runtime performance: For small target ranks r = 128 and 256, cuSOLVER fixed-rank range finding achieves the best performance.In this regime, Algorithm 2 is dominated by CPU panel factorization and construction of the T matrix.
VII. TOWARDS EFFICIENT LOW-RANK DEEP LEARNING
The paper connects adaptive randomized range finding to resource-efficient low-rank deep learning, replacing heuristic rank selection with tolerance-based determination while highlighting subspace-consistency issues across training steps.
- Randomized numerical linear algebra is presented as a promising tool for resource-efficient low-rank deep learning.
- GaLore represents optimizer states in a low-rank gradient subspace, but typically treats the rank r as a heuristically chosen hyperparameter.
- The adaptive range finder uses approximation tolerance σ as the hyperparameter, determining a suitable rank automatically.
- Adaptive methods can vary subspace dimensionality across consecutive training steps, making direct addition of optimizer states represented in different subspaces not meaningful.A linear transformation is applied to ensure subspace consistency.
VIII. CONCLUSION AND OUTLOOK
The proposed Householder-based blocked adaptive randomized range finder combines numerical robustness with efficient hybrid CPU–GPU execution. Its remaining challenges concern stopping-criterion accuracy, further optimization, distributed deep-learning integration, and comparison with reorthogonalized Gram–Schmidt methods.
- Conclusion: 0.407 s versus 9.91 s reduced the largest tested GH200 case relative to the CPU through overlapped CPU panel factorization and GPU matrix updates.The two-queue implementation demonstrates the runtime benefit of overlapping these operations.
- Outlook: Cancellation near the numerical noise floor can reduce the accuracy of the inexpensive Frobenius-norm stopping criterion.The paper proposes investigating more robust implementations, spectral-norm estimates, and power iterations for slowly decaying singular values.
- Outlook: Further optimization opportunities include fusing matrix updates with residual-norm computation and using alternative implementations or representations such as cuSOLVER, Kokkos, or compact WY.Determining when these variants are advantageous requires operation and data-movement analysis.
- Outlook: Applying the method to GaLore-like frameworks requires extensions for distributed matrix layouts arising from fully sharded data parallelism and model parallelism.
- Outlook: Reorthogonalized block Gram–Schmidt methods remain a promising alternative for adaptive randomized range finders, particularly in distributed settings.
- Conclusion: The method is presented as a numerically robust and efficient GPU-oriented approach to adaptive low-rank matrix approximation.