Performant Unified GPU Kernels for Portable Singular Value Computation Across Hardware and Precision

Abstract: This paper presents a portable, GPU-accelerated implementation of a QR-based singular value computation algorithm in Julia. The singular value ecomposition (SVD) is a fundamental numerical tool in scientific computing and machine learning, providing optimal low-rank matrix approximations. Its importance has increased even more in large-scale machine learning pipelines, including LLMs, where it enables low-rank adaptation (LoRA). The implemented algorithm is based on the classic two-stage QR reduction, consisting of successive matrix reduction to band form and bidiagonal form. Our implementation leverages Julia's multiple dispatch and metaprogramming capabilities, integrating with the GPUArrays and KernelAbstractions frameworks to provide a unified type and hardware-agnostic function. It supports diverse GPU architectures and data types, and is, to our knowledge, the first GPU-accelerated singular value implementation to support Apple Metal GPUs and half precision. Performance results on multiple GPU backends and data types demonstrate that portability does not require sacrificing performance: the unified function outperforms most linear algebra libraries (MAGMA, SLATE, rocSOLVER, oneMKL) for matrix sizes larger than 1024x1024, and achieves 80%-90% of the performance of cuSOLVER for large matrices.

• The paper presents a unified GPU-accelerated QR-based SVD method that offers portability across NVIDIA, AMD, Intel, and Apple platforms.
• It employs a two-stage algorithm with optimized panel factorization and fused trailing submatrix updates to maximize parallelism and minimize memory traffic.
• Performance benchmarks reveal the implementation outperforms vendor-specific libraries on large matrices while supporting low-precision FP16 for modern AI workloads.

Unified, Portable GPU Kernels for Singular Value Computation Across Hardware and Precision

Introduction

This paper presents a unified, high-performance GPU implementation of QR-based singular value computation, targeting portability across major hardware vendors (NVIDIA, AMD, Intel, Apple) and supporting multiple data precisions, including FP16. The implementation leverages Julia’s multiple dispatch and metaprogramming, integrating with GPUArrays and KernelAbstractions to provide a single, type- and hardware-agnostic API. The work addresses the limitations of existing SVD libraries, which are typically vendor-specific, lack support for Apple Metal, and do not provide unified support for low-precision arithmetic—a critical need for modern AI workloads such as LoRA in LLMs.

Algorithmic Approach and Implementation

The core algorithm is a two-stage QR-based SVD, following the reduction-to-band and reduction-to-bidiagonal paradigm. The first stage reduces the input matrix to band form using block Householder transformations, while the second stage further reduces the band matrix to bidiagonal and then diagonal form. The implementation focuses on maximizing parallelism and minimizing memory traffic, with custom GPU kernels for both panel factorization and trailing submatrix updates.

Figure 1: Phase 1 algorithm for reduction to band form, illustrating the sequence of RQ and LQ sweeps applied to diagonal tiles and their neighbors.

GPU Kernel Design

• Panel Factorization (GEQRT/TSQRT): Implemented as a single thread block per tile, maximizing L1 cache utilization and minimizing inter-thread communication. The kernel is sensitive to register pressure, with TILESIZE as a tunable parameter.
• Trailing Submatrix Update (UNMQR/TSMQR): Massively parallel, with each thread block processing a group of columns (COLPERBLOCK). The kernel is designed to maximize occupancy and memory coalescing.
• Fused Kernels: TSQRT and TSMQR operations are fused across tile rows, reducing kernel launch overhead and global memory traffic.

  Figure 2: Fused TSQRT and TSMQR kernels process all rows in a single launch, reducing memory traffic and synchronization overhead compared to traditional row-by-row processing.

Julia Abstractions

The implementation exploits Julia’s type inference and multiple dispatch to generate specialized LLVM code for each hardware and precision at compile time. GPUArrays.jl and KernelAbstractions.jl provide backend-agnostic interfaces, enabling the same high-level code to target CUDA, ROCm, oneAPI, and Metal backends.

Hyperparameter Tuning

Performance portability is achieved by abstracting hardware-sensitive parameters (TILESIZE, COLPERBLOCK, SPLITK) and tuning them per architecture and precision. Empirical results show up to 50% performance variation with different parameter choices, underscoring the necessity of autotuning for optimal results.

