Pivoted Cholesky Algorithm

• Pivoted Cholesky is a rank-revealing factorization that sequentially selects pivots based on the largest residual diagonal to construct a low-rank approximation.
• It admits an RKHS interpretation, where the pivot rule corresponds to greedy farthest point sampling and serves as an implicit Gram–Schmidt orthogonalization.
• Efficient implementations use diagonal updates and on-demand kernel evaluations, making the method practical for large-scale Gaussian process applications.

Searching arXiv for recent and foundational papers on pivoted Cholesky to ground the article. Pivoted Cholesky is a rank-revealing Cholesky-type factorization for symmetric positive semidefinite or positive definite matrices that constructs a low-rank approximation by selecting pivots sequentially from the residual diagonal. In kernel methods, it is a standard tool for approximating large kernel matrices without forming or storing the full matrix, and it is especially prominent in large-scale Gaussian process computation. Recent work has sharpened both its algebraic and geometric interpretations, showing that, in an RKHS, the pivot rule coincides with greedy farthest point sampling in the kernel metric and that the factor construction is an implicit Gram–Schmidt orthogonalization (Shabat, 7 Jan 2026).

1. Algebraic formulation

Let $A\in\mathbb R^{N\times N}$ be symmetric positive semidefinite, or let $K\in\mathbb R^{n\times n}$ be a symmetric positive-definite kernel matrix. The objective is a low-rank factorization

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

or, equivalently,

$$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$

The factor is built column by column. If $L_{:,k}$ denotes the $k$-th column, the residual after $t-1$ steps is

$R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$

with entrywise form

$R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$

The classical pivot rule chooses the next pivot index by maximizing the diagonal of the current residual: $p_t=\arg\max_i R_t(i,i).$ After the pivot is selected and conceptually swapped into the active position, the update is

$K\in\mathbb R^{n\times n}$0

Equivalently, one may work with a permutation matrix $K\in\mathbb R^{n\times n}$1 and write

$K\in\mathbb R^{n\times n}$2

The adaptive stopping rule is based on the largest residual diagonal: the factorization terminates when that quantity falls below a prescribed tolerance (Shabat, 7 Jan 2026, Loan et al., 2014).

A central invariant is the diagonal residual vector

$K\in\mathbb R^{n\times n}$3

After forming column $K\in\mathbb R^{n\times n}$4, the update is simply

$K\in\mathbb R^{n\times n}$5

This diagonal recursion is the basis of most efficient implementations. It also underlies the rank-revealing character of the algorithm: the selected pivots are precisely those indices whose unexplained diagonal mass is largest at each step (Loan et al., 2014).

2. RKHS geometry and orthogonalization

For kernel matrices, pivoted Cholesky admits a direct geometric interpretation in the reproducing kernel Hilbert space. Let $K\in\mathbb R^{n\times n}$6 be the feature map with

$K\in\mathbb R^{n\times n}$7

Then the residual diagonal satisfies

$K\in\mathbb R^{n\times n}$8

and is exactly the squared distance from $K\in\mathbb R^{n\times n}$9 to the span of the previously selected feature vectors. In the notation of the geometric derivation,

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

where

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

Maximizing $$A \approx L L^\top,\qquad L\in\mathbb R^{N\times r},\; r\ll N,$$2 therefore chooses the point whose feature is farthest from the current subspace in $$A \approx L L^\top,\qquad L\in\mathbb R^{N\times r},\; r\ll N,$$3. This is precisely greedy farthest point sampling in the RKHS metric (Shabat, 7 Jan 2026).

The same construction can be viewed as an implicit Gram–Schmidt process. At step $$A \approx L L^\top,\qquad L\in\mathbb R^{N\times r},\; r\ll N,$$4, one forms the residual of the pivot feature,

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

normalizes it to obtain

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

and then computes

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

for all $$A \approx L L^\top,\qquad L\in\mathbb R^{N\times r},\; r\ll N,$$8. Stacking the vectors $$A \approx L L^\top,\qquad L\in\mathbb R^{N\times r},\; r\ll N,$$9 yields an orthonormal basis $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$0 in feature space, while the Cholesky factor is the matrix of inner products $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$1. In this sense, pivoted Cholesky is not merely a matrix update scheme; it is an orthogonalization process carried out implicitly in feature space (Shabat, 7 Jan 2026).

This interpretation clarifies why the algorithm aligns naturally with greedy sampling, kernel interpolation, and posterior variance reduction. A plausible implication is that many apparently different kernel algorithms are rephrasings of the same residual-subspace geometry.

3. Efficient implementations and structure-preserving variants

A naïve implementation would recompute the full residual matrix at each step, incurring $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$2 work per iteration. The standard “lazy” variant avoids this by storing only the diagonal residuals

$$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$3

together with the previously computed factor columns. At iteration $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$4, it selects the pivot from the current diagonal residuals, computes the new column through

$$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$5

and updates

$$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$6

