Source-linked AI summary
Krylov-aware stochastic trace estimation
Tyler Chen, Eric Hallman
TL;DR
Estimating tr(f(A)) from matrix-vector products is often approached by treating f(A) as an implicit matrix, which can overlook its Krylov structure. This paper combines low-rank approximation and stochastic trace estimation with Krylov-aware computation, and reports favorable comparisons with existing methods.
Problem
Method
The algorithm combines low-rank approximation and quadratic trace estimation, using block Krylov structure to compute low-rank approximations of f(A) efficiently, including for multiple functions.
Source: §1 · p. 2Results
Numerical experiments show that the methods compare favorably to existing ones.
Source: §1 · p. 2Takeaways & Limitations
Using Krylov structure can reduce required matvecs and make low-rank approximation of matrix functions more efficient, including across multiple functions.
Source: §1 · p. 2Lanczos without reorthogonalization is susceptible to finite-precision effects, and its approximation bounds are often pessimistic.
Source: §2.2 · p. 6
Abstract
from arXiv · showhide
We introduce an algorithm for estimating the trace of a matrix function $f(\mathbf{A})$ using implicit products with a symmetric matrix $\mathbf{A}$. Existing methods for implicit trace estimation of a matrix function tend to treat matrix-vector products with $f(\mathbf{A})$ as a black-box to be computed by a Krylov subspace method. Like other recent algorithms for implicit trace estimation, our approach is based on a combination of deflation and stochastic trace estimation. However, we take a closer look at how products with $f(\mathbf{A})$ are integrated into these approaches which enables several efficiencies not present in previously studied methods. In particular, we describe a Krylov subspace method for computing a low-rank approximation of a matrix function by a computationally efficient projection onto Krylov subspace.
1. Introduction.
The paper studies estimating tr(f(A)) when a symmetric A is available only through matrix-vector products, and argues that using Krylov structure can make this estimation more efficient. For exp(−βA), changing β changes spectral decay and thus which trace-estimation approach is effective.
1. Introduction.
- Estimating tr(f(A)) with only matrix-vector products with symmetric A is the central problem, often split into trace estimation and Krylov approximation of products with f(A) [43].Source: §1 · p. 1
- Existing methods can treat f(A) as an arbitrary implicit matrix, missing its Krylov structure; some methods must be run separately for each implicit matrix [39].Other Krylov-based approaches can approximate traces for multiple functions simultaneously [48].Source: §1 · pp. 1–2
1.1. Motivating example.
- For exp(−βA), large β gives fast spectral decay while small β gives slow decay; low-rank approximation, quadratic trace estimation, or their combination suit different regimes.At intermediate β, the leading eigenvalue contributions can be computed by low-rank approximation and the remainder by quadratic trace estimation.Source: Figure 1.1 · p. 2, §1.1 · pp. 2–3
1. Introduction.
- The proposed algorithm combines low-rank approximation with quadratic trace estimation, exploiting block Krylov structure to reduce matvecs and support multiple functions f.The paper also presents adaptive and memory-limited variants; numerical experiments compare favorably with existing methods.Source: §1 · p. 2
1.2. Notation.
- The notation defines f(A) through the eigendecomposition A = UΛU^T, applying f to each diagonal entry of Λ.It also establishes conventions for vector and matrix symbols and norms.Source: §1.2 · p. 3
2. Background.
The background develops stochastic trace estimation with deflation, block-Krylov approximations to matrix functions, and an algorithm that combines these ideas. Its accuracy and practical use depend on Krylov approximation quality and finite-precision behavior.
2.1. Implicit trace estimation.
- Projecting random vectors onto the complement of Q’s column span yields vectors uniformly distributed on that complement’s unit sphere and supports an unbiased remainder estimator.For Gaussian vectors, normalizing by (d − b)/y_i^T y_i can improve performance when d − b is finite, particularly for a flat remainder spectrum [12, §2.3].Source: §2.1 · pp. 3–4
- If Q approximates B’s dominant eigenspace, the remainder can have smaller Frobenius norm, reducing the random-vector count needed for quadratic trace estimation.A sketch-based choice uses an orthonormal basis for BΩ; Algorithm 2.1 uses 2b + m matvecs with B.Source: §2.1 · p. 4
2.2. Krylov subspace methods for matrix function approximation.
- Block Lanczos builds an orthonormal basis for a block Krylov subspace through a symmetric block-tridiagonal recurrence.Under the stated full-dimension assumption, the matrix approximation is exact for polynomials of degree at most q − 1, and the quadratic-form approximation for degree at most 2q − 1.Source: §2.2 · pp. 4–5
- Bounds for Krylov matrix-function approximations can guide polynomial degree, but bounds are often pessimistic and finite-precision effects reduce convergence rates without reorthogonalization.Developing practical spectrum-dependent stopping bounds remains an active research area.Source: §2.2 · p. 6
2.3. A simple algorithm for the trace of a matrix function.
- Algorithm 2.3 uses Krylov approximations to construct a projection for estimating tr(f(A)); with sufficiently accurate f(A)Ω, its projection behaves almost identically to the exact-product version.The source notes that increasing q makes its projection converge to the one generated by Hutch++ with exact products, though that may be an inefficient use of computation.Source: §2.3 · p. 7
3. Krylov-aware function approximation.
The section presents a Krylov-aware low-rank approximation for stochastic trace estimation that exploits block-Krylov structure to reduce matrix-vector products. It discusses the approximation’s polynomial exactness, parameter trade-offs, and how the approach differs from prior work.
3. Krylov-aware function approximation.
- Using the block-Krylov relation, the method approximates the projection of f(A) onto Kq+1(A, Ω) with products involving matrices of size n × b.The approximation is exact for polynomials of degree at most 2n − 1 when Kq+n has dimension (q + n)b.Source: §3 · pp. 7–8
3.1. Main algorithm.
- Algorithm 3.1 combines low-rank approximation and quadratic trace estimation, exploiting block-Krylov structure to reduce matvecs and support multiple functions f.The algorithm’s stochastic trace-estimation steps can be run in parallel, including alongside the final block-Lanczos iterations.Source: §3.1 · pp. 8–9
- Algorithm 3.1 often performs significantly better than Algorithm 2.3 in the authors’ experiments.The theorem discussion also states that when dΔ is small relative to max(1, V), Algorithm 3.1 cannot perform significantly worse than Algorithm 2.3.Source: §3.1 · p. 10
3.2. Choice of parameters.
- With deflation-space size bounded, the authors favor small block size b and large Krylov depth q because the main expected savings are in matvecs, not matrix loads.The minimum matrix-load count is q + n, and the block-Krylov basis also incurs storage and reorthogonalization costs.Source: §3.2 · pp. 10–11
3.3. Relation to past work.
- Unlike the n = 1, m = 0 special case in, the approach allows any n, which matters when dominant eigenvalues of f(A) do not align with those of A.The paper also adapts the correlated-sample idea of to deflation only, preserving an unbiased trace estimator.Source: §3.3 · pp. 11–12
- The paper leaves a priori bounds for its Krylov-aware low-rank approximation to arbitrary f outside its scope.The authors identify understanding the theoretical behavior of these algorithms as further work.Source: §3.3 · p. 12
4. Variants.
The variants address parameter selection and memory constraints in the main algorithm. The adaptive variant chooses effort for variance reduction and trace sampling, while the low-memory variant uses restarting to limit stored Krylov basis vectors.
4.1. Adaptive variant.
- The adaptive variant addresses uncertainty about balancing variance reduction against quadratic trace estimation by selecting Lanczos iterations and estimator samples.It assumes a fixed block size and estimates the remainder norm to guide the sample count.Source: §4.1 · p. 12
- Gaussian sketches provide a probabilistic upper estimate of the remainder Frobenius norm, and the algorithm reuses samples to estimate that norm and the trace while choosing when to stop sampling.The norm estimate uses a chi-squared calibration tied to the user-specified failure tolerance.Source: §4.1 · p. 13
- The adaptive variant uses block-Lanczos estimates of computational cost to choose the variance-reduction iterations, stopping when the estimated objective indicates a local minimum.Matvec count is used as a cost proxy; the objective estimate is delayed by n iterations, whose work is already required by the algorithm.Source: §4.1 · pp. 13–14
4.2. Low-memory variant.
- The low-memory variant restarts block-Lanczos to approximate the dominant eigenspace of f(A) without storing the full Krylov basis.It focuses on implicitly restarted block-Lanczos, while noting that other techniques may suit specific situations.Source: §4.2 · pp. 14–15
- Restarting proceeds in phases so the number of restarts and filter polynomials can adapt to information gained, and the process can stop once the approximation is sufficiently desirable.A practical algorithm combining restarting with adaptive parameter selection is not specified because implementation depends on the problem and computing system.Source: §4.2 · pp. 15–16
5. Numerical experiments.
Across the experiments, the methods combine low-rank deflation with trace estimation, with performance depending on the spectrum, block size, and workload. The adaptive method reduced matvecs on two tests, while gains were limited on a well-conditioned log-determinant problem.
5.1. Quantum spin systems.
- For the spin chain, combining low-rank approximation and quadratic trace estimation gave more uniform accuracy across β than either approach alone at about the same cost.The combined method was less accurate at extreme β values but better at intermediate values; the reported errors are 90th percentiles from 100 trials.Source: §5.1 · pp. 16–17, Figure 5.1 · p. 17
- Restarting improved low-rank approximation quality while limiting storage costs, and sufficient restarts produced results comparable to or better than the same-cost combined approach.The figure compares zero, two, and four restarts.Source: §5.1 · pp. 17–18, Figure 5.1 · p. 17
5.2. Parameter selection for the inverse.
- For algebraically decaying spectra, restarting helped only mildly, whereas for geometric decay it enabled higher-quality projection approximations.Larger q, b, and r generally improved approximation quality; the figure omits cases with qb > 1024.Source: §5.2 · p. 18, Figure 5.2 · p. 19
- The black-box Krylov projection used by Algorithm 2.3 stopped improving once Lanczos approximations were accurate enough and was worse than the proposed projection at fixed q and b.The passage attributes the gap to Algorithm 2.3’s projection space being a strict subspace of Q̄_{q+1}.Source: §5.2 · p. 18
5.3. Adaptive algorithm.
- The adaptive experiments compared Algorithm 4.1 with A-Hutch++, whose remainder-estimation procedures were identical, so differences centered on deflation vectors and deflation cost.The tests covered synthetic data and two SuiteSparse matrices.Source: §5.3 · p. 18
5.3.1. Estrada index.
- On the Roget’s graph, Algorithm 4.1 used less than half as many matvecs as A-Hutch++ at p = 2, with its advantage reaching nearly 7× at p = 7.It used about 10 times as many deflation vectors but fewer matvecs for deflation; matvec counts do not capture orthogonalization and storage costs.Source: §5.3.1 · p. 19, Table 5.2 · p. 20
5.3.2. Log Determinant.
- On the well-conditioned thermomech TC matrix, Algorithm 4.1 and A-Hutch++ performed similarly because variance reduction was not particularly useful; the proposed method had lower overhead.The matrix had κ(A) ≈ 70.Source: §5.3.2 · p. 21, Table 5.4 · p. 21
5.3.3. Nuclear norm.
- Algorithm 4.1’s efficiency advantage on the synthetic test grew from about 2× at p = 2 to 7× at p = 7.It used about 20 times as many deflation vectors as A-Hutch++ while using fewer for estimating the remainder.Source: §5.3.3 · pp. 21–22, Table 5.5 · p. 22
6. Proofs.
The proofs establish polynomial exactness properties of the Krylov-based approximations and show that the trace estimator is unbiased.
6. Proofs.
- Approximation (2.6) is exact when f is a polynomial of degree at most q − 1.The proof uses linearity.Source: §6 · p. 22
- Approximation (2.7) is exact when f is a polynomial of degree at most 2q − 1.The proof partitions k ≤ 2q − 1 as k = 2j + ℓ, with j ≤ q − 1 and ℓ ∈ {0, 1}.Source: §6 · p. 23
- The approximation error is bounded using the best polynomial approximation to f over [λmin, λmax].The displayed bound is at most 2∥Z∥² max over that interval of |f(x) − p(x)|; the proof then optimizes over p.Source: §6 · p. 23
- The trace estimator is unbiased: E[est] = tr(f(A)).The proof uses normalized Gaussian vectors projected onto the complement of the deflation subspace.Source: §6 · p. 24
7. Conclusions and future work.
The conclusions report reduced matrix-vector product counts from block Krylov structure, while noting that computational-time benefits remain uncertain. They identify implementation costs and approximation-degree selection as areas for further study.
7. Conclusions and future work.
- Exploiting block Krylov structure can significantly reduce matrix-vector products, but its effect on runtime is unclear.Runtime depends on the cost of maintaining the orthonormal basis and how efficiently matrix-vector products with A run in parallel.Source: §7 · p. 25
- Selecting a smaller approximation degree after effective deflation may be possible, but the paper reports no formal results on this question.The cited paper [28] proposes using multiple degrees to hedge against overestimating n.Source: §7 · p. 25