Hierarchical Sparse Cholesky

• Hierarchical sparse Cholesky methods are algorithms that efficiently factorize large symmetric positive definite matrices using hierarchical partitioning, nested dissection, and conditional-independence.
• They employ first-order and second-order low-rank compression techniques to control fill-in and reduce computational iterations, enhancing scalability and accuracy.
• Applications span Gaussian process regression, PDE solvers, and Bayesian estimation, achieving significant memory savings and near-linear computational complexity in large-scale problems.

Hierarchical sparse Cholesky refers to a set of algorithms and data structures designed to efficiently compute Cholesky factorizations for large, structured, and often kernel- or PDE-derived symmetric positive definite (SPD) matrices. These methods maintain scalability and accuracy by exploiting intrinsic low-rank or conditional-independence structure through hierarchical partitioning, blockwise low-rank compression, sparsity, and recursive or multilevel organization. Hierarchical sparse Cholesky methods appear prominently in numerical linear algebra, Gaussian process regression, inverse problems, spatio-temporal filtering, and in Bayesian estimation of sparse precision matrices, where classical dense or unstructured sparse algorithms incur prohibitive computational or memory costs.

1. Hierarchical Partitioning, Nested Dissection, and Block Structure

A core component is hierarchical partitioning of the problem size, often via graph-based nested dissection. The adjacency graph of the SPD matrix $A \in \mathbb{R}^{n \times n}$ is recursively partitioned by small "separator" sets $S^{(\ell)}$ that break the graph into nearly equal-sized interiors $I_j^{(\ell)}$. This process is repeated through $\ell$ levels until base cases are reached. The matrix $A$ is correspondingly reordered into a block structure that isolates the separators and organizes interior-interior, separator-interior, and separator-separator couplings. This enables a block $2\times 2$ decomposition at each stage: $$A = \begin{pmatrix}A_{II} & A_{IS} \ A_{SI} & A_{SS}\end{pmatrix}$$ where $I$ is an interior, $S$ a separator, and subsequent block eliminations induce Schur complements whose fill-in is controlled by the hierarchy (Klockiewicz et al., 2020).

2. Low-Rank Compression and Second-Order Sparsification

Off-diagonal blocks appearing in the Schur complements, especially in matrices arising from elliptic differential operators or kernel methods, often have rapidly decaying singular values and are numerically low-rank. Low-rank approximations

$A_{SI} \approx U \Sigma V^T$

are computed, exploiting decay in spectral content.

Hierarchical sparse Cholesky methods deploy two levels of approximation:

• First-order (linear) schemes compress off-diagonal blocks and drop small fill-in terms, resulting in an error $S^{(\ell)}$0 if the low-rank approximation error is $S^{(\ell)}$1.
• Second-order (quadratic) schemes (such as second-order accurate hierarchical approximate factorizations) explicitly retain quadratic terms arising from residuals of the low-rank approximation, ensuring that the norm of the global factorization error scales as $S^{(\ell)}$2. This ensures, for a fixed low-rank approximation error, a squared improvement in downstream error, especially in preconditioning contexts. The second-order approach can halve the number of preconditioned conjugate gradient (PCG) iterations required for convergence relative to first-order schemes, with negligible increase in factorization complexity (Klockiewicz et al., 2020).

Supernodal and rank-structured Cholesky variants also apply randomized or pivoted low-rank approximation to off-diagonal and (optionally) diagonal blocks within the hierarchical partitioning to maximize computational and memory efficiency (Chadwick et al., 2015).

3. Hierarchical Vecchia and Screening-Induced Cholesky Sparsity

For Gaussian process (GP) kernel matrices and spatial statistical models, the hierarchical Vecchia ("HV") approximation generalizes conditional-independence assumptions to construct sparse Cholesky factors. In these methods, each variable conditions only on a small subset of ancestors defined hierarchically (knots and groups), inducing both theoretical and practical sparsity in Cholesky factors of the covariance and corresponding precision matrices. Ordering strategies such as maximin or coarse-to-fine orderings (whereby Dirac-type measurements come first, followed by derivatives) facilitate exponential decay of Cholesky factor entries away from the block diagonal, leading to provable and rapidly decaying fill-in (Chen et al., 2023, Jurek et al., 2020, Jurek et al., 2022).

A key result is that, for suitable orderings and locality-based sparsity masks $S^{(\ell)}$3, the incomplete Cholesky factorization with this sparsity achieves near-linear complexity in problem size $S^{(\ell)}$4, and the Kullback-Leibler divergence and Frobenius norm error between the exact and approximate inverses can be bounded at $S^{(\ell)}$5 for sparsity radius parameter $S^{(\ell)}$6 (Chen et al., 2023). Conditioning on derivative measurements and PDE constraints can be combined seamlessly within this framework, providing scalability to complex and nonlinear PDE systems.

4. Algorithmic Structure and Practical Implementation

A representative pseudocode for hierarchical sparse Cholesky (e.g., for the Sparsified Nested Dissection "spaND" or kernel-based variants) comprises the following major steps (Klockiewicz et al., 2020, Chen et al., 2023, Jurek et al., 2020):

1. Hierarchical partition or ordering of variables (nested dissection, coarse-to-fine, maximin, supernodal aggregation).
2. At each level or supernode: Exact elimination or low-rank compression of interior/interior-separator blocks; update Schur complement using low-rank or quadratic schemes.
3. For GP or kernel matrices, explicit construction of sparsity mask $S^{(\ell)}$7 based on spatial proximity or conditional sets.
4. Incomplete Cholesky factorization with prescribed sparsity mask; block- or supernode-wise assembly, exploiting reuse via aggregation.
5. Recursive elimination or update of remaining separator systems.
6. Parameter selection (accuracy $S^{(\ell)}$8, rank thresholds) for optimizing the factorization/solve tradeoff.

