Pivoted Cholesky Preconditioning

• The paper introduces a method that approximates an SPD matrix with a low-rank Cholesky factorization, enabling rapid convergence of iterative solvers.
• It employs greedy, randomized, and hybrid pivoting strategies to control error bounds and reduce computational complexity for large kernel matrices.
• The technique enhances scalability in Gaussian process regression and kernel methods by minimizing factorization cost and improving convergence properties.

Pivoted Cholesky preconditioning refers to the use of low-rank, diagonally-pivoted Cholesky decompositions as preconditioners for solving large, often ill-conditioned, symmetric positive-definite (SPD) linear systems. This technique is particularly significant in computational linear algebra, Gaussian process regression, kernel methods, and scalable numerical optimization, where the cost of direct factorization is prohibitive and rapid convergence of iterative solvers is essential.

1. Pivoted Cholesky Factorization: Algorithms and Variants

Pivoted Cholesky factorization constructs a low-rank approximation $A \approx LL^T$, where $A \in \mathbb{R}^{n \times n}$ is SPD, and $L \in \mathbb{R}^{n \times k}$ for $k \ll n$. Unlike the standard Cholesky, the pivoted version selects at each iteration the "pivot" index corresponding to the maximal (or, in randomized, high-probability) diagonal entry of the residual matrix.

Algorithmic principle:

Given $A^{(0)} = A$, for $i = 1, \dots, k$:

• Select $s_i$ according to the pivoting strategy (e.g., $s_i = \arg\max_j A^{(i-1)}_{jj}$ for greedy).
• Define the Cholesky vector:

$$\ell_i = \frac{A^{(i-1)}_{:,s_i}}{\sqrt{A^{(i-1)}_{s_i, s_i}}}$$

• Update the factor:

$$L_i = [\ell_1, \ldots, \ell_i], \qquad A^{(i)} = A^{(i-1)} - \ell_i \ell_i^T$$

• Continue until a rank $A \in \mathbb{R}^{n \times n}$0 approximation or the desired tolerance is achieved.

Pivot selection strategies:

• Greedy (deterministic): select the largest diagonal entry.
• Randomized: sample $A \in \mathbb{R}^{n \times n}$1 with probability proportional to $A \in \mathbb{R}^{n \times n}$2 with $A \in \mathbb{R}^{n \times n}$3 (trace-norm or Frobenius-norm contraction) (Steinerberger, 2024).
• Complete pivoting in the context of kernel functions selects the spatial location maximizing the current diagonal of the residual kernel (Jeong et al., 16 Sep 2025).

Hybrid methods such as cross approximation (CA) merge with diagonal pivoting in the kernel setting and have been demonstrated to be equivalent for SPD matrices (Liu et al., 2015).

Computational complexity:

The classic implementation costs $A \in \mathbb{R}^{n \times n}$4 or $A \in \mathbb{R}^{n \times n}$5 depending on whether column or row updating is employed. Using efficient data structures, such as binary trees or alias tables, pivot selection and diagonal updates can be performed in $A \in \mathbb{R}^{n \times n}$6 or even $A \in \mathbb{R}^{n \times n}$7 amortized time for randomized rules (Steinerberger, 2024).

2. Theoretical Guarantees and Convergence Rates

Norm contraction:

• For $A \in \mathbb{R}^{n \times n}$8 (trace-norm pivoting), at each step the expected reduction in $A \in \mathbb{R}^{n \times n}$9 is a universal contraction:

$L \in \mathbb{R}^{n \times k}$0

• For $L \in \mathbb{R}^{n \times k}$1 (Frobenius-norm pivoting), a similar contraction holds in the Frobenius norm:

$L \in \mathbb{R}^{n \times k}$2

(Steinerberger, 2024).

Greedy and complete pivoting:

• For spatial kernels $L \in \mathbb{R}^{n \times k}$3, if the diagonals are Lipschitz, complete pivoting ensures the residual after $L \in \mathbb{R}^{n \times k}$4 steps satisfies $L \in \mathbb{R}^{n \times k}$5 for $L \in \mathbb{R}^{n \times k}$6, improving to $L \in \mathbb{R}^{n \times k}$7 for kernels with Lipschitz derivatives (Jeong et al., 16 Sep 2025).
• In the matrix case, the maximal residual satisfies $L \in \mathbb{R}^{n \times k}$8, where $L \in \mathbb{R}^{n \times k}$9 is the next pivot. With smooth kernels and exponentially decaying eigenvalues, convergence is exponential in $k \ll n$0 (Liu et al., 2015).

Spectral bounds:

If $k \ll n$1 approximates $k \ll n$2, the preconditioned system $k \ll n$3 has eigenvalues in $k \ll n$4, $k \ll n$5 under mild assumptions (Jeong et al., 16 Sep 2025).

3. Pivoted Cholesky as a Preconditioner

The incomplete-cholesky factor $k \ll n$6 serves as a preconditioner for iterative solvers, especially conjugate gradient (CG), in solving $k \ll n$7. The preconditioned operator is: $k \ll n$8

Condition number and convergence:

• The condition number satisfies $k \ll n$9.
• For random or complete pivoting, the expected $A^{(0)} = A$0 contracts rapidly toward 1 as $A^{(0)} = A$1 grows (Steinerberger, 2024, Jeong et al., 16 Sep 2025).

Empirical evidence:

• In GPU-accelerated Gaussian process regression, a rank-5 pivoted Cholesky preconditioner can reduce CG iterations by a factor of 2–4, with overall wall-clock solve time halved (Gardner et al., 2018).
• For large kernel matrices, the preconditioner has been shown to yield equivalent accuracy to regularized or full Cholesky at a fraction of the computational cost (Liu et al., 2015, Roos et al., 28 Jul 2025).

4. Practical Implementation Considerations

Memory and cost:

• Factorization to rank $A^{(0)} = A$2 for an $A^{(0)} = A$3 matrix costs $A^{(0)} = A$4 in the basic approach and requires $A^{(0)} = A$5 memory for the factor and auxiliary vectors.
• In kernel methods, column/row updates, as well as Schur complement updates, are required only for the pivot columns and diagonals, facilitating online or batched settings (Gardner et al., 2018).

Preconditioned system application:

• Applying $A^{(0)} = A$6 in CG requires two triangular solves at $A^{(0)} = A$7 per iteration.
• When combined with the Woodbury matrix identity, fast solves and determinant computations are enabled for dense kernel plus diagonal regularization (Gardner et al., 2018).

Stopping criteria and stability:

Pivoted Cholesky naturally stops either at a target rank $A^{(0)} = A$8 or when the next pivot falls below a prescribed tolerance, maintaining a sharp entry-wise approximation guarantee and robustly handling near rank-deficiency without ad hoc regularization (Liu et al., 2015).

Comparison to regularization:

Pivoted Cholesky preconditioning achieves stability without perturbing the original matrix, in contrast to Tikhonov-type regularization, and avoids unnecessary increase in system condition number (Liu et al., 2015).

5. Role in Gaussian Process and Kernel Methods

In Gaussian process inference and kernel ridge regression, the underlying matrix $A^{(0)} = A$9 is often dense and ill-conditioned. The pivoted Cholesky preconditioner enables:

• Scalable exact inference via BBMM (Blackbox Matrix-Matrix Multiplication) and batched CG, reducing training and hyperparameter optimization wall-clock time by up to $i = 1, \dots, k$0 (Gardner et al., 2018).
• Efficient construction of low-rank GP approximations (e.g., Nystrom), with faster A-optimal (PCov) and data-fit aware (WPCov) variants enhancing uncertainty reduction and data fit (Roos et al., 28 Jul 2025).
• Stability in RBF interpolation, Karhunen–Loève expansions, and other kernel sum applications, with factorization cost and accuracy outperforming standard Cholesky or LU decompositions for large $i = 1, \dots, k$1 (Liu et al., 2015).

Notably, in high-dimensional or non-smooth kernel settings (e.g., Matérn-$i = 1, \dots, k$2, Green's functions), empirical performance matches the established algebraic rates without requiring $i = 1, \dots, k$3 regularity (Jeong et al., 16 Sep 2025).

6. Extensions, Limitations, and Empirical Performance

Extensions:

• Cross approximation based pivoted Cholesky methods, as well as hybrid strategies incorporating random and greedy elements, provide flexible trade-offs between theoretical guarantees and practical speed (Steinerberger, 2024, Liu et al., 2015).
• Enhanced variants (e.g., PCov, WPCov) focus directly on experimental design criteria, reducing uncertainty more efficiently in Bayesian nonparametrics (Roos et al., 28 Jul 2025).

Limitations:

• All variants require at minimum SPD entries; further efficiency improvements may need block or hierarchical matrix structures for extremely large $i = 1, \dots, k$4 (Liu et al., 2015).
• Non-orthogonality of $i = 1, \dots, k$5 can be corrected by post-QR steps when extracting eigenvectors.

Empirical performance:

• Large-scale benchmarks indicate order-of-magnitude improvements in solve time for eigen-decomposition and regression tasks, with negligible additional computational burden compared to regularized or full factorizations.
• With high-rank truncation, pivoted Cholesky solvers match machine precision accuracy in direct and iterative settings (Liu et al., 2015, Gardner et al., 2018).

7. Summary of Comparative Properties

Property	Pivoted Cholesky Preconditioner	Diagonal Regularization	Direct Cholesky
Cost (factorization)	$i = 1, \dots, k$6	$i = 1, \dots, k$7	$i = 1, \dots, k$8
Explicit rank truncation	Yes	No	No
Entrywise error control	Yes ($i = 1, \dots, k$9 diag)	No	Yes
Empirical CG speedup	Significant	Marginal	Not applicable
Regularization required	No	Yes	No

Pivoted Cholesky preconditioning achieves scalable, provably controlled approximation of SPD systems in large-scale kernel and regression contexts, yielding significant improvements in computational efficiency while guaranteeing both theoretical and practical convergence properties across a wide class of problems (Steinerberger, 2024, Jeong et al., 16 Sep 2025, Roos et al., 28 Jul 2025, Gardner et al., 2018, Liu et al., 2015).