---
title: Hutchinson Diagonal Estimator
url: https://www.emergentmind.com/topics/hutchinson-diagonal-estimator
type: topic
---

# Hutchinson Diagonal Estimator

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 [2208.03268].

## 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
\[
g^{(z)}\in\{-1,+1\}^n,\qquad
\Pr[g^{(z)}_i=+1]=\Pr[g^{(z)}_i=-1]=\tfrac12,
\]
with i.i.d. entries. Hutchinson’s diagonal estimator is
\[
\tilde d
=
\frac1m\sum_{z=1}^m g^{(z)}\odot\bigl(A\,g^{(z)}\bigr).
\]

Its basic identity is immediate:
\[
\mathbb E[\tilde d]=\operatorname{diag}(A).
\]
Equivalently, at the coordinate level,
\[
\mathbb E\!\left[g_i(Ag)_i\right]=A_{ii}.
\]
This is the diagonal analogue of the Girard–Hutchinson trace identity \(\mathbb E[g^\top A g]=\operatorname{Tr}A\), and it places diagonal estimation in the same matrix-free oracle model as stochastic trace estimation [2208.03268].

Operationally, one iteration requires one matrix–vector multiply \(A g\) and one Hadamard product. The estimator therefore remains applicable when \(A\) 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
\[
e = g\odot(A g)-\operatorname{diag}(A),
\]
then the per-sample error satisfies
\[
\|e\|_2^2
=
g^\top(\bar A^\top \bar A)g.
\]
Thus the squared error of a single diagonal probe is exactly a Hutchinson trace estimator applied to the positive semidefinite matrix
\[
B=\bar A^\top \bar A.
\]
Since
\[
\mathbb E[g^\top B g]=\operatorname{tr}(B)=\|\bar A\|_F^2,
\]
the Frobenius norm of the zero-diagonal part \(\bar A\) is the fundamental scale parameter for the analysis [2208.03268].

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 \(i\neq j\) fluctuate.

Coordinatewise variance statements make the same point. For the classical diagonal estimator,
\[
\operatorname{Var}\!\bigl[D_i^1(B)\bigr]=\sum_{j\neq i} B_{ij}^2,
\]
and in the related Girard–Hutchinson formulation,
\[
\operatorname{Var}[v^\top A v]=\sum_{i\neq j}|A_{ij}|^2\le \|A\|_F^2,
\qquad
\operatorname{Var}[v_i(Av)_i]\le \sum_{j\neq i}|A_{ij}|^2.
\]
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 [2508.04444].

## 3. Tight high-probability analysis

The central result of the modern analysis is a dimension-free concentration bound. There is a universal constant \(c>0\) such that for any \(\delta\in(0,1]\) and any \(m\ge 1\),
\[
\Pr\!\left[
\|\tilde d-\operatorname{diag}(A)\|_2
\le
c\,\sqrt{\frac{\ln(2/\delta)}{m}}\,\|\bar A\|_F
\right]
\ge 1-\delta.
\]
The failure probability enters only through \(\ln(1/\delta)\), and the estimate does not depend on \(n\) [2208.03268].

Earlier analyses, including Baston and Nakatsukasa (2022), gave a related guarantee with an extra factor of \(\sqrt{\ln n}\),
\[
\|\tilde d-\operatorname{diag}(A)\|_2
=
O\!\left(\sqrt{\frac{\ln(n/\delta)}{m}}\right)\,\|\bar A\|_F,
\]
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 \(\ln n\) 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 [2208.03268].

Another common misunderstanding is to read “independent of \(n\)” as “free of dimensional cost.” The theorem removes \(n\) from the probability bound, but each iteration still requires one matvec and \(O(n)\) 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,
\[
\|e\|_2^2 = g^\top(\bar A^\top\bar A)g = T(\bar A^\top\bar A),
\]
so the squared error is a trace-estimator random variable for a positive semidefinite matrix. Since \(T(B)\) is known to have exponential concentration, via the Hanson–Wright inequality or a moment-generating-function bound, each \(\|e\|_2^2\) is sub-exponential with mean \(\|\bar A\|_F^2\) [2208.03268].

