---
title: Pivoted Cholesky Algorithm
url: https://www.emergentmind.com/topics/pivoted-cholesky-algorithm
type: topic
---

# Pivoted Cholesky Algorithm

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

## 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
\[
L_{p_t,t}=\sqrt{R_t(p_t,p_t)},\qquad
L_{i,t}=\frac{R_t(i,p_t)}{L_{p_t,t}}.
\]
Equivalently, one may work with a permutation matrix \(P\) and write
\[
PAP^\top = LL^\top,
\qquad\text{or}\qquad
A\approx P^\top L(:,1\!:\!r)L(:,1\!:\!r)^\top P.
\]
The adaptive stopping rule is based on the largest residual diagonal: the factorization terminates when that quantity falls below a prescribed tolerance [2601.03706; 1411.6296].

A central invariant is the diagonal residual vector
\[
r_i^{(k-1)} = A_{ii}-\sum_{j=1}^{k-1}L_{ij}^2.
\]
After forming column \(k\), the update is simply
\[
r_i^{(k)} = r_i^{(k-1)} - L_{i,k}^2.
\]
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 [1411.6296].

## 2. RKHS geometry and orthogonalization

For kernel matrices, pivoted Cholesky admits a direct geometric interpretation in the reproducing kernel Hilbert space. Let \(\phi:X\to\mathcal H\) be the feature map with
\[
K(i,j)=\langle \phi(x_i),\phi(x_j)\rangle_{\mathcal H}.
\]
Then the residual diagonal satisfies
\[
R_t(i,i)=K(i,i)-\sum_{k<t}L_{i,k}^2
\]
and is exactly the squared distance from \(\phi(x_i)\) to the span of the previously selected feature vectors. In the notation of the geometric derivation,
\[
R_t(i,i)=\|\phi(x_i)-\operatorname{Proj}_{S_{t-1}}\phi(x_i)\|^2,
\]
where
\[
S_{t-1}=\operatorname{Span}\{\phi(x_{p_1}),\dots,\phi(x_{p_{t-1}})\}.
\]
Maximizing \(R_t(i,i)\) therefore chooses the point whose feature is farthest from the current subspace in \(\mathcal H\). This is precisely greedy farthest point sampling in the RKHS metric [2601.03706].

The same construction can be viewed as an implicit Gram–Schmidt process. At step \(t\), one forms the residual of the pivot feature,
\[
r_t=\phi(x_{p_t})-\sum_{k<t}\langle \phi(x_{p_t}),e_k\rangle e_k,
\]
normalizes it to obtain
\[
e_t=\frac{r_t}{\|r_t\|},
\]
and then computes
\[
L_{i,t}=\langle \phi(x_i),e_t\rangle
\]
for all \(i\). Stacking the vectors \(e_k\) yields an orthonormal basis \(Q\) in feature space, while the Cholesky factor is the matrix of inner products \(Q^\top\Phi\). In this sense, pivoted Cholesky is not merely a matrix update scheme; it is an orthogonalization process carried out implicitly in feature space [2601.03706].

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 \(O(n^2)\) work per iteration. The standard “lazy” variant avoids this by storing only the diagonal residuals
\[
d_i \equiv R_t(i,i),
\]
together with the previously computed factor columns. At iteration \(t\), it selects the pivot from the current diagonal residuals, computes the new column through
\[
c_i = K(i,p_t)-\sum_{k<t}L_{i,k}L_{t,k},
\qquad
L_{i,t}=c_i/L_{t,t},
\]
and updates
\[
d_i \leftarrow d_i - L_{i,t}^2.
\]
The resulting cost is \(O(nM)\) kernel evaluations, \(O(nM^2)\) arithmetic, and \(O(nM)\) memory, with no need to store the full \(n\times n\) matrix [2601.03706].

A closely related left-looking formulation requires only the current diagonal residuals and the new pivot column \(A(:,p)\). In applications where entries are expensive to compute, one evaluates \(A_{i,p}\) on demand and never forms the full trailing submatrix. After \(k\) steps, only \(O(Nk)\) entries of \(A\) are evaluated rather than \(O(N^2)\) [1411.6296]. 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 \(O(k^2 n)\) cost and sharp \(\|\cdot\|_\infty\) error control for SPSD matrices [1505.06195].

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 \(2.5\)–\(3\)\(\times\) for full-rank factorizations and about \(1.8\)–\(2.0\)\(\times\) in the low-rank, lazy-evaluation case [1411.6296].

