Stochastic Lanczos Quadrature for Computational Uncertainty in Linear Algebra

1 Introduction

• A truncated Lanczos approximation assigns zero to the orthogonal complement, which can hold substantial inverse mass when the matrix is ill-conditioned.The complement contains the smallest eigenvalues of A, which become the largest under inversion.Source: §1 · p. 2
• The method models the missing uncertainty with an exact rank-one coupling at the truncation boundary and an isotropic correction on the remaining complement bulk.The coupling comes from the Lanczos residual; P-SLQ estimates the bulk magnitude.Source: §1 · p. 2
• A single Lanczos run yields both a matrix-free inverse approximation and a Gaussian belief over the exact A−1b that contracts as the resolved subspace grows.The covariance separates resolved variance from variance remaining in the complement.Source: §1 · p. 3
• P-SLQ estimates the bulk magnitude with re-projected Lanczos probes that reuse one Krylov basis, at O(m + Nl) matrix-vector products rather than O(Nm).Its estimator is consistent, with bias decaying as O(ρ−2l) and variance O(N−1).Source: §1 · p. 3
• On FashionMNIST/CNN and CIFAR-10/ResNet-9, the corrected covariance achieves near-zero per-point KL to the exact posterior and matches exact NLL and ECE within standard errors; P-SLQ also beats some trace-estimation baselines at matched budgets.The trace-estimation comparison spans three conditioning regimes.Source: §1 · pp. 3–4

2.1 Krylov Methods for Matrix Functions

• Krylov error analyses study vector approximation, while this work uses the classical residual coupling to identify and model the inverse’s true complement block.The truncation assigns zero on the complement, which can carry substantial inverse spectral mass for ill-conditioned matrices.Source: §2.1 · p. 4

2.2 The Laplace Approximation in Deep Learning

• Common Hessian simplifications discard curvature and assign zero posterior variance in neglected directions; the method measures and corrects this uncertainty without structural assumptions on A.The matrix-free Lanczos approach avoids structural assumptions but still inherits the truncation problem.Source: §2.2 · p. 5

2.3 Stochastic Trace Estimation

• P-SLQ re-projects every Lanczos step into the fixed complement, addressing why naive SLQ can leave that subspace and corrupt the estimated spectral density.The complement projector is built before probes are drawn, so the compressed operator is fixed for them.Source: §2.3 · p. 5

2.4 Probabilistic Numerics and Computation-Aware Inference

• Unlike computation-aware linear-solve methods, this approach models the full complement and its exact rank-one boundary coupling for a matrix function.ULISSE instead debiases individual solves with Russian roulette, trading unbiasedness for variance that can grow as truncation is pushed earlier.Source: §2.4 · pp. 5–6

2.5 Low-Rank-Plus-Shift Approximations

• The isotropic bulk correction resembles scalar-shift approximations that preserve discarded trace, but it targets the coupling-free Krylov complement and estimates its magnitude matrix-free with P-SLQ.Related low-rank-plus-shift constructions include FITC and deflation preconditioning.Source: §2.5 · p. 6

3.1 Problem Statement

• For an SPD matrix accessed only through matrix-vector products, the objective is to quantify and restore the uncertainty discarded by a rank-m inverse approximation.The missing covariance has an exact rank-one boundary coupling from the Lanczos residual and an isotropic variance on the uncoupled bulk.Source: §3.1 · pp. 6–7

3.2 The Lanczos Algorithm

• Lanczos builds an orthonormal basis for the Krylov subspace from a reference vector using a three-term recurrence, producing the tridiagonal restriction of A.The basis is built once and reused for covariance approximation, sampling, and the later probe computations.Source: §3.2 · p. 7

3.3 Lanczos-Based Approximation of the Inverse

• Projecting A−1 onto the Krylov subspace yields the standard estimate Qm Tm−1 b from m matrix-vector products, but assigns zero to the complement.The paper quantifies that discarded variance and restores it with a matrix-free covariance factor.Source: §3.3 · p. 8