Blockwise, supernodal, and compressed storage data structures are crucial for managing metadata, low-rank factors, and iterative updates (Chadwick et al., 2015).

5. Complexity, Storage, and Theoretical Guarantees

Complexity depends on problem geometry and off-diagonal numerical ranks $S^{(\ell)}$9:

• 2D grid-like problems: $I_j^{(\ell)}$0 factorization, $I_j^{(\ell)}$1 application, $I_j^{(\ell)}$2 memory.
• 3D: $I_j^{(\ell)}$3 factorization, $I_j^{(\ell)}$4 application, $I_j^{(\ell)}$5 memory (Klockiewicz et al., 2020).
• Hierarchical Vecchia (GP/statistics): $I_j^{(\ell)}$6 space, $I_j^{(\ell)}$7 time for $I_j^{(\ell)}$8-dimensional grids (Chen et al., 2023).
• Supernodal low-rank Cholesky: $I_j^{(\ell)}$9 idealized overall cost, with memory $\ell$0, where $\ell$1 is a threshold on diagonal block sizes (Chadwick et al., 2015).

For scalable spatio-temporal filtering and smoothing, the forward-filter/backward-sampler with HV approximation achieves $\ell$2 time and $\ell$3 memory, where $\ell$4 is the maximum neighbor set size and $\ell$5 the number of time points (Jurek et al., 2022).

Theoretical error guarantees include provably exponential decay of KL divergence and operator norm error in the conditioning-set or sparsity radius size (Chen et al., 2023, Jurek et al., 2020, Jurek et al., 2022).

6. Advanced Variants: Dynamic Sparsity and Hierarchical Bayesian Estimation

Adaptive reuse and dynamic updating of hierarchical partitionings in the presence of changing sparsity patterns are addressed via hierarchical graph decomposition and local reordering strategies. The Parth framework uses binary tree decompositions and change-localizing "dirty" block tracking, delivering order-of-magnitude speedups in symbolic reordering for Cholesky solves in applications with temporally coherent, spatially localized changes. Theoretical complexity drops from $\ell$6 per reordering to $\ell$7 with $\ell$8 for localized changes (Zarebavani et al., 2024).

In Bayesian estimation of sparse precision matrices, explicit hierarchical sparse Cholesky parametrizations (e.g., S-Bartlett) assign distributions over Cholesky factors with known or latent sparsity structure, supporting arbitrary, non-decomposable patterns. The S-Bartlett decomposes the prior as a product of independent (or conditionally independent) Gamma and Normal densities, bypassing intractable constants inherent in G-Wishart priors and permitting MCMC inference through dual-averaging HMC (Mastrantonio et al., 11 Jun 2025).

7. Applications, Performance, and Extensions

Hierarchical sparse Cholesky methods are broadly applied in:

• Direct solvers and preconditioners for large sparse PDE and kernel systems (Klockiewicz et al., 2020, Chadwick et al., 2015).
• Efficient and accurate solution of nonlinear and spatio-temporal PDEs via Gaussian processes (Chen et al., 2023, Jurek et al., 2022, Jurek et al., 2020).
• Large-scale Bayesian smoothing and filtering, outperforming low-rank and standard Laplace schemes (Jurek et al., 2022, Jurek et al., 2020).
• Dynamic contact, remeshing, and other computational graphics scenarios with evolving sparsity (Zarebavani et al., 2024).
• Hierarchical Bayesian estimation with flexible, interpretable priors on the space of SPD matrices (Mastrantonio et al., 11 Jun 2025).

Numerical studies report memory reductions by factors of 3–10 and wall-clock time reductions by 2–10 over conventional sparse direct and multilevel methods, with convergence behavior and numerical stability preserved or improved (Chadwick et al., 2015, Klockiewicz et al., 2020). Hierarchical Vecchia-based solvers achieve near-linear complexity on high-dimensional spatio-temporal data, and the second-order accurate approaches halve the iteration count in iterative solvers for a fixed preconditioner accuracy (Klockiewicz et al., 2020, Chen et al., 2023).

The field continues to evolve, with extensions to Laplace approximations for non-Gaussian data, nonlinear state evolution (via extended filtering), adaptive hierarchical decomposition for dynamic sparsity, and statistical frameworks for hierarchical prior modeling.

• Klockiewicz et al., "Second Order Accurate Hierarchical Approximate Factorization of Sparse SPD Matrices" (Klockiewicz et al., 2020)
• Schaefer et al., "Sparse Cholesky Factorization for Solving Nonlinear PDEs via Gaussian Processes" (Chen et al., 2023)
• Katsfuss & Jurek, "Scalable Spatio-Temporal Smoothing via Hierarchical Sparse Cholesky Decomposition" (Jurek et al., 2022)
• Chadwick & Bindel, "An Efficient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization" (Chadwick et al., 2015)
• Song et al., "Adaptive Algebraic Reuse of Reordering in Cholesky Factorization with Dynamic Sparsity Pattern" (Zarebavani et al., 2024)
• Ghosh et al., "A new hierarchical distribution on arbitrary sparse precision matrices" (Mastrantonio et al., 11 Jun 2025)
• Jurek & Katzfuss, "Hierarchical sparse Cholesky decomposition with applications to high-dimensional spatio-temporal filtering" (Jurek et al., 2020)