The resulting cost is $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$7 kernel evaluations, $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$8 arithmetic, and $$K \approx L L^\top,\qquad L\in\mathbb R^{n\times M},\; M\ll n.$$9 memory, with no need to store the full $L_{:,k}$0 matrix (Shabat, 7 Jan 2026).

A closely related left-looking formulation requires only the current diagonal residuals and the new pivot column $L_{:,k}$1. In applications where entries are expensive to compute, one evaluates $L_{:,k}$2 on demand and never forms the full trailing submatrix. After $L_{:,k}$3 steps, only $L_{:,k}$4 entries of $L_{:,k}$5 are evaluated rather than $L_{:,k}$6 (Loan et al., 2014). This on-demand structure is particularly important when matrix entries arise from costly integral evaluations.

Diagonal-pivoted cross approximation provides another formulation of the same basic mechanism. In that setting, one maintains the diagonal of the residual, selects the largest diagonal pivot, forms the new column, and updates the diagonal. The tuned diagonal-pivoted CA algorithm attains $L_{:,k}$7 cost and sharp $L_{:,k}$8 error control for SPSD matrices (Liu et al., 2015).

When additional matrix symmetry is present, pivoted Cholesky can be made structure preserving. For matrices that are both symmetric and perfect-shuffle symmetric, a closed-form block diagonalization reduces the problem to one or two half-sized pivoted Cholesky factorizations. The paper reports a theoretical four-fold reduction in flop count, and also states that observed speedups are about $L_{:,k}$9–%%%%3$A\in\mathbb R^{N\times N}$3%%%%1 for full-rank factorizations and about $k$2–$k$3$k$4 in the low-rank, lazy-evaluation case (Loan et al., 2014).

4. Pivot rules and selection objectives

The classical pivot rule is greedy diagonal maximization. In Gaussian-process language, this corresponds to maximizing the remaining Schur-complement diagonal,

$k$5

so the standard rule is a greedy entropy maximization or D-optimal selection strategy (Roos et al., 28 Jul 2025).

Subsequent work has emphasized that pivot selection need not be tied to this single criterion.

Rule	Selection criterion	Stated target
Greedy diagonal pivoting	Largest residual diagonal	Entropy / determinant reduction
RPCholesky	Probability proportional to residual diagonal	Trace-norm contraction in expectation
$k$6 randomized rule	Probability proportional to squared diagonal	Frobenius-norm contraction in expectation
PCov / WPCov	Surrogates based on projected covariance or residual-weighted covariance	Sparse GP regression and preconditioning
$k$7-VFE	Exact one-step gain in variational free energy over sampled candidates	GP variational free energy

Randomly pivoted Cholesky samples the next pivot with probability proportional to the current residual diagonal. In its standard matrix form, it requires only $k$8 entry evaluations, $k$9 arithmetic, and $t-1$0 storage (Chen et al., 2022). A different randomized rule selects pivots with probability proportional to $t-1$1, yielding a per-step Frobenius-norm contraction guarantee analogous to the trace-norm guarantee obtained when sampling proportional to $t-1$2 (Steinerberger, 2024).

These randomized results also delimit a common misconception. Greedy largest-diagonal pivoting is not universally dominant: one paper notes that greedy can almost stall or pick “wrong” pivots when the diagonal is not representative of column norm, while randomized rules trade a deterministic one-step choice for robust contraction in expectation (Steinerberger, 2024).

For Gaussian-process workloads, alternative deterministic objectives have been proposed. PCov uses a trace-oriented projected covariance surrogate, and WPCov incorporates observations through a residual-weighted projected covariance criterion. The paper reports that PCov and WPCov typically require $t-1$3–$t-1$4 fewer CG iterations than variance pivoting, and that they reach the same negative log-marginal-likelihood with $t-1$5–$t-1$6 fewer inducing points, with extra overhead below $t-1$7 of total training time (Roos et al., 28 Jul 2025).

A more recent development is $t-1$8-VFE pivoted Cholesky, which maximizes the exact one-step change of the variational free energy under Cholesky-consistent rank-1 updates. Candidate pivots are sampled in batches proportional to the residual diagonal, and the total cost is $t-1$9 for an $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$0 matrix and target rank $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$1. The method guarantees monotonically non-decreasing functional values and, among Cholesky-consistent rank-1 updates, maximizes the per-step gain (Schaub et al., 1 Jun 2026).

For randomized methods, blocking and rejection sampling substantially alter runtime without changing the underlying approximation law. Accelerated RPCholesky uses block computations and rejection sampling to simulate the original RPCholesky while exploiting submatrix access; for kernel matrix approximation, the paper states that it can run over $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$2 faster, and that it precisely matches simple RPCholesky in accuracy (Epperly et al., 2024).

5. Error analysis and convergence theory

Several complementary error viewpoints exist. In diagonal-pivoted cross approximation, the residual

$R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$3

remains SPSD, and its maximum entrywise error is exactly the next pivot value: $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$4 If the spectrum decays exponentially, $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$5, then the diagonal-pivoted PCD error is $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$6 (Liu et al., 2015).