4 Computational Uncertainty on the Krylov Complement

• The resolved and unresolved subspaces couple only through a rank-one block with spectral norm βm, which vanishes exactly when Lanczos converges.The coupling links the last resolved direction to the first unresolved direction; setting βm = 0 treats the resolved subspace as invariant.Source: §4 · p. 9

4.1 Truncation Error and Complement Inverse

• The true inverse block on the complement is the inverse of rank-one-regularized S⊥, not the compressed operator’s inverse; omitting the coupling under-covers the boundary direction by a deficit governed by γ.Setting βm = 0 gives γ = 0 and recovers the compressed target.Source: §4.1 · p. 10

The exact block inverse and the coupling vectors

• The exact block inverse under rank-one coupling has four blocks governed by two coupling vectors, one on the resolved side and one on the complement side.The resolved vector uses Qm, Tm, and βm; the complement vector propagates the boundary direction through S⊥^-1.Source: The exact block inverse and the… · p. 10

The corrected Inverse Approximation: exact coupling, isotropic bulk

• The corrected inverse preserves the rank-one boundary coupling exactly and approximates only the coupling-free complement bulk isotropically.The bulk approximation preserves the complement trace; its isotropy is the remaining approximation.Source: The corrected Inverse Approximation:… · pp. 11–12

Sampling from the corrected inverse as covariance

• A structured square root of the corrected covariance yields Gaussian samples, with a two-dimensional coupling draw encoding the resolved–complement correlation and bulk noise restricted to the remaining complement.The bulk sampler removes both the resolved component and the coupling direction.Source: Sampling from the corrected inverse as… · p. 13

4.2 Derivation of the Complement Variance ω(m)

• The isotropic bulk variance is the mean inverse eigenvalue of S⊥ over the coupling-free bulk, excluding the coupling direction whose variance is retained exactly.When βm = 0, the coupling vectors and direction vanish and the formulas reduce to the decoupled case.Source: §4.2 · pp. 14–15

4.3 Justification of the Isotropic Complement Approximation

• The isotropic complement approximation is exact for a uniform complement spectrum, and its error is controlled by spectral spread; the complement condition number decreases with m and tends to 1 in deep networks as dominant outlier curvature is resolved.The complement condition number equals 1 exactly when all tail eigenvalues are equal.Source: §4.3 · p. 15
• In exact arithmetic, κ⊥(m) decreases from κ toward 1 as Krylov directions are resolved; finite-precision loss of orthogonality can cause small fluctuations.In deep networks, the decrease is rapid once dominant outlier eigenvalues leave the complement.Source: §4.3 · p. 17
• The corrected covariance’s bulk block D has interlacing eigenvalues whose exact mean is ω(m), while its cross-term c is separately bounded.The cross-term bound combines with the isotropic block error to bound the total covariance error.Source: §4.3 · p. 19
• Whenever β_m ≠ 0, the decoupled model underestimates total variance, whereas the coupled model assigns exactly the right amount.The two strictly positive variance terms are resolved-block inflation and complement trace deficit.Source: §4.3 · p. 21

4.4 Computational Uncertainty from the Exact Decomposition

