XNysTrace: Unbiased Nyström Trace Estimator

• XNysTrace is an unbiased estimator that combines Nyström low-rank approximation, residual correction, and leave-one-out symmetrization for PSD trace estimation.
• It reuses random probes exchangeably to achieve significantly lower error—often orders of magnitude better than Hutch++—within the same matvec budget.
• Ideal for PSD matrices with rapidly decaying eigenvalues, it balances low computational cost with high accuracy by exploiting spectral properties.

XNysTrace is an unbiased, exchangeable, Nyström-based stochastic trace estimator for symmetric positive semidefinite (PSD) matrices, introduced in "XTrace: Making the most of every sample in stochastic trace estimation" (Epperly et al., 2023). It addresses the implicit trace estimation problem: given a large matrix $A \in \mathbb{R}^{N\times N}$ that can be accessed only through matrix–vector products $w \mapsto Aw$, approximate $\operatorname{tr}(A)$ under a fixed budget of matvecs. Its defining feature is the combination of Nyström low-rank approximation, residual trace correction, and leave-one-out symmetrization, so that all random probes are reused in an exchangeable way. In the reported numerical experiments, this yields errors that are often orders of magnitude smaller than those of Girard–Hutchinson, Hutch++, and Nyström++ at the same matvec budget (Epperly et al., 2023).

1. Problem setting and PSD specialization

The relevant setting for XNysTrace is the implicit or stochastic trace estimation problem for a PSD matrix $A$, where the dominant cost is the number of matrix–vector products. The objective is: with a budget of $m$ matvecs, produce an accurate estimator of $\operatorname{tr}(A)$. This is a classical setting in randomized numerical linear algebra, but XNysTrace is specialized to the PSD case, where Nyström approximation is particularly efficient (Epperly et al., 2023).

The canonical baseline is the Girard–Hutchinson estimator,

$$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$

for isotropic random vectors $w_i$ satisfying $\mathbb{E}[w_i w_i^\ast]=I$. This estimator is unbiased, but its variance decays only at the Monte Carlo rate $O(m^{-1})$. Variance-reduced estimators improve on this by approximating a dominant low-rank component of $w \mapsto Aw$0 and applying stochastic estimation only to the residual trace.

For PSD matrices, Nyström approximation is attractive because a rank-$w \mapsto Aw$1 approximation can be formed using only $w \mapsto Aw$2 matvecs, whereas a randomized SVD typically requires $w \mapsto Aw$3 matvecs. The standard Nyström approximation based on a test matrix $w \mapsto Aw$4 is

$w \mapsto Aw$5

This economy in matvecs is one of the main reasons XNysTrace can exploit a fixed budget more aggressively than SVD-based competitors in rapidly decaying PSD spectra.

2. Estimator construction

XNysTrace replaces a fixed sketch/residual split by a leave-one-out symmetrization over all probes (Epperly et al., 2023). Given $w \mapsto Aw$6 and a matvec budget $w \mapsto Aw$7, the construction is:

1. Draw iid isotropic random vectors

$w \mapsto Aw$8

The theory uses $w \mapsto Aw$9. Reported experiments also use signs and normalized Gaussian or spherical vectors.

2. Form

$\operatorname{tr}(A)$0

3. For each $\operatorname{tr}(A)$1, remove the $\operatorname{tr}(A)$2-th column:

$\operatorname{tr}(A)$3

4. Build the leave-one-out Nyström approximation

$\operatorname{tr}(A)$4

5. Form the basic estimator

$\operatorname{tr}(A)$5

6. Average all basic estimators:

$\operatorname{tr}(A)$6

The associated empirical error estimate is

$\operatorname{tr}(A)$7

The estimator is unbiased because $\operatorname{tr}(A)$8 is independent of $\operatorname{tr}(A)$9, which depends only on $A$0. Conditionally on $A$1,

$A$2

and therefore

$A$3

Averaging preserves unbiasedness while reusing all probes.

3. Exchangeability and full sample reuse

The organizing principle behind XNysTrace is that if the probes $A$4 are exchangeable, then the minimum-variance unbiased estimator should be a symmetric function of them (Epperly et al., 2023). Standard deflation methods violate that principle by assigning different subsets of iid probes to different jobs. In Hutch++ or Nyström++, some vectors are designated as sketch vectors and others as residual probes; that asymmetry discards information that is available under the same sampling distribution.

XNysTrace enforces symmetry by cycling through all choices of the held-out residual probe and averaging the resulting leave-one-out estimators. Every probe contributes to the low-rank approximation in all iterations except one, and every probe also serves once as an independent residual probe. A common misconception is that XNysTrace is merely “Nyström + residual correction.” Its distinctive feature is Nyström approximation plus residual correction plus full sample reuse through symmetrization.

The distinction from related estimators can be summarized as follows:

Method	Matrix class	Sample usage
Girard–Hutchinson	general square matrix	quadratic-form probes only
Hutch++	general square matrix	fixed sketch/residual split
Nyström++	PSD	fixed Nyström/residual split
XTrace	general square matrix	exchangeable low-rank residual averaging
XNysTrace	PSD	exchangeable leave-one-out Nyström averaging

The $A$5 are not independent, because their leave-one-out Nyström approximations share nearly the same sample set. The analysis therefore controls not only the variance of each basic estimator but also their covariances. The favorable net variance comes from the fact that dominant eigenspaces are removed by the Nyström term, so the stochastic correction acts mainly on the spectral tail.

4. Implementation, rank-one updates, and normalization