## 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_{jj}-(LL^\top)_{jj}=\operatorname{Var}(f(x_j)\mid f(x_I)),
\]
so the standard rule is a greedy entropy maximization or D-optimal selection strategy [2507.20678].

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 |
| \(\beta=2\) 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 |
| \(\Delta\)-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+1)N\) entry evaluations, \(O(k^2N)\) arithmetic, and \(O(kN)\) storage [2207.06503]. A different randomized rule selects pivots with probability proportional to \(A_{ii}^2\), yielding a per-step Frobenius-norm contraction guarantee analogous to the trace-norm guarantee obtained when sampling proportional to \(A_{ii}\) [2404.11487].

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

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 \(10\)–\(30\%\) fewer CG iterations than variance pivoting, and that they reach the same negative log-marginal-likelihood with \(20\)–\(30\%\) fewer inducing points, with extra overhead below \(5\%\) of total training time [2507.20678].

A more recent development is \(\Delta\)-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 \(\mathcal O(snr^2)\) for an \(n\times n\) matrix and target rank \(r\). The method guarantees monotonically non-decreasing functional values and, among Cholesky-consistent rank-1 updates, maximizes the per-step gain [2606.01821].

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 \(40\times\) faster, and that it precisely matches simple RPCholesky in accuracy [2410.03969].

## 5. Error analysis and convergence theory

Several complementary error viewpoints exist. In diagonal-pivoted cross approximation, the residual
\[
R=\Phi-LL^\top
\]
remains SPSD, and its maximum entrywise error is exactly the next pivot value:
\[
\|R\|_\infty=\max_i R_{ii}=\gamma_{\ell+1}.
\]
If the spectrum decays exponentially, \(\sigma_k(\Phi)=O(e^{-\alpha k})\), then the diagonal-pivoted PCD error is \(O(e^{-\alpha \ell})\) [1505.06195].

For RPCholesky, the central result is near-optimality in expected trace-norm error. If \(A_r\) is an optimal rank-\(r\) psd approximation of \(A\), \(\eta=\operatorname{tr}(A-A_r)/\operatorname{tr}(A)\), and
\[
k\ge \frac{r}{\varepsilon}+r\log\!\Bigl(\frac{1}{\varepsilon\eta}\Bigr),
\]
then
\[
\mathbb E\,\operatorname{tr}(A-A^{(k)})\le (1+\varepsilon)\operatorname{tr}(A-A_r).
\]
The same paper presents an \(\eta\)-free variant and describes the approximation as nearly optimal [2207.06503].

Randomized partial Cholesky also admits simple one-step contraction laws. If the pivot is sampled with probability proportional to \(A_{ii}\), then
\[
\mathbb E\|R^{(t)}\|_{S_1}\le (1-1/n)^t\|A\|_{S_1}.
\]
If the pivot is sampled with probability proportional to \(A_{ii}^2\), then
\[
\mathbb E\|R^{(t)}\|_F^2\le (1-1/n)^t\|A\|_F^2.
\]
These results isolate trace-norm and Frobenius-norm contractivity as properties of the sampling law rather than of a single deterministic pivot rule [2404.11487].

A continuous convergence theory is available for kernels. For a symmetric positive definite Lipschitz continuous kernel \(K:\Omega\times\Omega\to\mathbb R\) on compact \(\Omega\subset\mathbb R^d\), 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,
\[
\|R_n\|_\infty = O(n^{-1/d}),
\]
and if \(K\) is \(C^{1,1}\), the rate improves to
\[
\|R_n\|_\infty = O(n^{-2/d}).
\]
The same framework identifies the residual with the Gaussian-process posterior covariance and the diagonal residual with the squared power function in \(P\)-greedy interpolation [2509.13582].

## 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 [2601.03706]. 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 [2306.03955].

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 \(\Delta\)-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 [2507.20678; 2606.01821].

In uncertainty quantification, diagonal-pivoted CA has been used to approximate a covariance matrix \(C\in\mathbb R^{56312\times 56312}\) with rank \(k=600\), followed by QR and a small SVD to produce an eigendecomposition. The reported timings are about \(12\,\mathrm{s}\) for QR on the low-rank factor and about \(47\,\mathrm{s}\) for the reduced SVD, compared with about \(38\,\mathrm{h}\) for direct dense eigendecomposition [1505.06195].

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 \(80000\) atomic orbitals, and the partitioned variant reduces peak memory by roughly \(1/K\) [1811.12890]. A separate application prunes overcomplete atomic-orbital basis sets by pivoted Cholesky of the overlap matrix; with \(\tau=10^{-6}\), the reported reductions are about \(9\%\) fewer functions in single-\(\zeta\), about \(21\%\) fewer in double-\(\zeta\), and about \(28\%\) fewer in triple-\(\zeta\) sets, while restoring numerical stability in problematic cases [1911.10372].

Pivoted Cholesky also interfaces with sparse inverse-Cholesky methods. One recent result shows that a partial pivoted Cholesky approximation of rank \(r\), 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 [2603.05709].

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.

Source: https://www.emergentmind.com/topics/pivoted-cholesky-algorithm