• Σ⊥_m is the exact Gaussian posterior covariance for A^-1v after matrix-free Lanczos observations, representing epistemic uncertainty from incomplete information.The observation basis is fixed by the reference vector and reused for every probe.Source: §4.4 · pp. 21–22
• The boundary residual adds rank-one uncertainty inside the resolved block, while the covariance null space is A range(Q_m), not range(Q_m), and the gap vanishes at convergence.The residual mediates both effects through β_m.Source: §4.4 · p. 23
• The covariance is positive semidefinite with rank d−m, null space A range(Q_m), and positive values on its orthogonal complement.Source: §4.4 · p. 23
• As in Cockayne et al. (2019a, Prop. 3), m observations determine m directions and leave d−m uncertain, but here uncertainty is assigned by the operator-determined Schur complement.Source: §4.4 · pp. 23–24
• The Gaussian belief has the m-step Lanczos estimate as its mean and the decomposition remainder as covariance, and concentrates on the exact value as m approaches d.Source: §4.4 · p. 24
• The decomposition supports posterior sampling with a matrix-free inverse approximation and trace estimation with computational uncertainty in A^-1v estimated by P-SLQ.Both interpretations use boundary vectors q_parallel and q_perp and an isotropic bulk term omega(m), recovered from one Lanczos run.Source: §4.4 · p. 24

5.1 Complement Invariance via Re-projection

• Per-step reprojection keeps Lanczos vectors in the Krylov complement, so the iteration represents the operator restricted to that subspace.Reorthogonalisation alone enforces orthogonality within the quadrature basis; complement invariance is the role of reprojection.Source: §5.1 · pp. 24–25

5.2 Consistency and Convergence of P-SLQ

• P-SLQ is consistent, with quadrature bias decaying exponentially in depth and estimator variance of order O(N^-1).The stated variance rate holds for every fixed quadrature depth; as depth grows, the variance bound recovers the pure Monte Carlo bound.Source: §5.2 · pp. 25–26

5.3 Estimating the Coupling-Corrected Magnitude

• A deterministic Sherman-Morrison correction uses boundary data to transfer P-SLQ’s consistency guarantees to the coupling-corrected trace.The additional cost is one boundary probe and one conjugate-gradient solve.Source: §5.2 · p. 28, §5.3 · p. 30

6.1 P-SLQ estimator tests

• Section 6.1 tests P-SLQ consistency, the role of reprojection, and trace-estimation accuracy against other methods at matched matrix-vector budgets.The comparisons include unbiased solver ULISSE and state-of-the-art trace-estimation techniques.Source: §6 · p. 32

6.2 Laplace approximation

• Section 6.2 applies the full pipeline to the Laplace approximation in Bayesian deep learning.Source: §6 · p. 32

• The validation is organized in two stages: estimator tests first, followed by the Laplace-approximation application.Source: §6 · p. 32

Consistency

• Bias decays geometrically with quadrature depth and relative variance as N^-1 across all three spectra.Bias decays fastest on the well-conditioned bulk-plus-outliers regime and slowest on the exponential regime, consistent with complement conditioning.Source: §6.1 · pp. 32–33, Figure 1 · p. 33

The role of reprojection

• Without per-step reprojection, relative bias grows with quadrature depth on the power-law spectrum, despite identical full reorthogonalisation.On exponential and bulk-plus-outlier spectra, the two variants give the same result up to noise.Source: §6.1 · p. 33, Figure 2 · p. 34

Trace estimation

• P-SLQ achieves the lowest relative RMSE across all tested spectral regimes at matched matrix-vector budgets.The comparison includes ULISSE and other trace estimators; the reported trade-off is a small, exponentially decaying bias for O(N^-1) variance independent of m.Source: §6.1 · p. 34, Figure 3 · p. 35

6.2 Application for Laplace Approximation in Bayesian deep learning

• The application compares corrected and exact predictive covariances for Laplace approximations on FashionMNIST and CIFAR-10.The study evaluates a convolutional network on FashionMNIST and a ResNet-9-style network without BatchNorm on CIFAR-10.Source: §6.2 · p. 35, §6.2.1 · p. 35

6.2.1 Setup

• The FashionMNIST model has 11,878 parameters and permits forming and inverting the full Laplace posterior; the CIFAR-10 model has 672,890 parameters and does not.Results are reported as mean ± standard error over five seeds, with only stochastic components varying across seeds.Source: §6.2.1 · p. 35
• The CIFAR-10 generalized Gauss-Newton matrix is assembled from 200 training points and rescaled, while its prior precision is fitted on the unscaled operator.The KFAC surrogate is fitted on the same 200 points without rescaling, making its curvature small relative to the evaluated operator; the text characterizes this as conservative for low-rank KFAC.Source: §6.2.1 · pp. 35–36