A naive implementation would recompute $A$6 Nyström approximations from scratch. XNysTrace avoids that. It can be implemented with exactly $A$7 matvecs with $A$8 and $A$9 post-processing arithmetic, so its matvec complexity is comparable to that of Hutch++ (Epperly et al., 2023).

Let

$m$0

Compute a thin QR factorization

$m$1

and define

$m$2

Then the Nyström approximation based on all columns is

$m$3

If $m$4 is obtained from $m$5 by deleting row and column $m$6, then

$m$7

The key identity is the rank-one update

$m$8

where $m$9 is the $\operatorname{tr}(A)$0-th standard basis vector in $\operatorname{tr}(A)$1. This permits inexpensive leave-one-out updates from full-sample quantities.

For numerical stability, the MATLAB implementation perturbs $\operatorname{tr}(A)$2 before QR or Cholesky steps via

$\operatorname{tr}(A)$3

The reported implementation also recommends a normalized variant using spherically symmetric vectors. If

$\operatorname{tr}(A)$4

then

$\operatorname{tr}(A)$5

and one uses

$\operatorname{tr}(A)$6

This normalization reduces variance arising from random probe lengths, especially when the spectrum contains flat regions.

5. Unbiasedness, variance bounds, and spectral dependence

Theorem 1.1 states that XNysTrace is unbiased: $\operatorname{tr}(A)$7 A central variance estimate is Proposition 5.3: for PSD $\operatorname{tr}(A)$8 and XNysTrace with $\operatorname{tr}(A)$9 standard normal probes,

$$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$0

This ties trace-estimation error directly to the Frobenius error of the leave-one-out Nyström approximation (Epperly et al., 2023).

The more detailed theorem refines this into an explicit spectrum-sensitive bound in terms of a best rank-$$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$1 approximation $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$2, the trace norm tail $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$3, the Frobenius tail $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$4, and the gap $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$5, for $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$6, $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$7, and $$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$8. The structural implication is that XNysTrace inherits its accuracy from how well Nyström compresses the dominant PSD eigenspaces under the available probe budget.

For exponentially decaying eigenvalues,

$$\widehat{\operatorname{tr}}_{\mathrm{GH}} = \frac{1}{m}\sum_{i=1}^m w_i^\ast A w_i,$$9

the reported asymptotic comparison is

$w_i$0

$w_i$1

$w_i$2

Accordingly, XTrace improves the exponent from $w_i$3 to $w_i$4 relative to Hutch++, and XNysTrace improves it further to $w_i$5, described in the paper as a $w_i$6 exponential-rate improvement over Hutch++.

At the level of variance order, Girard–Hutchinson remains at $w_i$7, whereas Hutch++, XTrace, and XNysTrace attain $w_i$8. The paper also states that for all three improved methods, $w_i$9 matvecs suffice to achieve variance at most $\mathbb{E}[w_i w_i^\ast]=I$0, with XNysTrace having the strongest spectrum-sensitive guarantee.

6. Empirical behavior, favorable regimes, and limitations

The experimental study compares Hutch, LRA, Hutch++, Nyström++, XTrace, and XNysTrace on synthetic PSD matrices of the form $\mathbb{E}[w_i w_i^\ast]=I$1 with $\mathbb{E}[w_i w_i^\ast]=I$2, using flat, polynomial-decay, exponential-decay, and step spectra, and reporting average relative error over 1000 trials (Epperly et al., 2023). Variance-reduced methods outperform Hutch except on the flat spectrum. XNysTrace is frequently the most accurate method; on exponential decay it converges at about $\mathbb{E}[w_i w_i^\ast]=I$3 the rate of XTrace and Nyström++, and $\mathbb{E}[w_i w_i^\ast]=I$4 the rate of Hutch++. On the step spectrum, XNysTrace and Nyström++ can perform poorly because the relevant error bound depends on the slowly decaying nuclear norm tail, though normalization improves behavior.

A second experiment estimates the partition function

$\mathbb{E}[w_i w_i^\ast]=I$5

for the transverse-field Ising model on 18 sites, where $\mathbb{E}[w_i w_i^\ast]=I$6. In that setting matvecs with $\mathbb{E}[w_i w_i^\ast]=I$7 dominate runtime. With only $\mathbb{E}[w_i w_i^\ast]=I$8 matvecs, all variance-reduced methods are already about 5 orders of magnitude better than Hutch. As $\mathbb{E}[w_i w_i^\ast]=I$9 increases, XTrace and XNysTrace outperform Hutch++; at $O(m^{-1})$0, XTrace is about $O(m^{-1})$1 more accurate than Hutch++, and XNysTrace is about $O(m^{-1})$2 more accurate. The explanation given is a sharp eigenvalue drop at indices $O(m^{-1})$3–$O(m^{-1})$4,

$O(m^{-1})$5

which XNysTrace can exploit earlier because its effective approximation rank is roughly $O(m^{-1})$6, compared with roughly $O(m^{-1})$7 for XTrace and $O(m^{-1})$8 for Hutch++.

The paper also reports that the posterior error estimate

$O(m^{-1})$9

tracks the true error well in the partition-function experiment, within about a factor of $w \mapsto Aw$00.

The method is most favorable when $w \mapsto Aw$01 is PSD, effectively low rank, and has rapidly decaying eigenvalues so that a relatively small dominant eigenspace carries most of the trace. It is less favorable when the spectrum has a long slowly decaying tail, when $w \mapsto Aw$02 remains large, or when post-processing cost matters more than matvec count. The supplementary runtime discussion therefore notes that XNysTrace’s post-processing can be higher than XTrace’s, so when matvecs are cheap, XTrace may be preferable. On PSD problems with strong spectral decay, however, XNysTrace is usually the more accurate of the two.