Performance Evaluation

Comparison with State-of-the-Art Libraries

The unified implementation is benchmarked against MAGMA, SLATE, cuSOLVER, rocSOLVER, and oneMKL across a range of hardware and matrix sizes.

Figure 3: Runtime ratio of the unified API to SLATE and MAGMA; the unified function outperforms SLATE for all data sizes and MAGMA for $n > 1024$.

Figure 4: Runtime ratio of the unified function to vendor libraries; the unified function outperforms rocSOLVER (MI250), cuSOLVER (RTX4060), and oneMKL (Intel PVC) for large matrices, and approaches cuSOLVER performance (80–90%) on A100/H100 for large $n$.

Key findings:

• Outperforms SLATE and MAGMA for all but the smallest matrices ($n > 1024$).
• Exceeds rocSOLVER and oneMKL for all tested sizes on AMD and Intel hardware.
• Achieves 80–90% of cuSOLVER performance on high-end NVIDIA GPUs for large matrices; outperforms cuSOLVER on consumer GPUs.
• First GPU-accelerated SVD for Apple Metal and first to support FP16 on GPU.

Portability Across Hardware and Precision

Figure 5: Runtime of the unified function for singular values across NVIDIA H100, AMD MI250, Apple M1, and Intel PVC for FP16 and FP32. FP16 on NVIDIA matches FP32 due to upcasting; Apple Metal and Intel show similar scaling.

The unified API demonstrates consistent scaling and performance across all tested hardware and precisions. On NVIDIA, FP16 and FP32 runtimes are identical due to lack of scalar FP16 support on CUDA cores, but FP16 enables larger problem sizes. Apple Metal and Intel results confirm the generality of the approach.

Kernel Runtime Breakdown

Figure 6: Relative runtime of panel factorization, trailing submatrix update, reduction to bidiagonal, and reduction to diagonal. As $n$ increases, the first stage (reduction to band) dominates, with the trailing update becoming the primary bottleneck.

The analysis reveals that for large matrices, the trailing submatrix update becomes the dominant cost, especially on GPUs with fewer SMs. This highlights the importance of optimizing BLAS3-like operations for scalability.

Numerical Accuracy

The implementation maintains high accuracy across all tested precisions and matrix sizes, with relative errors matching or closely tracking those of cuSOLVER. No precision-specific rescaling is applied, but future work may incorporate mixed-precision and rescaling strategies for improved robustness.

Implications and Future Directions

Practical Implications

• Unified API: Enables seamless deployment of SVD in heterogeneous environments, reducing code duplication and maintenance.
• FP16 Support: Facilitates low-memory, high-throughput SVD for AI workloads, including LoRA in LLMs.
• Apple Metal Support: Fills a critical gap for scientific computing on Apple hardware.
• Open Source: The implementation is available for integration into larger frameworks and as a drop-in replacement for vendor-specific kernels.

Theoretical Implications

• Performance Portability: Demonstrates that high performance can be achieved with a unified algorithmic implementation, provided that hardware-sensitive parameters are tunable.
• Abstraction Efficacy: Validates the use of high-level language abstractions (Julia) for portable, high-performance GPU computing.

Future Work

• Full SVD (Singular Vectors): Extension to compute singular vectors, enabling complete SVD functionality.
• Automated Tuning: Integration of autotuning for kernel parameters based on matrix size and hardware.
• Out-of-Core and Multi-GPU: Support for larger-than-memory problems and distributed execution.
• Non-Square and Tall-Skinny Matrices: Specialized algorithms for broader applicability.
• Integration with Task-Based Runtimes: Composability with distributed and pipelined scientific workflows.

Conclusion

This work establishes that a unified, portable GPU implementation of QR-based singular value computation is feasible and performant across major hardware vendors and data precisions. The approach leverages Julia’s abstraction capabilities to deliver a single API that matches or exceeds the performance of state-of-the-art, vendor-specific libraries for large matrices, while providing new capabilities such as FP16 and Apple Metal support. The results underscore the importance of hyperparameter tuning for performance portability and set the stage for further advances in unified, high-performance linear algebra for heterogeneous and AI-centric computing environments.