Implicit Trace Estimation Methods

• Implicit trace estimation is the process of approximating a matrix or operator’s trace using oracle access methods like quadratic forms and matrix-vector products.
• Randomized probing, structure-aware sampling, and graph-based techniques are employed to control variance and enhance accuracy in trace approximation.
• Modern methods such as Hutch++ and XTrace leverage low-rank deflation and exchangeable estimators to accelerate convergence and optimize query complexity.

Implicit trace estimation is the problem of approximating $\operatorname{tr}(A)$, or more generally $\operatorname{tr}(f(A))$, when the target matrix or operator is not available entrywise and can only be accessed through oracle operations such as quadratic forms $x^\dagger A x$, matrix-vector products $w \mapsto Aw$, or operator-function products $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$ (Fitzsimons et al., 2016, Epperly et al., 2023, Zvonek et al., 2023). The standard mechanism is randomized probing: choose vectors or functions whose second-moment structure implies an identity such as $\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$, or $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$ in normalized quadratic-form models, and then reduce error by averaging, deflation, or structure-aware sketching (Fitzsimons et al., 2016, Roosta-Khorasani et al., 2013). The subject spans worst-case variance analysis, high-probability sample complexity, low-rank variance reduction, Krylov and Nyström methods, graph-structured probing, dynamic updates, and operator-theoretic generalizations.

1. Oracle models and problem formulations

The finite-dimensional formulations split into two closely related access models. In the quadratic-form model, an algorithm queries a vector $x$ and receives $x^T A x$ or $x^\dagger A x$; this is the model used in the query-complexity analysis of trace estimation and in the mutually unbiased bases construction (Wimmer et al., 2014, Fitzsimons et al., 2016). In the matrix-vector model, the algorithm applies $\operatorname{tr}(f(A))$0 to chosen probes $\operatorname{tr}(f(A))$1 and then forms $\operatorname{tr}(f(A))$2; this is the setting of Girard–Hutchinson, Hutch++, XTrace, and related randomized sketching methods (Epperly et al., 2023, Jiang et al., 2021).

The same formulation extends naturally to matrix functions. When $\operatorname{tr}(f(A))$3 is symmetric and $\operatorname{tr}(f(A))$4 is applied spectrally, the target becomes

$\operatorname{tr}(f(A))$5

but $\operatorname{tr}(f(A))$6 is still unavailable explicitly, so the algorithm combines randomized trace identities with Krylov approximations to $\operatorname{tr}(f(A))$7 (Chen et al., 2022, Frommer et al., 2023). A further extension replaces matrices by trace-class integral operators $\operatorname{tr}(f(A))$8, whose trace is

$\operatorname{tr}(f(A))$9

while access is only through operator-function products (Zvonek et al., 2023).

Dynamic variants replace a single implicit object by a sequence $x^\dagger A x$0 or, more generally, states $x^\dagger A x$1 that change slowly over time. The trace-estimation task then becomes one of maintaining accurate estimates while reusing information from previous time steps (Dharangutte et al., 2021, Gokhale et al., 22 Jun 2026). This suggests that implicit trace estimation is best viewed as a family of oracle-identification problems whose common invariant is the extraction of trace information without explicit diagonal access.

2. Classical Monte Carlo estimators and fundamental complexity

The classical estimators are all quadratic-form averages, but they differ in probe distribution, variance, and high-probability guarantees. In the matrix-vector model, the Girard–Hutchinson estimator uses isotropic random vectors $x^\dagger A x$2 with $x^\dagger A x$3, so that

$x^\dagger A x$4

is unbiased (Epperly et al., 2023). In the real SPSD setting, sharper sample-size guarantees are known: with $x^\dagger A x$5, Hutchinson satisfies the probabilistic relative-error guarantee if

$x^\dagger A x$6

while the Gaussian estimator satisfies it if

$x^\dagger A x$7

(Roosta-Khorasani et al., 2013).

