Randomized Cholesky (RCHOL)

• Randomized Cholesky (RCHOL) is a class of algorithms that approximates symmetric positive definite matrices via hierarchical partitioning and randomized low-rank compression.
• It employs randomized range finding, power iterations, and QR factorization to compress off-diagonal blocks while maintaining numerical stability and reducing computational cost.
• RCHOL enhances scalability in applications such as PDE solvers, spatio-temporal filtering, and state-space inference by drastically lowering memory usage and speeding up computations.

Randomized Cholesky (RCHOL) refers to a class of algorithms for the factorization and approximation of large, sparse, or structured symmetric positive definite (SPD) matrices using randomized methods within a hierarchical or block-structured Cholesky framework. These algorithms are designed to improve computational scalability and memory efficiency by leveraging low-rank structure in the matrix blocks, stochastic compression, and hierarchical partitioning. RCHOL methods have seen broad impact in PDE solvers, spatial statistics, and large-scale Kalman or state-space inference, where exact dense Cholesky decompositions are infeasible.

1. Mathematical Foundations and Algorithmic Structure

The essential goal of randomized or rank-structured Cholesky algorithms is to decompose an SPD matrix $A\in\mathbb{R}^{n\times n}$ such that $A \approx LL^\top$, where $L$ is lower-triangular, with structure exploited for algorithmic and storage efficiency. In hierarchical RCHOL, $L$ is assembled blockwise, along a supernodal or nested dissection partitioning of the matrix.

For each supernode or separator block, the off-diagonal coupling blocks (which connect to the rest of the matrix) are approximated in low rank:

$L_j^O \approx V_j U_j^\top$

where $U_j$ is orthonormal, $V_j$ is computed via efficient (often randomized) range finding, and the approximation maintains symmetry and positive-definiteness of updates to the Schur complement. The randomized aspect enters in the use of random test matrices (for example, $G \sim \mathcal{N}(0, I)$) to probe the range and estimate the dominant low-rank components of the off-diagonal block efficiently, followed by projection (e.g., via power iteration and QR factorization) (Chadwick et al., 2015).

Let $C_j$ be the columns of supernode $j$, and $A \approx LL^\top$0 the set of rows below $A \approx LL^\top$1 where fill may occur. The compressed representation only stores the factors $A \approx LL^\top$2, reducing memory and arithmetic costs compared to explicit storage of large dense blocks. Diagonal blocks within large supernodes may themselves be recursively partitioned and compressed using hierarchical matrix techniques.

2. Hierarchical Matrix Representation and Compression

RCHOL algorithms rely on a hierarchical block partitioning, usually induced by a fill-reducing ordering such as nested dissection. The matrix is split recursively into interiors and interfaces at each level; the hierarchy induces a block structure where separators (interfaces) are compressed using low-rank approximations (Chadwick et al., 2015, Klockiewicz et al., 2020). The process is as follows:

1. At the $A \approx LL^\top$3-th level, the variables are divided into interiors $A \approx LL^\top$4 and interface (boundary) blocks $A \approx LL^\top$5.
2. The Schur complement on $A \approx LL^\top$6 is computed after eliminating $A \approx LL^\top$7 by block Cholesky.
3. Off-diagonal blocks between separator blocks and their neighbors are approximated in low rank using randomized projections, with the worst-case spectral error controlled by the chosen rank and truncation threshold $A \approx LL^\top$8.
4. The compression of off-diagonal blocks is performed via randomized sampling of the range, QR factorization, and subsequent projection, in order to represent large separators with minimal memory overhead.

When combined with a further hierarchical subdivision within each diagonal block (i.e., “H-matrix” structure), this yields recursive blockwise compression suitable for problems with multiscale couplings, such as those from elliptic PDEs or spatial Gaussian processes (Chadwick et al., 2015).

3. Randomized Compression: Algorithms and Guarantees

Randomized Cholesky compression exploits the rapid decay in singular values of off-diagonal Schur complements to construct compact low-rank representations. The algorithm proceeds as follows (Chadwick et al., 2015, Klockiewicz et al., 2020):

• Generate a random Gaussian matrix $A \approx LL^\top$9.
• Multiply: $L$0 (performed via left-looking traversal over previously factored supernodes).
• Apply a fixed number $L$1 of power iterations to improve spectral decay, alternating multiplication by $L$2 and $L$3.
• Orthonormalize the samples to obtain $L$4.
• Project to compute $L$5.
• Store $L$6 without explicit construction of the full off-diagonal block.

Error analysis shows that the norm of the approximation error depends on the tail singular values of the compressed block and can be controlled to any prescribed tolerance $L$7 by adjusting the sketch dimension and power iterations (Klockiewicz et al., 2020). In the full “second-order” scheme, quadratic error terms in the Schur complement involving the discarded fine components are retained, leading to a factorization error of $L$8 in operator norm, as opposed to $L$9 for simple truncation. This higher order accuracy leads to approximately halved Preconditioned Conjugate Gradient (PCG) iteration counts and sharper spectral preconditioners without significant additional cost (Klockiewicz et al., 2020).

