A GPU-resident Memory-Aware Algorithm for Accelerating Bidiagonalization of Banded Matrices

TL;DR

This work tackles the memory-bound bottleneck of the banded-to-bidiagonal reduction stage in two-stage SVD pipelines by introducing the first fully GPU-resident bulge-chasing algorithm, implemented in a portable Julia framework. It combines a bandwidth-tiling strategy with a memory-aware kernel design and a hardware-aware performance model to guide hyperparameter tuning across NVIDIA, AMD, Intel, and Apple GPUs and FP16/FP32/FP64 precisions. The approach delivers order-of-magnitude speedups over CPU libraries (PLASMA, SLATE) and scales with matrix bandwidth and size, enabling large-scale, fully GPU-based SVD pipelines. This work demonstrates that, with principled algorithm design and performance portability, memory-bound linear algebra can be accelerated on contemporary GPUs, reshaping the design of future GPU-accelerated SVD workflows.

Abstract

The reduction of a banded matrix to a bidiagonal form is a crucial step in the Singular Value Decomposition (SVD), a cornerstone of scientific computing and AI. Despite being a highly parallel algorithm, it was previously believed to be unsuitable for GPU computation because it is memory bandwidth-bound. Recent developments in GPU hardware, including larger L1 memory per Streaming Multiprocessor/Compute Unit, have changed that. We present the first GPU algorithm for reducing a banded matrix to bidiagonal form as part of the NextLA$.$jl open-source software package. Our algorithm is based on previous CPU-based multicore parallel cache-efficient bulge chasing algorithms and adapted to optimize for GPU throughput. We leverage Julia Language's Array abstractions and KernelAbstractions to implement a single hardware- and data precision-agnostic function on NVIDIA, AMD, Intel, and Apple Metal GPUs for half, single, and double precision, and examine performance optimization across hardware architectures and data precision. We also develop a hardware-aware performance model and identify key hyperparameters, such as inner tilewidth and block concurrency, that govern optimal GPU execution for bandwidth-bound workloads. We demonstrate highly parallel bandwidth-bound algorithm on the GPU can outperform CPU-based implementations: the GPU algorithm outperforms multithreaded CPU High-Performance libraries PLASMA and SLATE as of matrix size 1024 x 1024 and by a factor over 100 for matrices of 32k x 32k. In addition, the performance of the algorithm increases linearly with matrix bandwidth size, making faster reduction of larger matrix bandwidths now also possible. With this work, we break memory bandwidth barriers, as well as matrix bandwidth barriers, resulting in orders-of-magnitude faster algorithms for the reduction of banded matrices to bidiagonal form on the GPU.

Figures (7)

• Figure 1: Illustration of the process of successive reduction of the bandwidth of the matrix by inner tilewidth for inner tilewidth two and total bandwidth 6. The proposed strategy increases memory locality and cache reuse, crucial for the memory-bound algorithm on the GPU, where limited memory and memory bandwidth are available.
• Figure 2: Visualization of the bandwidth-tiled, sweep-based GPU implementation of band-to-bidiagonal reduction using Householder reflectors. Each row sweep (e.g., Row 1, Row 2) processes a bandwidth tile, chasing bulges (purple dots) down the matrix by applying orthogonal transformations to annihilate matrix elements (red dots). Excess bulges generated during one row’s sweep are deferred (gray squares) and picked up in subsequent sweeps, enabling pipelined execution. The right panel illustrates how each thread block cooperatively computes transformations across rows using chunked tiling, where each thread handles a row and loops over multiple elements. This design enables fine-grained parallelism, improved memory locality, and GPU-resident execution across bandwidth tiles.
• Figure 3: Relative error of singular values computed via GPU-based reduction to bidiagonal form, followed by LAPACK BDSDC in double precision. The boxplots show errors over 30 trials per precision (FP64/FP32/FP16) across matrix sizes and bandwidths. Results confirm stable accuracy with increasing bandwidth and only moderate error growth with size.
• Figure 4: Performance of the GPU-based band reduction algorithm across kernel hyperparameters, hardware, and precisions. Parallel coordinate plots show color-coded relative runtimes for combinations of three hyperparameter max blocks, tile width, and threads per block on NVIDIA H100 (FP32, FP64) and AMD MI300X (FP32), at bandwidths 32 and 128. Yellow lines indicate faster configurations. While tile width—optimally matching cache line sizes (32 for FP32, 16 for FP64)—has the greatest impact, both max blocks and threads per block also significantly influence performance, especially at wider bandwidths.
• Figure 5: Performance gains from GPU architectural advancements, showing the relative performance loss of using older GPU architecture for varying matrix sizes and bandwidths across AMD and NVIDIA, and highlighting the impact of larger L1 and L2 caches.