GPU-Accelerated Cholesky Factorization of Block Tridiagonal Matrices

TL;DR

This work addresses the bottleneck of solving alternating block-tridiagonal linear systems in real-time applications by introducing a GPU-accelerated Cholesky framework that leverages nested-dissection–style multi-stage permutations to achieve $O(\log_2(N)\,n^3)$ complexity. It combines partitioning and operation-level parallelism (via CUDA streams and atomic updates) with fused and tiled kernels to exploit dense block structure and maximize GPU throughput. Comprehensive GPU experiments show substantial speedups over sparse solvers (e.g., $>100\times$ vs QDLDL) and CPU implementations (up to $40\times$ with BLASFEO) and even outperform NVIDIA’s CUDSS by about $2\times$ for both factorization and solve. The results demonstrate the practical impact for long-horizon, real-time optimization in robotics, MPC, and Kalman filtering, and the approach is open-source for further adoption and extension.

Abstract

This paper presents a GPU-accelerated framework for solving block tridiagonal linear systems that arise naturally in numerous real-time applications across engineering and scientific computing. Through a multi-stage permutation strategy based on nested dissection, we reduce the computational complexity from $\mathcal{O}(Nn^3)$ for sequential Cholesky factorization to $\mathcal{O}(\log_2(N)n^3)$ when sufficient parallel resources are available, where $n$ is the block size and $N$ is the number of blocks. The algorithm is implemented using NVIDIA's Warp library and CUDA to exploit parallelism at multiple levels within the factorization algorithm. Our implementation achieves speedups exceeding 100x compared to the sparse solver QDLDL, 25x compared to a highly optimized CPU implementation using BLASFEO, and more than 2x compared to NVIDIA's CUDSS library. The logarithmic scaling with horizon length makes this approach particularly attractive for long-horizon problems in real-time applications. Comprehensive numerical experiments on NVIDIA GPUs demonstrate the practical effectiveness across different problem sizes and precisions. The framework provides a foundation for GPU-accelerated optimization solvers in robotics, autonomous systems, and other domains requiring repeated solution of structured linear systems. The implementation is open-source and available at https://github.com/PREDICT-EPFL/socu.

Figures (22)

• Figure 1: Submatrix-Cholesky
• Figure 2: Theoretical speed-up of parallel factorization of $\Psi_p$ compared to sequential factorization of $\Psi$, using $p$ number of threads. The red line boundary divides actual speed-ups from slowdowns.
• Figure 3: Elimination trees for different permutations of $A$ for $N=20$.
• Figure 4: Speed-up of parallel factorization of $\Phi_1$ compared to sequential factorization of $\Psi$. The red line boundary divides actual speed-ups from slowdowns.
• Figure 5: Visualization of the Cholesky factorization of $\Phi_\infty$ for $N=20$ using Algorithm \ref{['alg:factorization_multi_stage']}. Visualized is the parallel iteration for $s=2$, i.e., the second iteration. The full arrows operations updating data, while the dashed arrow shows a deferred operation.

Theorems & Definitions (2)

• Proposition 1
• proof