Hutchinson Diagonal Estimator

• The Hutchinson diagonal estimator is a stochastic method that approximates the diagonal of a matrix using randomized Rademacher vectors and matrix–vector products.
• Its error is governed by the off-diagonal entries, yielding a dimension-free high-probability bound with the ℓ₂-error scaling as √(ln(2/δ)/m) times the Frobenius norm of the off-diagonals.
• The method efficiently operates in scenarios with large, sparse, or implicitly defined matrices, making it a practical tool for matrix-free inference and related applications.

Searching arXiv for relevant papers on the Hutchinson diagonal estimator and closely related methods. arxiv_search query: "Hutchinson diagonal estimator" max_results: 10 The Hutchinson diagonal estimator is a stochastic method for approximating the diagonal of a matrix from matrix–vector products alone. For a matrix $A\in\mathbb R^{n\times n}$, possibly non-symmetric, it replaces explicit access to $\operatorname{diag}(A)$ by randomized probing with i.i.d. Rademacher vectors and averages the entrywise products $g\odot(Ag)$. The estimator is unbiased, its error is governed by the off-diagonal mass of $A$, and a tight high-probability analysis shows that its $\ell_2$-error scales as $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$ with no dependence on the ambient dimension $n$, improving earlier bounds by a $\log(n)$ factor (Dharangutte et al., 2022).

1. Formal definition

Let

$$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$

For $z=1,\dots,m$, draw independent Rademacher vectors

$\operatorname{diag}(A)$0

with i.i.d. entries. Hutchinson’s diagonal estimator is

$\operatorname{diag}(A)$1

Its basic identity is immediate: $\operatorname{diag}(A)$2 Equivalently, at the coordinate level,

$\operatorname{diag}(A)$3

This is the diagonal analogue of the Girard–Hutchinson trace identity $\operatorname{diag}(A)$4, and it places diagonal estimation in the same matrix-free oracle model as stochastic trace estimation (Dharangutte et al., 2022).

Operationally, one iteration requires one matrix–vector multiply $\operatorname{diag}(A)$5 and one Hadamard product. The estimator therefore remains applicable when $\operatorname{diag}(A)$6 is large, sparse, or only implicitly defined through a matvec oracle.

2. Error decomposition and dependence on off-diagonal structure

The estimator’s randomness is driven entirely by off-diagonal interactions. If

$\operatorname{diag}(A)$7

then the per-sample error satisfies

$\operatorname{diag}(A)$8

Thus the squared error of a single diagonal probe is exactly a Hutchinson trace estimator applied to the positive semidefinite matrix

$\operatorname{diag}(A)$9

Since

$g\odot(Ag)$0

the Frobenius norm of the zero-diagonal part $g\odot(Ag)$1 is the fundamental scale parameter for the analysis (Dharangutte et al., 2022).

This characterization clarifies why diagonal recovery is easy when the off-diagonal mass is small and hard when it is large. The diagonal itself does not contribute to the stochastic error, because the estimator reproduces diagonal entries exactly in expectation and only the terms with $g\odot(Ag)$2 fluctuate.

Coordinatewise variance statements make the same point. For the classical diagonal estimator,

$g\odot(Ag)$3

and in the related Girard–Hutchinson formulation,

$g\odot(Ag)$4

These formulas show that stochastic diagonal estimation is not controlled by a global spectral quantity alone; it is controlled by the off-diagonal geometry of the rows and columns being probed (Tsyganov et al., 6 Aug 2025).

3. Tight high-probability analysis

The central result of the modern analysis is a dimension-free concentration bound. There is a universal constant $g\odot(Ag)$5 such that for any $g\odot(Ag)$6 and any $g\odot(Ag)$7,

$g\odot(Ag)$8

The failure probability enters only through $g\odot(Ag)$9, and the estimate does not depend on $A$0 (Dharangutte et al., 2022).

Earlier analyses, including Baston and Nakatsukasa (2022), gave a related guarantee with an extra factor of $A$1,

$A$2

so the newer result removes the ambient-dimension term entirely. In this sense, the analysis is tight with respect to dimension dependence.

A recurrent misconception is that such a bound should follow by estimating each coordinate separately and applying a union bound. That strategy is precisely what incurs the $A$3 penalty. The dimension-free result shows that the estimator must instead be treated as a genuinely vector-valued random object, with concentration controlled at the level of its full Euclidean norm rather than coordinatewise maxima (Dharangutte et al., 2022).

Another common misunderstanding is to read “independent of $A$4” as “free of dimensional cost.” The theorem removes $A$5 from the probability bound, but each iteration still requires one matvec and $A$6 entrywise work, so computational cost remains linear in the vector dimension per sample.

4. Proof architecture

The proof strategy proceeds by connecting diagonal estimation to trace estimation and then avoiding coordinatewise decoupling. For one sample,

$A$7

so the squared error is a trace-estimator random variable for a positive semidefinite matrix. Since $A$8 is known to have exponential concentration, via the Hanson–Wright inequality or a moment-generating-function bound, each $A$9 is sub-exponential with mean $\ell_2$0 (Dharangutte et al., 2022).

The essential difficulty is then to control

$\ell_2$1

without paying a union-bound penalty. The analysis uses three ingredients.

