Variational Free Energy Pivot Selection for Pivoted Cholesky

Abstract: Pivoted Cholesky factorizations construct low-rank approximations of symmetric positive definite matrices by sequentially selecting pivots from the residual diagonal. Classical greedy and randomized rules, such as randomly pivoted Cholesky, target the algebraic trace-norm error of the residual. In many applications, however, the matrix enters a nonlinear matrix functional whose value, not the trace-norm error, determines solution quality, and residual-based rules ignore this structure. We derive a pivot rule that maximizes the exact one-step change of such a functional under Cholesky-consistent rank-1 updates, for a functional combining log-determinant, quadratic, and trace terms. This functional arises as the variational free energy in Gaussian process regression, where the matrix is a kernel matrix. The resulting per-step gain admits a closed-form additive decomposition into complexity, data-fit, and trace contributions, and is used directly as a pivot-selection criterion. We refer to the resulting method as $Δ$-VFE pivoted Cholesky. At each iteration, the criterion is evaluated on a batch of $s$ candidate pivots sampled proportionally to the residual diagonal via incremental Woodbury updates, at a total cost of $\mathcal{O}(snr<sup>2)$ for an $n\times n$ matrix and target rank $r$. This matches the asymptotic complexity of randomly pivoted Cholesky up to the batch factor $s$. Cholesky-consistent rank-1 updates yield monotonically non-decreasing functional values, and the proposed rule maximizes the per-step gain among them. Numerical experiments show improved objective values and predictive accuracy at low to moderate ranks compared to classical and randomly pivoted Cholesky, while preserving trace-norm approximation quality.

• The paper derives an exact per-pivot variational free energy gain that combines complexity, data-fit, and trace terms, enabling objective-aware Cholesky selection with monotonic improvement guarantees.
• The proposed randomized batch method runs in O(snr²) time and O(nr + r²) space, reduces to RPCholesky when s=1, and uses incremental Woodbury updates for efficiency.
• Experiments on Abalone and QM7 show up to 57.79% smaller VFE gaps, improved low-rank prediction accuracy, and trace errors matching or outperforming established baselines.

Motivation and problem statement

Pivoted Cholesky factorizations construct low-rank approximations $K \approx LL^\top$ of symmetric positive definite matrices by greedily or randomly selecting pivots from the residual diagonal. Classical greedy pivoted Cholesky and its randomized variant RPCholesky target the algebraic trace-norm error $\mathrm{tr}(K - \widetilde K)$, a criterion that is well understood and carries provable guarantees. In Gaussian process regression (GPR), however, the quality of an approximation is governed not by the trace error but by the variational free energy (VFE) lower bound on the log marginal likelihood, which combines log-determinant, quadratic data-fit, and trace terms. Residual-diagonal pivot rules explicitly optimize only the third of these components. The paper addresses this mismatch by deriving a pivot rule that maximizes the exact one-step change of the VFE under rank-1 updates that are consistent with the incremental Cholesky factorization.

The per-pivot gain formula

The central theoretical result is a closed-form additive decomposition of the exact change in the VFE induced by a rank-1 update $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$. Writing $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$, $a = \bm u^\top \Sigma_i^{-1}\bm u$, and $b = \bm u^\top \Sigma_i^{-1}\bm y$, the gain is

$$\Delta\mathcal L = -\frac{1}{2}\left(\log(1+a) - \frac{b^2}{1+a} - \frac{1}{\sigma_\varepsilon^2}\|\bm u\|_2^2\right).$$

The three terms correspond exactly to the complexity, data-fit, and trace penalty components of the VFE. Restricting to Cholesky-consistent updates $\bm u_p = \bm t_p/\sqrt{d_p^{(i)}}$, where $\bm t_p$ is the residual column at candidate pivot $p$, yields computable scalars $\mathrm{tr}(K - \widetilde K)$0 and $\mathrm{tr}(K - \widetilde K)$1 expressible purely in terms of quantities already maintained by pivoted Cholesky.

Two structural observations sharpen the interpretation. First, the trace term $\mathrm{tr}(K - \widetilde K)$2 equals the exact one-step decrease in $\mathrm{tr}(K - \widetilde K)$3, namely $\mathrm{tr}(K - \widetilde K)$4, whereas greedy pivoted Cholesky uses the diagonal entry $\mathrm{tr}(K - \widetilde K)$5 only as a tractable proxy for this quantity. Second, using the identity $\mathrm{tr}(K - \widetilde K)$6, the data-fit term is a scaled inner product between the candidate direction and the current predictive residual — a directional quantity, in contrast to the pointwise residual-magnitude criterion of Schreiter et al. This places the method in a three-level hierarchy: diagonal-based rules (greedy, RPCholesky), exact trace-reduction maximizers, and finally the full variational objective.

A notable consequence is that monotonicity holds for any admissible pivot: since $\mathrm{tr}(K - \widetilde K)$7, every Cholesky-consistent update satisfies $\mathrm{tr}(K - \widetilde K)$8, so greedy pivoted Cholesky and RPCholesky also produce monotonically non-decreasing VFE sequences. The proposed rule simply selects, at each step, the admissible pivot maximizing this guaranteed gain.

Algorithmic realization and complexity

Evaluating the gain for all $\mathrm{tr}(K - \widetilde K)$9 candidates would cost $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$0 overall, which is infeasible. The algorithm instead samples a batch $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$1 of $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$2 candidates proportionally to the residual diagonal — the same distribution as RPCholesky, so the method reduces exactly to RPCholesky at $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$3 — and evaluates $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$4 only within the batch. Repeated application of $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$5 is handled via Woodbury identities, with the auxiliary matrix $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$6 updated incrementally by a Schur-complement block formula rather than refactored. For $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$7, total runtime is $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$8 with space complexity $$\widetilde K_{i+1} = \widetilde K_i + \bm u\bm u^\top$$9, matching RPCholesky up to the batch factor $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$0. This assumes kernel columns are formed on demand; precomputing $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$1 raises costs to $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$2 time and $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$3 space.