The principal estimator families and their representative one-sample properties are summarized below.

Estimator family	Probe construction	Representative property
Fixed basis	$x^\dagger A x$8 uniform in a fixed orthonormal basis; estimate $x^\dagger A x$9	$w \mapsto Aw$0
Hutchinson	i.i.d. Rademacher entries	$w \mapsto Aw$1
Gaussian	i.i.d. $w \mapsto Aw$2 entries	$w \mapsto Aw$3
Unit-vector	columns of $w \mapsto Aw$4, with or without replacement	without replacement has smaller variance

Classical theory also gives matrix-sensitive sufficient conditions. For Hutchinson, the bound depends on

$w \mapsto Aw$5

so diagonal matrices are exact in one sample. For Gaussian probing, performance depends on

$w \mapsto Aw$6

which is small when the spectrum is not too skewed. For unit-vector estimators, the relevant quantity is

$w \mapsto Aw$7

so nearly constant diagonals are favorable (Roosta-Khorasani et al., 2013).

Query-complexity lower bounds show that these Monte Carlo rates are not artifacts of specific probe distributions. In the quadratic-form oracle model, any estimator requires $w \mapsto Aw$8 queries to have variance at most $w \mapsto Aw$9, and any estimator requires

$(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$0

queries to achieve a $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$1-multiplicative approximation with probability at least $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$2 (Wimmer et al., 2014). The same work gives an exact optimum within the class of linear nonadaptive unbiased estimators: sample $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$3 random orthogonal unit vectors and output

$(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$4

This establishes a sharp distinction between asymptotic optimality of sample complexity and variance optimality of specific probe designs.

3. Probe geometry, orthogonality, and structure-aware sampling

A major line of work improves trace estimation not by changing the estimator form, but by choosing probes with stronger geometric properties. The most explicit example is the estimator based on mutually unbiased bases (MUBs). If $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$5 is a collection of mutually unbiased bases, with

$(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$6

then sampling a basis uniformly, sampling a vector uniformly within it, and estimating $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$7 preserves unbiasedness while improving both worst-case single-shot variance and randomness cost (Fitzsimons et al., 2016).

For prime or prime-power dimensions, $(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$8 is available, and the exact MUB variance is

$(Fu)(x)=\int_\Omega f(x,y)u(y)\,dy$9

In the positive semidefinite case this yields

$\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$0

while Hutchinson has

$\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$1

and Gaussian probing has

$\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$2

(Fitzsimons et al., 2016). The same estimator requires only $\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$3 random bits per probe, versus $\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$4 random bits for Hutchinson and effectively infinite precision for exact Gaussian sampling.

A complementary structure-aware route is stochastic probing based on graph colorings. For sparse symmetric $\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$5 and matrix functions $\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$6 with distance-decaying entries,

$\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$7

one first partitions vertices by a distance-$\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$8 coloring and then replaces each deterministic probing vector on a color class by a random Rademacher vector supported on that class (Frommer et al., 2023). The estimator remains unbiased, and for Rademacher probes its variance decomposes as

$\mathbb{E}[w^*(Aw)]=\operatorname{tr}(A)$9

Under banded, lattice, or constant-sign assumptions, the analysis shows a scaling improvement from linear in $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$0 for deterministic probing to square-root in $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$1 for stochastic probing (Frommer et al., 2023). Under a constant-sign condition on the off-diagonal entries within color classes, using one stochastic vector per color is always better than deterministic probing with the same coloring.

These results collectively suggest that probe geometry matters at three levels: isotropy controls unbiasedness, orthogonality controls cross-covariance, and problem-specific locality or unbiased-basis structure can reduce variance far below that of generic iid probes.

4. Low-rank deflation, Hutch++, and exchangeable estimators

Modern variance-reduction methods are built around deflation. Hutch++ decomposes the trace into a low-rank part, approximated by a randomized sketch, and a residual part, estimated stochastically. For symmetric $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$2, the estimator has the form

$\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$3

and for symmetric PSD matrices it achieves relative error with only $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$4 matrix-vector products, improving on Hutchinson’s $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$5 behavior (Persson et al., 2021). Adaptive Hutch++ chooses the split between approximation and residual estimation from the observed matrix structure, while Nyström++ replaces randomized SVD by a one-pass Nyström approximation for PSD matrices and retains Hutch++-style asymptotics (Persson et al., 2021).

The nonadaptive complexity of this paradigm is now understood sharply. For PSD matrices, NA-Hutch++ can be implemented with

$\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$6

matrix-vector products, and any nonadaptive algorithm must use at least

$\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$7

queries (Jiang et al., 2021). This identifies sketch-based, mergeable, and parallelizable trace estimation as essentially optimal in the high-probability nonadaptive regime.

XTrace and XNysTrace refine the same deflation principle by enforcing exchangeability. Rather than splitting samples into separate “approximation” and “residual” roles, they form leave-one-out estimators in which every sample contributes symmetrically to both tasks (Epperly et al., 2023). The resulting estimators are unbiased and have variance decaying as $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$8, in contrast to Hutchinson’s $\mathbb{E}[x^\dagger A x]=\operatorname{tr}(A)/n$9. For exponentially decaying eigenvalues $x$0,

$x$1

and

$x$2

versus $x$3 for Hutch++ (Epperly et al., 2023). On a partition-function computation for a quantum spin model, with $x$4 matvecs the variance-reduced methods achieved errors about five orders of magnitude smaller than Hutchinson, and at $x$5, XTrace was roughly $x$6 more accurate than Hutch++, while XNysTrace was about $x$7 more accurate (Epperly et al., 2023).

Subsequent work on XTrace shows that averaging over right-orthogonal rotations offers only slight practical benefits, whereas replacing the single power-iteration sketch by the full Krylov space $x$8 can lead to significant improvements depending on the spectrum (Hallman, 2 Dec 2025). A plausible implication is that the main unresolved issue in deflation-based estimators is no longer whether low-rank information helps, but how aggressively that information should be recycled without destroying unbiasedness or computational balance.

5. Matrix functions, Krylov structure, and continuous operators

For matrix functions, the dominant issue is that products with $x$9 are themselves approximate. The Krylov-aware approach addresses this by refusing to treat $x^T A x$0-matvecs as a black box. If block-Lanczos on $x^T A x$1 generates $x^T A x$2, then the estimator

$x^T A x$3

with

$x^T A x$4

reuses the same Krylov structure for both deflation and residual correction (Chen et al., 2022). The underlying approximation is exact for polynomials of degree $x^T A x$5 in $x^T A x$6 and of degree $x^T A x$7 in $x^T A x$8, and the method offers fewer matvecs, better low-rank approximations of $x^T A x$9, reuse for multiple functions $x^\dagger A x$0, and compatibility with parallel blocked matvecs and restart schemes (Chen et al., 2022).

The Block Krylov perspective also clarifies upper and lower limits. Approximating $x^\dagger A x$1 by a degree-$x^\dagger A x$2 polynomial $x^\dagger A x$3 identifies Krylov depth with polynomial degree. For $x^\dagger A x$4 and $x^\dagger A x$5 on $x^\dagger A x$6, the required degree is

$x^\dagger A x$7

and hence the same order of Krylov steps suffices (Yu, 28 Jun 2025). However, lower bounds for Wishart inputs show that for $x^\dagger A x$8 with $x^\dagger A x$9, constant-factor estimation requires $\operatorname{tr}(f(A))$00 matrix-vector queries (Yu, 28 Jun 2025). This demonstrates that the efficiency of Krylov-based trace estimation is constrained jointly by approximation theory and information-theoretic query limits.

In the infinite-dimensional setting, ContHutch++ extends Hutch++ from matrices to trace-class integral operators. Continuous Hutchinson replaces Gaussian vectors by Gaussian-process probes $\operatorname{tr}(f(A))$01, yielding

$\operatorname{tr}(f(A))$02

which becomes asymptotically unbiased for $\operatorname{tr}(f(A))$03 as the squared-exponential kernel approaches the identity in the limit $\operatorname{tr}(f(A))$04 (Zvonek et al., 2023). ContHutch++ then combines a continuous randomized range finder with projected residual estimation: $\operatorname{tr}(f(A))$05 Its main significance is not only the Hutch++-style $\operatorname{tr}(f(A))$06 relative-error scaling, but also the avoidance of spectral artifacts introduced by discretization, including non-converged eigenvalues, spectral pollution, and spectral invisibility (Zvonek et al., 2023).

6. Dynamic estimation, application-specific adaptations, and limitations

Dynamic trace estimation exploits temporal smoothness. If $\operatorname{tr}(f(A))$07, DeltaShift updates

$\operatorname{tr}(f(A))$08

and achieves total matrix-vector complexity

$\operatorname{tr}(f(A))$09

(Dharangutte et al., 2021). In the natural regime $\operatorname{tr}(f(A))$10, this is a quadratic improvement over repeated Hutchinson. A broader framework for slowly varying sequences replaces the older $\operatorname{tr}(f(A))$11 dependence by the path-length-style bound

$\operatorname{tr}(f(A))$12

and can estimate local changes on the fly with nearly no added cost in certain cases (Gokhale et al., 22 Jun 2026).

Application-specific reformulations show how elastic the trace-estimation viewpoint has become. In finite-width neural tangent kernel analysis, the empirical NTK is treated as a matrix-free operator

$\operatorname{tr}(f(A))$13

and Hutch++ is used to estimate $\operatorname{tr}(f(A))$14, $\operatorname{tr}(f(A))$15, alignment, and effective rank without explicitly forming the kernel (Hazelden, 13 Nov 2025). Because

$\operatorname{tr}(f(A))$16

one-sided estimators using only reverse-mode or only forward-mode automatic differentiation are also possible (Hazelden, 13 Nov 2025). This suggests that implicit trace estimation is increasingly functioning as an interface between randomized numerical linear algebra and matrix-free automatic differentiation.

Other adaptations are more heuristic. Machine-learned probing vectors replace random noise by a small trained set $\operatorname{tr}(f(A))$17, producing estimators of the form

$\operatorname{tr}(f(A))$18

with optional bias correction; in numerical experiments, $\operatorname{tr}(f(A))$19 learned probing vectors were similar in precision to $\operatorname{tr}(f(A))$20 random noise vectors, but the estimator can be biased and depends on representative training data (Yoon, 2016). Approximate-inverse interpolation instead fits the diagonal of $\operatorname{tr}(f(A))$21 from the diagonal of an inexpensive approximate inverse and can serve either as a standalone estimator or as variance reduction for Monte Carlo (Wu et al., 2015).

Not all favorable asymptotics survive in practice. A lattice-QCD study of $\operatorname{tr}(f(A))$22 with $\operatorname{tr}(f(A))$23-improved Wilson fermions found no noise reduction from Hutch++ or XTrace at a moderate number of sources; up to about $\operatorname{tr}(f(A))$24 inversions, neither method beat the Girard–Hutchinson baseline, and the deflated contribution remained around $\operatorname{tr}(f(A))$25 (Cotellucci et al., 2023). Likewise, stochastic probing can outperform Hutch++ and XTrace when $\operatorname{tr}(f(A))$26 is not well approximated by a low-rank matrix, whereas low-rank methods become superior when the spectrum decays rapidly (Frommer et al., 2023). The central misconception corrected by these results is that a better asymptotic rate automatically yields a better estimator at modest budgets. In implicit trace estimation, oracle model, spectral decay, locality, and update structure are all first-order determinants of performance.