The essential difficulty is then to control
\[
\left\|\frac1m\sum_{z=1}^m e_z\right\|_2
\]
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
\[
\left\|\sum_z e_z\right\|_2^2
\]
has a moment-generating function bounded like a sub-exponential random variable with parameter \(m\,\|\bar A\|_F^2\).

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
\[
Z(B)=g^\top B g-\operatorname{tr}(B),
\]
namely
\[
\mathbb E[e^{\lambda Z(B)}]\le e^{C\lambda^2\|B\|_F^2}
\quad\text{for }\ |\lambda|\le c/\|B\|_2,
\]
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 [2208.03268].

## 5. Sample complexity, oracle cost, and tightness

The high-probability theorem yields a direct sample-complexity prescription. To guarantee
\[
\|\tilde d-\operatorname{diag}(A)\|_2\le \epsilon\,\|\bar A\|_F
\]
with probability at least \(1-\delta\), it is sufficient to choose
\[
m \gtrsim \frac{\ln(2/\delta)}{\epsilon^2}.
\]
Beyond the choice of \(m\), the estimator has no additional parameters [2208.03268].

Its computational footprint is correspondingly simple. Each iteration requires one matrix–vector multiply \(A g\), with cost \(T_{mv}\), and \(O(n)\) work to form the Hadamard product. The total cost is therefore
\[
m\,(T_{mv}+O(n)).
\]
For large sparse or implicitly defined matrices, this makes the method inexpensive relative to explicit diagonal formation.

The dependence on \(\ln(1/\delta)\) and \(\epsilon^{-2}\) is not an artifact of the proof. A \(2\times 2\) example establishes tightness:
\[
A=\begin{pmatrix}0 & 1\\ 0 & 0\end{pmatrix}.
\]
For this matrix, the estimator yields \(\tilde d=(S/m,0)\), where \(S\) is the sum of \(m\) independent Rademacher products \(g_1g_2\). By the tightness of the Chernoff bound, achieving \(|S/m|\le \epsilon\) with probability at least \(1-\delta\) requires
\[
m=\Omega\!\left(\frac{\ln(1/\delta)}{\epsilon^2}\right).
\]
Accordingly, the rate \(O(\ln(1/\delta)/\epsilon^2)\) is optimal for Hutchinson’s estimator in general [2208.03268].

## 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 \(O(1/m^2)\) |
| TwINEst / TwINEst++ | Diagonal estimation on \(AA^\top\) plus exact row-norm evaluation | Matrix-free estimation of \(\|A\|_{2\to\infty}\) and \(\|A\|_{1\to2}\) |
| Tensor estimator | Hadamard product of \(N-1\) probe vectors with a tensor–vector query | Unbiased diagonal and trace estimation for \(N\)-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 \(\operatorname{diag}(A)\) 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-\(i\) approximation \(A_{(i)}\), and the exchangeability argument reduces the mean-square error from the standard \(O(1/m)\) Monte Carlo rate to \(O(1/m^2)\). When the spectrum decays exponentially, the stated bound becomes exponential in \(m\) [2301.07825].

TwINEst and TwINEst++ adapt the diagonal estimator to norm estimation in the matrix-free setting. The key identity is
\[
\|A\|_{2\to\infty}^2=\max_i \operatorname{diag}(AA^\top)_i.
\]
A noisy diagonal estimate of \(AA^\top\) is formed first, then the algorithm identifies an index \(j=\arg\max_i D_i\), computes the exact row norm \(\|A_j\|_2=\|A^\top e_j\|_2\), 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 \(\Delta/2\) [2508.04444].

The tensor generalization extends the same idea beyond matrices. For an \(N\)-order tensor \(\mathcal T\), one draws \(N-1\) independent probe vectors \(g^{(1)},\dots,g^{(N-1)}\), computes
\[
v=\mathcal T\times_1 g^{(1)}\times_2 g^{(2)}\cdots\times_{N-1} g^{(N-1)},
\]
forms the entrywise product
\[
y=(g^{(1)}*\cdots * g^{(N-1)})*v,
\]
and obtains \(\mathbb E[y_i]=a_{i,\dots,i}\). When \(N=2\), the construction reduces exactly to the Bekas–Saad diagonal estimator and the classical Hutchinson trace estimator [2510.22157].

A conceptually different extension appears in quantum information. There, random-phase states and diagonal state \(t\)-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 \(k=1,2,3\) moments, it preserves the unbiased trace and diagonal estimators and their variance formulas [2401.04176].

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.

Source: https://www.emergentmind.com/topics/hutchinson-diagonal-estimator