---
title: Randomized Cholesky (RCHOL)
url: https://www.emergentmind.com/topics/randomized-cholesky-rchol
type: topic
---

# Randomized Cholesky (RCHOL)

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) [1507.05593].

Let $C_j$ be the columns of supernode $j$, and $R_j$ the set of rows below $C_j$ where fill may occur. The compressed representation only stores the factors $(U_j, V_j)$, 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 [1507.05593][2007.00789]. The process is as follows:

1. At the $l$-th level, the variables are divided into interiors $I_l$ and interface (boundary) blocks $B_l$.
2. The Schur complement on $B_l$ is computed after eliminating $I_l$ 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 $\varepsilon$.
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 [1507.05593].

## 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 [1507.05593][2007.00789]:

- Generate a random Gaussian matrix $G \in \mathbb{R}^{|R_j|\times (r_j + p)}$.
- Multiply: $B = L_j^O G$ (performed via left-looking traversal over previously factored supernodes).
- Apply a fixed number $s$ of power iterations to improve spectral decay, alternating multiplication by $L_j^O$ and $(L_j^O)^\top$.
- Orthonormalize the samples to obtain $U_j$.
- Project to compute $V_j$.
- Store $L_j^O \approx V_j U_j^\top$ 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 $\varepsilon$ by adjusting the sketch dimension and power iterations [2007.00789]. 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 $O(\varepsilon^2)$ in operator norm, as opposed to $O(\varepsilon)$ 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 [2007.00789].

## 4. Computational Complexity and Storage

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

| Algorithm         | Factorization Time             | Memory       | Comments                         |
|-------------------|-------------------------------|--------------|-----------------------------------|
| Direct Cholesky   | $O(n^3)$                      | $O(n^2)$     | Prohibitive for large $n$         |
| Low-rank FFBS     | $O(T n N^2)$                  | $O(n N)$     | Can miss fine scale dependence    |
| Randomized/Hierarchical Cholesky | $O(n r^2 \log n)$ [1507.05593], $O(T n N^2)$ [2207.09384] | $O(n r \log n)$ [1507.05593], $O(n N)$ [2207.09384] | 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$ ensures compatibility with distributed and parallel computing architectures. Reporting from [1507.05593] demonstrates memory reductions of over an order of magnitude versus direct Cholesky (e.g., $0.96$ GB vs $7.6$ GB for $n \approx 4\times 10^5$), 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 [2207.09384][2006.16901]. 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 $O(n N)$ storage and $O(n N^2 T)$ total computation.

Empirical results show that modest neighborhood sizes (e.g., $N=40$–$80$) retain predictive accuracy within a few percent of the full dense FFBS (as measured by CRPS), with typical speedups of $20\times$ and superior accuracy to low-rank approximations at comparable computational cost [2207.09384].

## 6. Advanced Error Bounds and Theoretical Guarantees

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

- Given any chosen sparsity mask $S$, the hierarchical sparse Cholesky factorization yields the optimal KL approximation to $N(0, \Sigma)$ within the class of distributions induced by sparse Cholesky factors with pattern $S$ [2207.09384].
- For hierarchical Vecchia patterns, error in the log-determinant or inverse-trace decays algebraically or exponentially as the neighborhood size $N$ increases, under standard regularity (e.g., Matérn covariance smoothness) [2207.09384][2006.16901].
- In randomized low-rank block compression, the norm of the error introduced at each interface block is $O(\varepsilon^2)$ (second-order scheme) and the overall impact on PCG iteration counts is a halving compared to first-order truncations [2007.00789].

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 [1507.05593].

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 [1507.05593].

---

**Key references:** "An Efficient Solver for Sparse Linear Systems Based on Rank-Structured Cholesky Factorization" [1507.05593], "Second Order Accurate Hierarchical Approximate Factorization of Sparse SPD Matrices" [2007.00789], "Scalable Spatio-Temporal Smoothing via Hierarchical Sparse Cholesky Decomposition" [2207.09384], and "Hierarchical sparse Cholesky decomposition with applications to high-dimensional spatio-temporal filtering" [2006.16901].

Source: https://www.emergentmind.com/topics/randomized-cholesky-rchol