Relation to prior work

The scalar per-candidate gain was previously derived by Cao et al. in their CholQR algorithm for sparse GP inducing-point optimization, via an augmented QR representation interleaved with local-search swaps and hyperparameter optimization. The contribution here is a re-derivation within the fixed-hyperparameter pivoted Cholesky framework, making the additive decomposition explicit and connecting the rule to the randomized pivoted-Cholesky literature through the shared sampling distribution. Compared to ridge leverage score sampling, the data-fit term $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$4 differs structurally: it is directional, built from the current low-rank approximation rather than the full kernel, depends on $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$5, and evolves incrementally. Relative to Bach–Jordan's Cholesky with side information, the gain here is exact rather than approximated with look-ahead strategies.

Numerical results

Experiments on Abalone ($\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$6, Gaussian kernel) and QM7 ($\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$7, Laplacian kernel), averaged over 15 randomized runs, compare greedy pivoted Cholesky, RPCholesky, and $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$8-VFE pivoted Cholesky up to rank 2048.

• VFE objective: On Abalone at $\Sigma_i = \widetilde K_i + \sigma_\varepsilon^2 I$9, $a = \bm u^\top \Sigma_i^{-1}\bm u$0-VFE reduces the relative gap $a = \bm u^\top \Sigma_i^{-1}\bm u$1 by approximately 57.79% versus greedy pivoted Cholesky and 26.91% versus RPCholesky; the advantage over greedy remains above 70% throughout $a = \bm u^\top \Sigma_i^{-1}\bm u$2. On QM7, gains over both baselines are near 50% at $a = \bm u^\top \Sigma_i^{-1}\bm u$3, peaking at 57.04% over greedy at $a = \bm u^\top \Sigma_i^{-1}\bm u$4. Notably, even at $a = \bm u^\top \Sigma_i^{-1}\bm u$5 no method has converged on QM7, where $a = \bm u^\top \Sigma_i^{-1}\bm u$6-VFE still attains the smallest gap (230.66% vs. 360.89% and 317.43%).
• Trace-norm error: Despite optimizing a different objective, $a = \bm u^\top \Sigma_i^{-1}\bm u$7-VFE matches or improves upon RPCholesky's relative trace error throughout, with gaps reaching 15.88% (Abalone, $a = \bm u^\top \Sigma_i^{-1}\bm u$8) and 59.06% over greedy (QM7, $a = \bm u^\top \Sigma_i^{-1}\bm u$9). This indicates the functional-aware gains do not come at the cost of standard approximation quality.
• Prediction: On Abalone, RMSE improvements over baselines reach roughly 5.5–6% at $b = \bm u^\top \Sigma_i^{-1}\bm y$0; on QM7, MAE reductions reach 20.43% (vs. greedy) and 38.37% (vs. RPCholesky) at $b = \bm u^\top \Sigma_i^{-1}\bm y$1. At very small ranks ($b = \bm u^\top \Sigma_i^{-1}\bm y$2) greedy occasionally wins, so the advantage is specific to low-to-moderate ranks.
• Batch size: On QM7 with $b = \bm u^\top \Sigma_i^{-1}\bm y$3, all batch sizes beat RPCholesky on all metrics; the relative VFE gap decreases monotonically in $b = \bm u^\top \Sigma_i^{-1}\bm y$4 from $b = \bm u^\top \Sigma_i^{-1}\bm y$5 onward but with diminishing returns, and $b = \bm u^\top \Sigma_i^{-1}\bm y$6 already captures most of the achievable gain while yielding strictly lower trace errors than RPCholesky across the displayed range. This supports treating $b = \bm u^\top \Sigma_i^{-1}\bm y$7 as a small constant in practice.

Limitations and open questions

Several caveats are stated in the paper itself. The empirical validation covers two benchmark datasets and two kernels; generalization to other regimes is not established. The batch size $b = \bm u^\top \Sigma_i^{-1}\bm y$8 is treated empirically rather than analyzed theoretically, and no convergence-rate analysis for the resulting pivot sequence is provided. The method operates in the fixed-hyperparameter regime, so the interaction between objective-aware pivoting and hyperparameter optimization — where the data-fit term's dependence on $b = \bm u^\top \Sigma_i^{-1}\bm y$9 could couple with kernel learning — remains unexamined. Finally, the derivation is specific to the VFE functional; extension to other nonlinear matrix functionals is posed as an open direction rather than solved.

Conclusion

This paper derives an exact, closed-form per-step gain for the variational free energy under Cholesky-consistent rank-1 updates and uses it as a greedy pivot-selection criterion implementable at $$\Delta\mathcal L = -\frac{1}{2}\left(\log(1+a) - \frac{b^2}{1+a} - \frac{1}{\sigma_\varepsilon^2}\|\bm u\|_2^2\right).$$0 cost via incremental Woodbury updates. The rule provably produces monotonically non-decreasing VFE values, subsumes RPCholesky as the $$\Delta\mathcal L = -\frac{1}{2}\left(\log(1+a) - \frac{b^2}{1+a} - \frac{1}{\sigma_\varepsilon^2}\|\bm u\|_2^2\right).$$1 special case, and empirically improves both the target objective and predictive accuracy at low-to-moderate ranks without degrading trace-norm approximation quality. The main open questions concern convergence theory, principled batch-size selection, and interaction with hyperparameter learning.