Dense Matrices Are Alike; Sparse Matrices Are Sparse in Their Own Way: A Structure-Adaptive Tile Cholesky Factorization
Abstract
Sparse direct Cholesky solvers fix one data structure for an entire matrix, but symmetric positive definite systems range from nearly dense to irregular, sometimes mixing both within one matrix. We let the data structure follow the sparsity structure, across matrices and across tiles within a matrix. Before numerical work starts, a lightweight selector captures the sparsity pattern of the Cholesky factor and routes the matrix, on one static shared-memory schedule, to one of three regimes: dense, sparse, or an intermediate semisparse regime. Dense tiles are stored in full; the semisparse regime uses a new active-column tile that keeps only the columns the factorization will touch, so banded or arrowhead-shaped structures common in spatial models still reach BLAS-3 efficiency. The approach suits the integrated nested Laplace approximation (INLA): one sparsity pattern factorized thousands of times with different values, so the one-time analysis cost is amortized and every factorization saving compounds. We evaluate the technique on 60 SPD matrices, comparing against MUMPS, PaStiX, CHOLMOD, symPACK, and Intel oneMKL PARDISO on Intel Xeon and AMD EPYC nodes; the selector alone achieves the lowest total factorization time in every regime. Summed over the suite, it beats the best fixed single-structure mode by 1.6 to 2.6x, and every alternative by 1.8 to 12.5x on Intel and 2.6 to 10.1x on AMD, with the largest gains on the most expensive factorizations. It trades more one-time analysis for less time per factorization, pulling ahead by the third factorization of a given pattern. As a first GPU extension, the dense route on one NVIDIA A100, with the factor resident on the device, runs 1.2 to 6.3x faster than on the faster CPU node, the margin widening with factor size. Solver, Python/R/Julia interfaces, benchmark suite, and results are open at https://github.com/esmail-abdulfattah/sTiles.
Related
- Adaptive multi-resolution Gaussian processes: Scalable exact inference with naturally data-sparse covariance matrices
- Right Precision, Right Place: Runtime Adaptations for Mixed-Precision Cholesky Factorization
- Fast nonparametric estimation of sparse precision matrices with a priori sparsity structure
- The ϕ-conjugation of quaternionic matrices and generalized Autonne-Takagi factorization
- An Alternating Minimization Algorithm for Federated Sparse Orthogonal Nonnegative Matrix Factorization
BibTeX
@misc{fattah2026dense,
title = {{Dense Matrices Are Alike; Sparse Matrices Are Sparse in Their Own Way: A Structure-Adaptive Tile Cholesky Factorization}},
author = {Esmail Abdul Fattah and Hatem Ltaief and Håvard Rue and David E. Keyes},
howpublished = {arXiv},
year = {2026},
month = sep,
doi = {10.48550/arxiv.2609.29765},
eprint = {2609.29765},
archivePrefix = {arXiv},
url = {https://arxiv.org/abs/2609.29765},
note = {CC BY-NC-ND 4.0},
}
RIS
TY - UNPB
TI - Dense Matrices Are Alike; Sparse Matrices Are Sparse in Their Own Way: A Structure-Adaptive Tile Cholesky Factorization
AU - Esmail Abdul Fattah
AU - Hatem Ltaief
AU - Håvard Rue
AU - David E. Keyes
PY - 2026
DA - 2026/09/24
PB - arXiv
DO - 10.48550/arxiv.2609.29765
UR - https://arxiv.org/abs/2609.29765
AB - Sparse direct Cholesky solvers fix one data structure for an entire matrix, but symmetric positive definite systems range from nearly dense to irregular, sometimes mixing both within one matrix. We let the data structure follow the sparsity structure, across matrices and across tiles within a matrix. Before numerical work starts, a lightweight selector captures the sparsity pattern of the Cholesky factor and routes the matrix, on one static shared-memory schedule, to one of three regimes: dense, sparse, or an intermediate semisparse regime. Dense tiles are stored in full; the semisparse regime uses a new active-column tile that keeps only the columns the factorization will touch, so banded or arrowhead-shaped structures common in spatial models still reach BLAS-3 efficiency. The approach suits the integrated nested Laplace approximation (INLA): one sparsity pattern factorized thousands of times with different values, so the one-time analysis cost is amortized and every factorization saving compounds. We evaluate the technique on 60 SPD matrices, comparing against MUMPS, PaStiX, CHOLMOD, symPACK, and Intel oneMKL PARDISO on Intel Xeon and AMD EPYC nodes; the selector alone achieves the lowest total factorization time in every regime. Summed over the suite, it beats the best fixed single-structure mode by 1.6 to 2.6x, and every alternative by 1.8 to 12.5x on Intel and 2.6 to 10.1x on AMD, with the largest gains on the most expensive factorizations. It trades more one-time analysis for less time per factorization, pulling ahead by the third factorization of a given pattern. As a first GPU extension, the dense route on one NVIDIA A100, with the factor resident on the device, runs 1.2 to 6.3x faster than on the faster CPU node, the margin widening with factor size. Solver, Python/R/Julia interfaces, benchmark suite, and results are open at https://github.com/esmail-abdulfattah/sTiles.
ER -