4. Computational Complexity and Storage

The hierarchical and randomized framework achieves near-optimal complexity for large systems:

Algorithm	Factorization Time	Memory	Comments
Direct Cholesky	$L$0	$L$1	Prohibitive for large $L$2
Low-rank FFBS	$L$3	$L$4	Can miss fine scale dependence
Randomized/Hierarchical Cholesky	$L$5 (Chadwick et al., 2015), $L$6 (Jurek et al., 2022)	$L$7 (Chadwick et al., 2015), $L$8 (Jurek et al., 2022)	Richer factors; more accurate

The randomized low-rank updates reduce both time and memory associated with off-diagonal blocks, while batch-wise or streaming construction of $L$9 ensures compatibility with distributed and parallel computing architectures. Reporting from (Chadwick et al., 2015) demonstrates memory reductions of over an order of magnitude versus direct Cholesky (e.g., $L_j^O \approx V_j U_j^\top$0 GB vs $L_j^O \approx V_j U_j^\top$1 GB for $L_j^O \approx V_j U_j^\top$2), while also achieving superior solve times compared to Jacobi, ICC, and AMG preconditioning for large-scale PDE discretizations.

5. Applications in Spatio-Temporal Statistics and Filtering

The RCHOL paradigm is crucial in scalable inference for high-dimensional Gaussian and non-Gaussian state-space models, prominently in large-scale spatio-temporal smoothing and filtering. Through the hierarchical Vecchia approximation, the sparsity structure derived from conditional independence directly determines the sparsity mask for the Cholesky factor (Jurek et al., 2022, Jurek et al., 2020). This approach underpins efficient implementations of the forward-filter-backward-sampler (FFBS) for high-dimensional time series and spatial models.

Key steps include:

• Construction of a multi-level Vecchia tree/graph, with a max–min variable ordering.
• Computation of a sparsity-masked Cholesky factor respecting the hierarchical conditional sets.
• All filtering, smoothing, and simulation steps performed with only $L_j^O \approx V_j U_j^\top$3 storage and $L_j^O \approx V_j U_j^\top$4 total computation.

Empirical results show that modest neighborhood sizes (e.g., $L_j^O \approx V_j U_j^\top$5–$L_j^O \approx V_j U_j^\top$6) retain predictive accuracy within a few percent of the full dense FFBS (as measured by CRPS), with typical speedups of $L_j^O \approx V_j U_j^\top$7 and superior accuracy to low-rank approximations at comparable computational cost (Jurek et al., 2022).

6. Advanced Error Bounds and Theoretical Guarantees

Theoretical properties of RCHOL with hierarchical or randomized block compression schemes include:

• Given any chosen sparsity mask $L_j^O \approx V_j U_j^\top$8, the hierarchical sparse Cholesky factorization yields the optimal KL approximation to $L_j^O \approx V_j U_j^\top$9 within the class of distributions induced by sparse Cholesky factors with pattern $U_j$0 (Jurek et al., 2022).
• For hierarchical Vecchia patterns, error in the log-determinant or inverse-trace decays algebraically or exponentially as the neighborhood size $U_j$1 increases, under standard regularity (e.g., Matérn covariance smoothness) (Jurek et al., 2022, Jurek et al., 2020).
• In randomized low-rank block compression, the norm of the error introduced at each interface block is $U_j$2 (second-order scheme) and the overall impact on PCG iteration counts is a halving compared to first-order truncations (Klockiewicz et al., 2020).

These features allow practitioners to explicitly trade off accuracy versus cost by selecting rank or neighborhood size and truncation parameters.

7. Practical Considerations and Implementation Issues

RCHOL algorithms require integration with fill-reducing permutation strategies (e.g., nested dissection), careful management of data structures for block and off-diagonal storage, and tuning of block sizes and rank thresholds for practical performance. In typical usage, interfaces with highly optimized BLAS/LAPACK numerical kernels are retained for diagonal steps, while random projection and low-rank updates are parallelizable and scalable to very large problem sizes (Chadwick et al., 2015).

The approach is robust to a wide array of SPD matrices, including those with challenging spectral properties and variable coefficients, offering a “black-box” preconditioner that outperforms standard ICC or AMG for large-scale scientific computing tasks. Numerical benchmarks confirm scalability to millions of unknowns where direct methods become intractable (Chadwick et al., 2015).

---

Key references: "An Efficient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization" (Chadwick et al., 2015), "Second Order Accurate Hierarchical Approximate Factorization of Sparse SPD Matrices" (Klockiewicz et al., 2020), "Scalable Spatio-Temporal Smoothing via Hierarchical Sparse Cholesky Decomposition" (Jurek et al., 2022), and "Hierarchical sparse Cholesky decomposition with applications to high-dimensional spatio-temporal filtering" (Jurek et al., 2020).