6.2.2 Metrics

• Predictive fidelity is assessed with per-point Gaussian KL divergence and log-trace, alongside test negative log-likelihood, 15-bin expected calibration error, and Brier score.Singular approximate covariance blocks yield infinite KL and are excluded, with their excluded fraction reported.Source: §6.2.2 · p. 36, Predictive performance · p. 36

6.2.3 Subspace baselines

• Rank-s projection baselines assign zero variance to the remaining directions and under-cover predictive uncertainty; the CIFAR-10 full-GGN optimum is omitted because it exceeds the compute budget.On CIFAR-10, the comparison is therefore against feasible methods, not the optimum reference ceiling.Source: §6.2.3 · pp. 36–37

6.2.4 Predictive performance

• On FashionMNIST at s = 1000, Σm nearly matches the exact full Laplace on NLL (0.2499 versus 0.2498) and Brier score (0.1313 versus 0.1312), while KL separates posterior fidelity more strongly.On CIFAR-10, NLL, ECE, and Brier differences are within one or two standard errors for almost every method, while KL values differ by four orders of magnitude.Source: §6.2.4 · pp. 37, 39, 40, Table 1 · p. 38, Table 2 · p. 39

6.2.5 Fidelity: KL divergence and log-trace

• At s = 1000, FashionMNIST Σm reaches KL = 0.0160 versus 0.9755 for full-rank KFAC; on CIFAR-10 at s = 200, it reaches 0.0087 versus 0.3907 for KFAC.The full-rank comparison is reported over five seeds; KL values below 10^-4 are at the resolution of the CG reference.Source: §6.2.6 · p. 41, Table 3 · p. 41

6.2.6 Key Observations

• On CIFAR-10, the low-rank-plus-shift baseline is more accurate at s = 20, 50, and 100, but Σm overtakes it at s = 200; on FashionMNIST, Σm leads at every s.The CIFAR-10 GGN is subset-built and far below full rank, making the inverse low-rank plus a multiple of the identity, the form assumed by the shift baseline.Source: §6.2.6 · p. 41, §6.2.5 · p. 41
• Σm is full rank, so the rank-s projection trace ceiling does not bind it; its trace sits at the ceiling on CIFAR-10 and slightly above it on FashionMNIST.Every projection baseline sits below the ceiling on both datasets.Source: §6.2.5 · p. 41, §6.2.6 · p. 42
• The subset-Diagonal baseline has infinite KL because it selects blocks with numerically zero variance, so it is excluded from the tables.The authors report that Faller and Martin (2025) observed the same behavior and attributed it to many parameters having gradients near zero.Source: §6.2.6 · p. 42

7 Conclusion

• The corrected covariance retains the Lanczos-residual boundary coupling exactly, matches the total discarded variance, and approximates only the coupling-free complement bulk isotropically.Its error is controlled by the computable diagnostic κ⊥.Source: §7 · p. 42
• A single Lanczos run yields both a matrix-free approximation of A^-1 and a Gaussian belief over A^-1v that contracts toward the exact answer as computation proceeds.P-SLQ supplies the complement trace using a reused Krylov basis, with exponentially decaying bias and variance of order O(N^-1).Source: §7 · p. 42
• At equal matrix-vector budgets, P-SLQ had the lowest relative RMSE on all four tested matrices, with bias decaying exponentially in quadrature depth and variance of order O(N^-1).The tested set included a kernel matrix not constructed from a chosen spectrum.Source: §7 · p. 42
• When κ⊥ remains large, one isotropic variance cannot represent the complement and the approximation worsens.The authors identify this as the case when the complement spectrum decays quickly.Source: §7 · p. 43