Lazy Pivoted Cholesky

• Lazy Pivoted Cholesky is a method for computing low-rank approximations of SPSD matrices by selectively evaluating only essential matrix elements, which reduces computational and memory costs.
• It employs a greedy pivot selection based on current residual diagonals, enabling efficient iterative updates and ensuring controlled uniform error bounds.
• The approach integrates techniques like farthest-point sampling and Gram-Schmidt orthogonalization, leading to significant runtime improvements in applications such as kernel methods and quantum chemistry.

Lazy Pivoted Cholesky refers to a family of algorithms for computing low-rank approximations of large symmetric positive semidefinite (SPSD) matrices, where only a subset of matrix elements—typically the diagonals and selected columns—are computed or updated at each iteration, often on demand. This approach enables efficient decomposition both in computational time and memory, particularly for matrices arising from kernel methods, quantum chemistry electron repulsion integrals, or large-scale probabilistic models. Lazy Pivoted Cholesky merges concepts from pivoted Cholesky, maximal-volume cross approximation, and greedy farthest-point sampling, with the distinctive feature of avoiding full matrix access except when strictly required.

1. Mathematical Formulation and Decomposition Objective

Let $A \in \mathbb{R}^{n \times n}$ denote an SPSD matrix, which may represent, for example, electron repulsion integrals in quantum chemistry, kernel matrices in machine learning, or covariance matrices in Gaussian process modeling. The goal is to construct a low-rank decomposition,

$$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$

such that the residual $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$ does not exceed a prespecified threshold $\tau > 0$ (Folkestad et al., 2018, Liu et al., 2015, Shabat, 7 Jan 2026).

A defining property of the lazy variant is that it advances the decomposition relying only on the current residual diagonal, or a restricted set of matrix elements, thereby minimizing matrix evaluations and storage. For kernel matrices $K_{ij}=k(x_i, x_j)$, the decomposition operates strictly via kernel evaluations without constructing $K$ in full (Shabat, 7 Jan 2026). Similarly, in quantum chemistry, only a small subset of electron repulsion integrals is needed for each step (Folkestad et al., 2018).

2. Algorithmic Structure and Lazy Pivot Selection

The lazy pivoted Cholesky method proceeds in iterative stages:

1. Diagonal Initialization: Initialize a residual diagonal $d_i = A_{ii}$ for all $i$.
2. Greedy Pivot Selection: At each iteration, select the pivot index $i^* = \arg\max_{i} d_i$, i.e., the index with the largest residual diagonal entry. Only these diagonals are updated after each step. In practical lazy implementations, a max-heap or priority queue over $d_i$ is maintained so that pivot selection operates in $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$0 time, and only a subset of $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$1 are updated at each step (Kaminetz et al., 5 Mar 2026, Hemmi et al., 2022).
3. Rank-1 Schur Complement Update: For each non-pivoted index $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$2, update only the necessary row/column of $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$3 or compute $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$4 as required. For the next Cholesky column $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$5 with pivot $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$6:

$$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$7

Then update $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$8 (Liu et al., 2015, Folkestad et al., 2018).

1. Termination: The process repeats until the maximal diagonal residual falls below the threshold $$A \approx L L^\top, \qquad L \in \mathbb{R}^{n \times r}, \qquad r \ll n$$9.

The algorithm may be adapted so that, at each iteration, only a small batch of candidate indices (with residuals close to the current maximum, controlled by a span parameter $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$0) is considered to further reduce work (Folkestad et al., 2018).

A "partitioned" variant splits the active diagonal index set into batches, applies the procedure within each batch, and then merges (Folkestad et al., 2018).

3. Geometric and Theoretical Underpinnings

The lazy selection of pivots has a rigorous geometric and algebraic interpretation. In the kernel and RKHS setting, each step of the algorithm corresponds to greedy farthest-point sampling (FPS): at each iteration, the point with the largest squared distance from the span of previously selected features is chosen as the next pivot. The update

$\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$1

is the squared residual in the RKHS, and the algorithm precisely implements greedy FPS according to the kernel-induced metric (Shabat, 7 Jan 2026).

Additionally, the Cholesky factor construction via the Schur-complement update implicitly performs Gram-Schmidt orthogonalization on the selected feature vectors. Formally, the process is equivalent to a column-pivoted QR of the (generally infinite-dimensional) feature matrix $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$2 (Shabat, 7 Jan 2026).

For cross-approximation of SPSD matrices, the maximal diagonal pivot at each iteration is justified by the property that all off-diagonal entries are bounded by the geometric mean of the diagonal entries due to SPSD structure. This ensures that the pivoted column yields the maximal possible update in the Schur complement and underpins the uniform entrywise error bounds (Liu et al., 2015).

4. Complexity, Error Bounds, and Practical Efficiency

Lazy Pivoted Cholesky is designed for computational scalability:

• Time Complexity:

For $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$3-rank approximation and $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$4-dimensional matrices, the cost per iteration is typically $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$5 for diagonal updates and column construction, giving total complexity $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$6 in basic lazy implementations (Liu et al., 2015). If a max-heap or approximate priority queue is used for pivot selection and only a small subset of diagonals are updated, the per-pivot cost may be reduced to $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$7, yielding near $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$8 overall (Kaminetz et al., 5 Mar 2026, Hemmi et al., 2022). For kernel matrices, the cost is $\max_{p, q} |A_{pq} - (LL^\top)_{pq}|$9 for $\tau > 0$0 pivots in an $\tau > 0$1-point dataset (Shabat, 7 Jan 2026).

• Memory Usage:

Only $\tau > 0$2 entries of $\tau > 0$3 are stored, with no need to form the full $\tau > 0$4 matrix. In large datasets, the peak memory can be reduced further with partitioning and screening (Folkestad et al., 2018).

• Error Guarantees:

If the process is terminated at a threshold $\tau > 0$5 (i.e., all diagonal residuals $\tau > 0$6), the uniform norm $\tau > 0$7, and the error on all off-diagonals is at most $\tau > 0$8 via the Cauchy-Schwarz bound. For trace and spectral norm bounds, the residual is controlled as

$\tau > 0$9

where $K_{ij}=k(x_i, x_j)$0 is the $K_{ij}=k(x_i, x_j)$1-th eigenvalue of $K_{ij}=k(x_i, x_j)$2 (Kaminetz et al., 5 Mar 2026).

• Empirical Performance:

On quantum chemistry matrices (e.g., $K_{ij}=k(x_i, x_j)$3 basis functions), lazy Cholesky matches conventional pivoted Cholesky in accuracy ($K_{ij}=k(x_i, x_j)$4) but achieves over an order of magnitude improvement in runtime and can reduce peak memory by factors of 3–6 using partitioning (Folkestad et al., 2018). In kernel learning and DPP MAP inference, lazy pivoted Cholesky (often termed lazy+fast greedy) ran up to 15–20× faster than non-lazy alternatives with identical outputs (Hemmi et al., 2022).

5. Extensions, Hybrid Schemes, and Structure Exploitation

Lazy Pivoted Cholesky is extended and adapted to various structured domains:

• Partitioned and Two-Phase Schemes:

For matrices decomposed into multiple blocks (e.g., batches of orbitals in quantum chemistry), partitioned lazy pivoting followed by a merge phase is used to further minimize memory and computation (Folkestad et al., 2018).

• Structure-Preserving Factorization:

For matrices with additional symmetries, such as perfect-shuffle symmetry (as in certain tensor unfoldings for electron integrals), lazy pivoted Cholesky is performed after explicit block-diagonalization, often halving the size of the matrix to be factorized (Loan et al., 2014). This reduces the cost by up to a factor of 4, with approximately proportional reductions in storage and required matrix evaluations.

• Hybrid with Sparse Approximation (Vecchia):

The residual after lazy partial Cholesky may itself be approximated by a sparse inverse-Cholesky (Vecchia) factorization. The sum then constitutes a generalized Vecchia approximation with improved global (low-rank) and local (sparse) representation (Kaminetz et al., 5 Mar 2026).

• Kernel and MAP Inference Applications:

For large-scale kernel methods and determinantal point process (DPP) MAP inference, the lazy approach provides a geometric (RKHS-based) and algorithmic foundation for submodular maximization, leveraging properties of farthest-point sampling and submodularity of the log-determinant (Shabat, 7 Jan 2026, Hemmi et al., 2022).

6. Implementation and Best Practices

The minimalist implementation requires:

• Maintenance of a (possibly partial) residual diagonal,
• On-demand evaluation of matrix columns/entries,
• Update of basis coefficients and diagonals per the Cholesky recursion,
• Termination when the prescribed error threshold is met.

Python pseudocode for kernel decomposition is given in (Shabat, 7 Jan 2026), while a rigorous two-phase basis-selection and projection procedure is described in quantum chemistry contexts (Folkestad et al., 2018).

In practice, performance depends on the efficacy of diagonal screening, correctness of heap or priority queue management, and the sparsity and decay properties of the spectrum of $K_{ij}=k(x_i, x_j)$5. For matrices with fast eigenvalue decay, convergence to high accuracy is typically achieved in $K_{ij}=k(x_i, x_j)$6 iterations (Liu et al., 2015).

7. Limitations and Comparative Context

The primary limitations of the method stem from the non-optimality of the greedy pivot sequence: lazy heuristics (partial updates, approximate heaps) may infrequently select suboptimal pivots, slightly degrading spectral tail bounds. However, for most practical purposes, the error guarantees and runtime efficiency make lazy pivoted Cholesky an attractive approach for large SPSD matrices (Kaminetz et al., 5 Mar 2026).

Relative to random sampling/Nystrom-type approximations, lazy pivoted Cholesky provides deterministic uniform-norm control and automatic rank-adaptation, at the cost of more involved pivot management. In problems where matrix entries are especially costly to compute, such as high-order integral evaluations, the lazy evaluation paradigm delivers significant savings by restricting access to the minimal required subset (Loan et al., 2014, Folkestad et al., 2018).

A hybrid approach with Vecchia greatly broadens the class of matrices that can be accurately approximated in practice, uniting low-rank and sparsity-exploiting representations (Kaminetz et al., 5 Mar 2026).