First, it symmetrizes the sum and compares its moments to those of a scalar sum of i.i.d. sub-Gaussian variables. The data summarize this step as a coupling of inner products to absolute values.

Second, it shows that the squared norm

$\ell_2$2

has a moment-generating function bounded like a sub-exponential random variable with parameter $\ell_2$3.

Third, it applies a scalar tail bound for sub-exponential variables to conclude the vector-valued high-probability estimate.

The argument draws explicitly on a lemma attributed to Vershynin for the trace-estimator error

$\ell_2$4

namely

$\ell_2$5

together with a moment-comparison plus symmetrization step in the style of Yurinskii’s vector Bernstein argument, and standard sub-Gaussian and sub-exponential tail inequalities (Dharangutte et al., 2022).

5. Sample complexity, oracle cost, and tightness

The high-probability theorem yields a direct sample-complexity prescription. To guarantee

$\ell_2$6

with probability at least $\ell_2$7, it is sufficient to choose

$\ell_2$8

Beyond the choice of $\ell_2$9, the estimator has no additional parameters (Dharangutte et al., 2022).

Its computational footprint is correspondingly simple. Each iteration requires one matrix–vector multiply $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$0, with cost $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$1, and $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$2 work to form the Hadamard product. The total cost is therefore

$\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$3

For large sparse or implicitly defined matrices, this makes the method inexpensive relative to explicit diagonal formation.

The dependence on $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$4 and $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$5 is not an artifact of the proof. A $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$6 example establishes tightness: $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$7 For this matrix, the estimator yields $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$8, where $\sqrt{\ln(2/\delta)/m}\,\|\bar A\|_F$9 is the sum of $n$0 independent Rademacher products $n$1. By the tightness of the Chernoff bound, achieving $n$2 with probability at least $n$3 requires

$n$4

Accordingly, the rate $n$5 is optimal for Hutchinson’s estimator in general (Dharangutte et al., 2022).

6. Variants, extensions, and adjacent research directions

The estimator has become a template for several distinct lines of work: variance reduction for matrix diagonals, matrix-free norm estimation, tensor generalization, and quantum constructions that reproduce the same first- and second-moment identities.

Direction	Mechanism	Stated result
XDIAG	Leave-one-out low-rank control variate with exchangeability	Unbiased; variance $n$6
TwINEst / TwINEst++	Diagonal estimation on $n$7 plus exact row-norm evaluation	Matrix-free estimation of $n$8 and $n$9
Tensor estimator	Hadamard product of $\log(n)$0 probe vectors with a tensor–vector query	Unbiased diagonal and trace estimation for $\log(n)$1-order tensors
Quantum diagonal state designs	Random-phase-style states generated by diagonal Hamiltonians	Reproduces unbiased Girard–Hutchinson trace/diagonal estimators

XDIAG, developed in the XTrace framework, keeps the same target $\log(n)$2 but replaces plain Monte Carlo averaging by a leave-one-out low-rank control variate. In the summarized analysis, each basic estimate remains unbiased because the probe used in the residual term is independent of the leave-$\log(n)$3 approximation $\log(n)$4, and the exchangeability argument reduces the mean-square error from the standard $\log(n)$5 Monte Carlo rate to $\log(n)$6. When the spectrum decays exponentially, the stated bound becomes exponential in $\log(n)$7 (Epperly et al., 2023).

TwINEst and TwINEst++ adapt the diagonal estimator to norm estimation in the matrix-free setting. The key identity is

$\log(n)$8

A noisy diagonal estimate of $\log(n)$9 is formed first, then the algorithm identifies an index $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$0, computes the exact row norm $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$1, and returns that exact quantity. TwINEst++ further adds a Hutch++-style low-rank correction. The stated oracle-complexity results show exact recovery with high probability once the diagonal estimation error is below the row-gap threshold $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$2 (Tsyganov et al., 6 Aug 2025).

The tensor generalization extends the same idea beyond matrices. For an $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$3-order tensor $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$4, one draws $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$5 independent probe vectors $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$6, computes

$$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$7

forms the entrywise product

$$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$8

and obtains $$\operatorname{diag}(A)=(A_{11},A_{22},\dots,A_{nn})^\top\in\mathbb R^n, \qquad \bar A_{ii}=0,\quad \bar A_{ij}=A_{ij}\ \text{for }i\neq j.$$9. When $z=1,\dots,m$0, the construction reduces exactly to the Bekas–Saad diagonal estimator and the classical Hutchinson trace estimator (Verma et al., 25 Oct 2025).

A conceptually different extension appears in quantum information. There, random-phase states and diagonal state $z=1,\dots,m$1-designs are used to reproduce the same moment identities underlying the classical Girard–Hutchinson formulas. The summarized result constructs a diagonal state 3-design from real-time evolution under 2-local Hamiltonians, with randomness arising from stochastic durations. Because the ensemble matches the necessary $z=1,\dots,m$2 moments, it preserves the unbiased trace and diagonal estimators and their variance formulas (Shen et al., 2024).

Taken together, these developments show that the Hutchinson diagonal estimator is not merely a special-purpose Monte Carlo device for matrix diagonals. It is a structural primitive for matrix-free inference, one whose unbiasedness survives in several generalized oracle models, while its concentration behavior can be sharpened, variance-reduced, or repurposed depending on the surrounding algorithmic objective.