For RPCholesky, the central result is near-optimality in expected trace-norm error. If $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$7 is an optimal rank-$R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$8 psd approximation of $R_t = K - \sum_{k=1}^{t-1} L_{:,k}L_{:,k}^\top,$9, $R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$0, and

$R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$1

then

$R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$2

The same paper presents an $R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$3-free variant and describes the approximation as nearly optimal (Chen et al., 2022).

Randomized partial Cholesky also admits simple one-step contraction laws. If the pivot is sampled with probability proportional to $R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$4, then

$R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$5

If the pivot is sampled with probability proportional to $R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$6, then

$R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$7

These results isolate trace-norm and Frobenius-norm contractivity as properties of the sampling law rather than of a single deterministic pivot rule (Steinerberger, 2024).

A continuous convergence theory is available for kernels. For a symmetric positive definite Lipschitz continuous kernel $R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$8 on compact $R_t(i,j)=K(i,j)-\sum_{k=1}^{t-1}L_{i,k}L_{j,k}.$9, the residual of the continuous pivoted Cholesky algorithm with any pivoting strategy is bounded by a constant multiple of the fill distance of the pivots. Under complete pivoting,

$p_t=\arg\max_i R_t(i,i).$0

and if $p_t=\arg\max_i R_t(i,i).$1 is $p_t=\arg\max_i R_t(i,i).$2, the rate improves to

$p_t=\arg\max_i R_t(i,i).$3

The same framework identifies the residual with the Gaussian-process posterior covariance and the diagonal residual with the squared power function in $p_t=\arg\max_i R_t(i,i).$4-greedy interpolation (Jeong et al., 16 Sep 2025).

6. Applications, extensions, and adjacent formalisms

The most visible contemporary use of pivoted Cholesky is large-scale kernel approximation. Low-rank approximations of large kernel matrices are described as ubiquitous in machine learning, particularly for scaling Gaussian processes to massive datasets (Shabat, 7 Jan 2026). In kernel quadrature, randomly pivoted Cholesky yields node sets that achieve quadrature error rates comparable to more expensive schemes based on continuous volume sampling, thinning, and recombination, while remaining Mercer-free and adaptable to complicated geometries with arbitrary kernels (Epperly et al., 2023).

In Gaussian-process inference, pivoted Cholesky is used both as a sparse inducing-point mechanism and as a preconditioner for iterative solvers. The GP-specific selection rules PCov, WPCov, and $p_t=\arg\max_i R_t(i,i).$5-VFE shift the emphasis from residual trace alone to nonlinear objectives such as posterior uncertainty, trace penalties, data fit, log-determinants, and variational free energy (Roos et al., 28 Jul 2025, Schaub et al., 1 Jun 2026).

In uncertainty quantification, diagonal-pivoted CA has been used to approximate a covariance matrix $p_t=\arg\max_i R_t(i,i).$6 with rank $p_t=\arg\max_i R_t(i,i).$7, followed by QR and a small SVD to produce an eigendecomposition. The reported timings are about $p_t=\arg\max_i R_t(i,i).$8 for QR on the low-rank factor and about $p_t=\arg\max_i R_t(i,i).$9 for the reduced SVD, compared with about $K\in\mathbb R^{n\times n}$00 for direct dense eigendecomposition (Liu et al., 2015).

In quantum chemistry, pivoted Cholesky has been adapted to electron repulsion integrals by separating pivot selection from vector construction, using aggressive screening and optional partitioning. The reported algorithm is applied to systems with up to $K\in\mathbb R^{n\times n}$01 atomic orbitals, and the partitioned variant reduces peak memory by roughly $K\in\mathbb R^{n\times n}$02 (Folkestad et al., 2018). A separate application prunes overcomplete atomic-orbital basis sets by pivoted Cholesky of the overlap matrix; with $K\in\mathbb R^{n\times n}$03, the reported reductions are about $K\in\mathbb R^{n\times n}$04 fewer functions in single-$K\in\mathbb R^{n\times n}$05, about $K\in\mathbb R^{n\times n}$06 fewer in double-$K\in\mathbb R^{n\times n}$07, and about $K\in\mathbb R^{n\times n}$08 fewer in triple-$K\in\mathbb R^{n\times n}$09 sets, while restoring numerical stability in problematic cases (Lehtola, 2019).

Pivoted Cholesky also interfaces with sparse inverse-Cholesky methods. One recent result shows that a partial pivoted Cholesky approximation of rank $K\in\mathbb R^{n\times n}$10, combined with a Vecchia approximation of the residual, is exactly a Vecchia approximation of the original matrix with an augmented sparsity pattern. This places partial Cholesky within a broader inverse-factor framework rather than treating it as an isolated low-rank technique (Kaminetz et al., 5 Mar 2026).

Taken together, these developments position pivoted Cholesky as a family of residual-driven factorizations rather than a single algorithmic template. Its classical form is a deterministic, diagonal-greedy, low-rank Cholesky construction; its modern variants reinterpret pivoting geometrically, randomize it for norm-specific guarantees, adapt it to nonlinear objectives, and fuse it with structure-preserving, sparse, or domain-specific computational strategies.