# Updates for Hutch++

At SOSA 2021, I presented the paper Hutch++: Optimal Stochastic Trace Estimation (with Cameron Musco, Christopher Musco, and David P. Woodruff). This webpage contains some short updates that may be useful to people.

# Hutch++ for Undergrads

I believe that the analysis of both Hutchinson's Estimator and Hutch++ is simple and approachable. This section contains the analysis an advanced undergraduate should hopefully be able to understand.

## Definition 1: Hutchinson's Estimator Definition

Fix $\boldsymbol{ A} \in \mathbb R^{n \times n}$ be a real matrix, and fix $\ell\in\mathbb N$. Let $\mathbf{ x}_1,\ldots,\mathbf{ x}_\ell \in\mathbb R^n$ be vectors drawn with i.i.d. $\mathcal N(0,1)$ entries. Then, Hutchinson's Estimator is

$H_\ell(\boldsymbol{ A}) \;{\vcentcolon=}\; \frac1\ell \sum_{i=1}^\ell \mathbf{ x}_i^\intercal \boldsymbol{ A} \mathbf{ x}_i$

We can easily verify two properties of this classical estimator for symmetric $\boldsymbol{ A}$:

## Lemma 1: Hutchinson's Mean and Variance Lemma

$\mathbb E[H_\ell(\boldsymbol{ A})] = \text{tr}(\boldsymbol{ A})$ and $\text{Var}[H_\ell(\boldsymbol{ A})]=\frac2\ell\| \boldsymbol{ A}\|_{ F}^2$

Proof. The expectation is easy to verify:

$\mathbb E[H_\ell(\boldsymbol{ A})] = \mathbb E[\mathbf{ x}^\intercal\boldsymbol{ A}\mathbf{ x}] = \sum_{i=1}^n \sum_{j=1}^n \mathbb E[[\boldsymbol{ A}]_{i,j}\mathbf{ x}_i\mathbf{ x}_j] = \sum_{i=1}^n \boldsymbol{ A}_{i,i} = \text{tr}(\boldsymbol{ A})$

In order to analyze the variance, we present the analysis of Lemma 9 from Avron, Toledo (2011). Let $\boldsymbol{ A}=\boldsymbol{ U}\boldsymbol{ \Lambda}\boldsymbol{ U}^\intercal$ be the eigendecomposition of $\boldsymbol{ A}$, and let $\mathbf{ y} \;{\vcentcolon=}\; \boldsymbol{ U}^\intercal\mathbf{ x}$. Then, for $\mathbf{ x}\sim\mathcal N(0,\boldsymbol{ I})$, we have $\mathbf{ y} \sim\mathcal N(0,\boldsymbol{ U}\boldsymbol{ U}^\intercal) = \mathcal N(0,\boldsymbol{ I})$. Therefore,

\begin{aligned} \text{Var}[H_\ell(\boldsymbol{ A})] &= \frac1\ell \text{Var}[\mathbf{ x}^\intercal\boldsymbol{ A}\mathbf{ x}] \\ &= \frac1\ell \text{Var}[\mathbf{ x}^\intercal\boldsymbol{ U}\boldsymbol{ \Lambda}\boldsymbol{ U}^\intercal\mathbf{ x}] \\ &= \frac1\ell \text{Var}[(\boldsymbol{ U}^\intercal\mathbf{ x})^\intercal\boldsymbol{ \Lambda}(\boldsymbol{ U}\mathbf{ x})] \\ &= \frac1\ell \text{Var}[\mathbf{ y}^\intercal\boldsymbol{ \Lambda}\mathbf{ y}] \\ &= \frac1\ell \text{Var}[\sum_{i=1}^n \lambda_i y_i^2] \\ &= \frac1\ell \sum_{i=1}^n \text{Var}[\lambda_i y_i^2] \\ &= \frac1\ell \sum_{i=1}^n 2\lambda_i^2 \\ &= \frac2\ell \| \boldsymbol{ A}\|_{ F}^2 \end{aligned}

Note we cannot just use linearity of trace naively on the sum $\mathbf{ x}^\intercal\boldsymbol{ A}\mathbf{ x}=\sum_{i=1}^n\sum_{j=1}^n x_ix_j \boldsymbol{ A}_{i,j}$ since linearity of variance only holds for independent samples.

$\blacksquare \, \,$

This variance is in term of the Frobenius norm of $\boldsymbol{ A}$. So, we expect that Hutchinson's Estimator satisfies $H_\ell(\boldsymbol{ A}) \in \text{tr}(\boldsymbol{ A}) \pm \frac{\sqrt2}{\sqrt\ell} \| \boldsymbol{ A}\|_{ F}$ with high probability. Typically in computer science we are given an error tolerance $\varepsilon$, and want to know how many matrix-vector products we need to compute. (i.e. how small can we make $\ell$ while still being accurate?) In order to build a guarantee here, we assume $\boldsymbol{ A}$ is positive semi-definite (PSD), which in turn implies that $\| \boldsymbol{ A}\|_{ F} \leq \text{tr}(\boldsymbol{ A})$. With this in mind, we can prove the standard result for Hutchinson's Estimator:

## Lemma 2: Hutchinson's Estimator Lemma

Fix $\boldsymbol{ A} \in \mathbb R^{d \times d}$ be a PSD matrix. Then, $\text{Var}[H_\ell(\boldsymbol{ A})]\leq\frac{2}{\ell}\text{tr}^2(\boldsymbol{ A})$

So, if we want to ensure $H_\ell(\boldsymbol{ A})\in(1\pm\varepsilon)\text{tr}(\boldsymbol{ A})$ with good probability, we need the standard deviation to have $\frac{\sqrt 2}{\sqrt \ell}\text{tr}(\boldsymbol{ A}) \leq \varepsilon \text{tr}(\boldsymbol{ A})$, so we need $\ell=O(\frac{1}{\varepsilon^2})$ samples.

## Part 2: Hutchinson's versus the Top Few Eigenvalues

This was a section in the first draft of the paper, but was unfortunately cut. This is a slow, methodical intuition for the geometry Hutch++ takes advantage of.

Hutchinson's Estimator is powerfull, and gives a nice rate of convergence. However, the proof of Lemma 2 has a suspicious step. If we write out three lines of analysis:

\begin{aligned} |\text{tr}(\boldsymbol{ A}) - H_\ell(\boldsymbol{ A})| &\leq {\textstyle\frac{ \sqrt 2}{ \sqrt \ell}} \| \boldsymbol{ A}\|_F & \hspace{1cm} & (\text{Standard Deviation})\\ &\leq {\textstyle\frac{ \sqrt 2}{ \sqrt \ell}} \text{tr}(\boldsymbol{ A}) & & (\| \boldsymbol{ A}\|_F \leq \text{tr}(\boldsymbol{ A})) \\ &= \varepsilon \ \text{tr}(\boldsymbol{ A}) & & (\ell = O({\textstyle\frac{ 1}{ \varepsilon}})) \\ \end{aligned}

Looking at these inequalities,

• The analysis of line 1 is always tight, since we know the variance of $H_\ell(\boldsymbol{ A})$ exactly.

• The analysis of line 3 is always tight, since we just set $\ell$ to the smallest value that gets error $\varepsilon$.

• The analysis of line 2 is not always tight.

The second line is tight only if $\| \boldsymbol{ A}\|_F \approx \text{tr}(\boldsymbol{ A})$. That is, the analysis above only tells us that Hutchinson's Estimator needs $O(\frac1\varepsilon)$ samples if $\boldsymbol{ A}$ is the kind of matrix that has $\| \boldsymbol{ A}\|_F \approx \text{tr}(\boldsymbol{ A})$. So, what kind of matrix is that?

## Claim 1: Eigenvalue Intuition Claim

If $\| \boldsymbol{ A}\|_F\approx\text{tr}(\boldsymbol{ A})$, then $\boldsymbol{ A}$ has very few large eigenvalues, followed by much smaller eigenvalues.

Try to show this yourself. If you're having trouble, look at the argument below.

Proof. First, we rewrite $\| \boldsymbol{ A}\|_F$ and $\text{tr}(\boldsymbol{ A})$ explicitly in terms of the eigenvalues of $\boldsymbol{ A}$. Let $\mathbf{ v} = [\lambda_1, \lambda_2, \ldots, \lambda_d]$ be the eigenvalues of $\boldsymbol{ A}$. Then, recall that

$\| \boldsymbol{ A}\|_F = \| \mathbf{ v}\|_2 \hspace{0.75cm} and \hspace{0.75cm} \text{tr}(\boldsymbol{ A}) = {\textstyle \sum_{i=1}^d \lambda_i} = \| \mathbf{ v}\|_1$

where the last equality follows from $\boldsymbol{ A}$ being PSD. So, we are given that $\| \mathbf{ v}\|_2 \approx \| \mathbf{ v}\|_1$: this is now just a statement about vectors in $\mathbb R^d$!

We now need to show that if a vector $\mathbf{ v}$ has $\| \mathbf{ v}\|_2 \approx \| \mathbf{ v}\|_1$, then $\mathbf{ v}$ is nearly sparse: it has only a few large values.

There are many intuition, and we will look at the extreme, and assume that $\| \mathbf{ v}\|_2=\| \mathbf{ v}\|_1$ exactly. Prove that $\mathbf{ v}$ has at most 1 nonzero element in it.

$\blacksquare \, \,$

We now finish the discussion with a comparison of two ideas:

• So, Hutchinson's Estimator only needs to use many samples (i.e. $O(\frac1\varepsilon)$ samples) if $\boldsymbol{ A}$ has very special structure: it has a small number of large eigenvalues.

• If a matrix has a small number of large eigenvalues, the trace must be well approximated by the sum of those eigenvalues.

With this, we conclude the fundamental intuition behind Hutch++:

## Claim 2: Hutch++ Intuition Claim

If $\boldsymbol{ A}$ is the kind of matrix that is hard for Hutchinsons's Estimator to handle, then $\text{tr}(\boldsymbol{ A})$ is well approximated by the top few eigenvalues of $\boldsymbol{ A}$

This leads us to pick the following rough design of an algorithm:

1. Find a good low-rank approximation $\tilde\boldsymbol{ A}_k$

2. Notice that $\text{tr}(\boldsymbol{ A}) = \text{tr}(\tilde \boldsymbol{ A}_k) + \text{tr}(\boldsymbol{ A} - \tilde\boldsymbol{ A}_k)$

3. Compute $\text{tr}(\tilde \boldsymbol{ A}_k)$ exactly

4. Approximate $\text{tr}(\boldsymbol{ A}-\tilde\boldsymbol{ A}_k)$ with Hutchinson's Estimator

5. Return $\text{Hutch++}(\boldsymbol{ A}) = \text{tr}(\tilde\boldsymbol{ A}_k) + H_\ell(\boldsymbol{ A}-\tilde\boldsymbol{ A}_k)$

In the next section, we state a formal version of Claim 1 (see Lemma 3), show how to compute such a matrix $\tilde\boldsymbol{ A}_k$ (see $\boldsymbol{ Q}\boldsymbol{ Q}^\intercal\boldsymbol{ A}$ in Theorem 1), and bound the complexity of the Hutch++ estimator (see Theorem 2).

## Part 3: The Variance of Hutch++

Special thanks to Joel A. Tropp for working most of this out, and encouraging us to look into sharpening the constants in the analysis of Hutch++.

By expressing the variance of Hutch++ with all of its constants laid bare, and by using a very simple analysis, this analysis will hopefully allow practitioners to easily the exact parameters in Hutch++ code.

Before bounding the variance of Hutch++, we include the proof of one lemma and import another theorem:

## Lemma 3: L2/L1/L0 Bounds Lemma

Let $\boldsymbol{ A}\in\mathbb R^{d \times d}$ be a PSD matrix. Let $\boldsymbol{ A}_k$ be the best rank-k approximation to $\boldsymbol{ A}$. Then, $\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F} \leq \frac{1}{2\sqrt k} \text{tr}(\boldsymbol{ A})$.

Proof. Let $\mathbf{ v}=[\lambda_1 \ldots \lambda_n]$ be the vector of eigenvalues of $\boldsymbol{ A}$, such that $\lambda_1 \geq \lambda_2 \geq ... \ \lambda_n \geq 0$. Then $\mathbf{ v}_k \;{\vcentcolon=}\; [\lambda_1 \ldots \lambda_k \ 0 \ldots 0]$ is the vector of eigenvalues of $\boldsymbol{ A}_k$. With this notation, we just need to prove that $\| \mathbf{ v}-\mathbf{ v}_k\|_{ 2} \leq \frac{1}{2\sqrt k} \| \mathbf{ v}\|_{ 1}$. The rest of this proof is just Lemma 7 from this paper.

Hölder's inequality states that $\sum_{i=1}^n a_i b_i \leq \| \mathbf{ a}\|_{ p} \cdot \| \mathbf{ b}\|_{ q}$ for any vector $\mathbf{ a},\mathbf{ b}$ so long as $\frac1p + \frac1q = 1$. If we take $\mathbf{ a}=\mathbf{ b}=\mathbf{ v}-\mathbf{ v}_k$, $p=1$, and $q=1$, we find that $\| \mathbf{ v}-\mathbf{ v}_k\|_{ 2}^2 \leq \| \mathbf{ v}-\mathbf{ v}_k\|_{ 1} \ \| \mathbf{ v}-\mathbf{ v}_k\|_{ \infty}$. Further, note that by our construction of $\mathbf{ v}_k$, $\| \mathbf{ v}\|_{ 1} = \| \mathbf{ v}_k\|_{ 1} + \| \mathbf{ v}-\mathbf{ v}_k\|_{ 1}$. Then,

$\frac{\| \mathbf{ v}-\mathbf{ v}_k\|_{ 2}}{\| \mathbf{ v}\|_{ 1}} \leq \frac{\sqrt{\| \mathbf{ v}-\mathbf{ v}_k\|_{ 1} \ \| \mathbf{ v}-\mathbf{ v}_k\|_{ \infty}}}{\| \mathbf{ v}_k\|_{ 1} + \| \mathbf{ v}-\mathbf{ v}_k\|_{ 1}}$

Note that $\| \mathbf{ v}-\mathbf{ v}_k\|_{ \infty} = \lambda_{k+1}$ and $\| \mathbf{ v}_k\|_{ 1} = \sum_{i=1}^k \lambda_i \geq k\lambda_{k+1}$. Then, if we let $\gamma \;{\vcentcolon=}\; \| \mathbf{ v}-\mathbf{ v}_k\|_{ 1}$, we find

$\frac{\sqrt{\| \mathbf{ v}-\mathbf{ v}_k\|_{ 1} \ \| \mathbf{ v}-\mathbf{ v}_k\|_{ \infty}}}{\| \mathbf{ v}_k\|_{ 1} + \| \mathbf{ v}-\mathbf{ v}_k\|_{ 1}} \leq \frac{\sqrt{\gamma \ \lambda_{k+1}}}{k\lambda_{k+1} + \gamma}$

We don't have a good way to bound $\gamma = \| \mathbf{ v}-\mathbf{ v}_k\|_{ 1}$, so we instead notice that the function $\gamma \mapsto \frac{\sqrt{a\gamma}}{b+\gamma}$ is maximized at $\gamma=b$ (check with calculus). So, $\frac{\sqrt{a\gamma}}{b+\gamma} \leq \frac{\sqrt{ab}}{2b}$ for all $\gamma$:

$\frac{\sqrt{\gamma \ \lambda_{k+1}}}{k\lambda_{k+1} + \gamma} \leq \frac{\sqrt{(k\lambda_{k+1}) \lambda_{k+1}}}{2k\lambda_{k+1}} = \frac{1}{2\sqrt k}$

So, overall, we have shown that

$\frac{\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F}}{\text{tr}(\boldsymbol{ A})} = \frac{\| \mathbf{ v}-\mathbf{ v}_k\|_{ 2}}{\| \mathbf{ v}\|_{ 1}} \leq \frac{1}{2\sqrt k}$
$\blacksquare \, \,$

We also import the following theorem, which controls the expected error from our low-rank approximation, and which we do not prove:

## Theorem 1: Expected Projection Error Theorem

Sample $\boldsymbol{ S}\in\mathbb R^{n \times (k+p)}$ with $\mathcal N(0,1)$ entries. Let $\boldsymbol{ Q}$ be any orthonormal basis for $\boldsymbol{ A}\boldsymbol{ S}$ (e.g. a QR decomposition). Then,

$\mathbb E[\| (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A}\|_{ F}^2] \leq (1+{\textstyle\frac{ k}{ p-1}})\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F}^2$

Proof. This is almost exactly the statement of Theorem 10.5 in , but not exactly. Theorem 10.5 bounds the expected Frobenius norm (not squared). Instead, at the bottom of page 57 and top of page 58 from this version of the paper, we see that the proof of Theorem 10.5 actually proves

$\left(\mathbb E [\| (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A}\|_{ F}^2]\right)^{1/2} \leq \ldots \leq (1+\frac{k}{p-1})\| \boldsymbol{ \Sigma}_2\|_{ F}^2$

Where, $\| \boldsymbol{ \Sigma}_2\|_{ F}^2=\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F}^2$, completing the importing of this theorem.

$\blacksquare \, \,$

Now, we proceed to the analysis of the variance of Hutch++:

## Theorem 2: Variance of Hutch++ Theorem

Fix parameters $k$ and $\ell$. Let $p=2k+1$ and construct $\boldsymbol{ Q}\in\mathbb R^{n \times k+p}$ as in Theorem 1. Then let

$\text{Hutch++}(\boldsymbol{ A}) \;{\vcentcolon=}\; \text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q}) + H_\ell((\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})$

Then,

$\mathbb E[\text{Hutch++}(\boldsymbol{ A})] = \text{tr}(\boldsymbol{ A})$ $\text{Var}[\text{Hutch++}(\boldsymbol{ A})] \leq \frac{1}{k\ell}\text{tr}^2(\boldsymbol{ A})$

Proof.

We first look at the expectation by using the Tower Law, the expectation of Hutchinson's estimator, the cyclic property of trace, and the linearity of trace: \begin{aligned} \mathbb E[\text{Hutch++}(\boldsymbol{ A})] &= \mathbb E[\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q}) + H_\ell((\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})] \\ &= \mathbb E[\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q})] + \mathbb E\left[\mathbb E[H_\ell((\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})|\boldsymbol{ Q}]\right] \\ &= \mathbb E[\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q})] + \mathbb E[\text{tr}((\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})] \\ &= \mathbb E[\text{tr}(\boldsymbol{ A}\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)] + \mathbb E[\text{tr}(\boldsymbol{ A}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal\boldsymbol{ A})] \\ &= \mathbb E[\text{tr}(\boldsymbol{ A}\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)+\text{tr}(\boldsymbol{ A}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal\boldsymbol{ A})] \\ &= \text{tr}(\boldsymbol{ A}) \end{aligned} We now move onto the Variance bound. First, recall the Conditional Variance Formula, which says for any random variables $X,Y$,

$\text{Var}[Y] = \mathbb E[\text{Var}[Y|X]] + \text{Var}[\mathbb E[Y|X]]$

Taking $Y=\text{Hutch++}(\boldsymbol{ A})$ and $X=\boldsymbol{ Q}$, we can bound the variance of Hutch++:

$\text{Var}[\text{Hutch++}] = \mathbb E[\text{Var}[\text{Hutch++} | \boldsymbol{ Q}]] + \text{Var}[\mathbb E[\text{Hutch++} | \boldsymbol{ Q}]]$

The second term above is always zero: \begin{aligned} \text{Var}[\mathbb E[\text{Hutch++}(\boldsymbol{ A}) | \boldsymbol{ Q}]] &= \text{Var}\left[\mathbb E[\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q}) + H_\ell( (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A}) | \boldsymbol{ Q}]\right] \\ &= \text{Var}\left[\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q}) + \text{tr}((\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A}) \right] \\ &= \text{Var}\left[\text{tr}(\boldsymbol{ A})\right] \\ &= 0 \end{aligned} So, we just need to bound the first term from the conditional variance formula: \begin{aligned} \mathbb E[\text{Var}[\text{Hutch++} | \boldsymbol{ Q}]] &= \mathbb E[\text{Var}[\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q})+H_\ell( (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})|Q]] \\ &= \mathbb E[\text{Var}[H_\ell( (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})|Q]] \\ &\leq {\textstyle\frac{ 2}{ \ell}} \mathbb E[\| (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A}\|_{ F}^2] \\ &\leq {\textstyle\frac{ 2}{ \ell}} (1+\frac{k}{p-1})\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F}^2 \\ \end{aligned} Where the first inequality uses the variance of Hutchinson's Estimator, and the second uses Theorem 1. Recall that we set $p=k+1$. After substituting and simplifying the expression, we find

$\mathbb E[\text{Var}[\text{Hutch++} | \boldsymbol{ Q}]] \leq {\textstyle\frac{ 4}{ \ell}}\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F}^2$

And by Lemma 3, we complete the proof:

${\textstyle\frac{ 4}{ \ell}}\| \boldsymbol{ A}-\boldsymbol{ A}_k\|_{ F}^2 \leq {\textstyle\frac{ 4}{ \ell}} \cdot {\textstyle\frac{ 1}{ 4\ell}}\text{tr}^2(\boldsymbol{ A}) = {\textstyle\frac{ 1}{ k\ell}}\text{tr}^2(\boldsymbol{ A})$
$\blacksquare \, \,$

## Part 4: Practitioner Advice

The variance in Theorem 2 still has two parameters in it, which leave a bit of ambiguity for practitioners. So, we quickly analyze the choice of $k$ and $\ell$ that minimizes our bound on the variance.

Formally, suppose we are allowed to compute exactly $m$ matrix-vector products with $\boldsymbol{ A}$. Then,

• Computing $\boldsymbol{ Q}$ uses $k+p=2k+1$ products

• Computing $\text{tr}(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q})$ uses $k+p=2k+1$ products

• Computing $H_\ell( (\boldsymbol{ I}-\boldsymbol{ Q}\boldsymbol{ Q}^\intercal)\boldsymbol{ A})$ uses $\ell$ products

So, we have $m=2(2k+1)+\ell=4k+2+\ell$. We can then verify that $k=\frac{m-2}{8}$ and $\ell=\frac{m}{2}-1$ minimizes the variance. Notably, this is equivalent to setting $\ell=4k$, which looks more intuitive. This produces a final variance of $\text{Var}[\text{Hutch++}(\boldsymbol{ A})]\leq\frac{16}{(m-2)^2}\text{tr}^2(\boldsymbol{ A})$.

# Nyström-Hutch++

One additional approach to implementing Hutch++ for PSD matrices considers the Nyström method, where we use the approximation

$\tilde\boldsymbol{ A}_k = (\boldsymbol{ A}\boldsymbol{ Q})(\boldsymbol{ Q}^\intercal\boldsymbol{ A}\boldsymbol{ Q})^{-1}(\boldsymbol{ A}\boldsymbol{ Q})^{\intercal}$

Note that we can compute $\tilde\boldsymbol{ A}_k=\boldsymbol{ Y}(\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}\boldsymbol{ Y}^\intercal$ where $\boldsymbol{ Y}=\boldsymbol{ A}\boldsymbol{ Q}$ contains all the matrix-vector products with $\boldsymbol{ A}$. We can compute $\text{tr}(\tilde\boldsymbol{ A}_k)$ efficiently:

\begin{aligned} \text{tr}(\tilde\boldsymbol{ A}_k) &= \text{tr}(\boldsymbol{ Y}(\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}\boldsymbol{ Y}^\intercal) \\ &= \text{tr}((\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}(\boldsymbol{ Y}^\intercal\boldsymbol{ Y})) \end{aligned}

The matrix inversion is only for a $(k+p) \times (k+p)$ sized matrix, which is easy. We can then define the NYS-Hutch++ estimator as

\begin{aligned} \text{NYS-Hutch++}(\boldsymbol{ A}) &\;{\vcentcolon=}\; \text{tr}((\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}(\boldsymbol{ Y}^\intercal\boldsymbol{ Y})) + H_\ell(\boldsymbol{ A}-\tilde\boldsymbol{ A}_k) \\ \end{aligned}

We can write out this Hutchinson's Estimator step as

$H_\ell(\boldsymbol{ A}-\tilde\boldsymbol{ A}_k) = \frac{1}{\ell}\left(\text{tr}(\boldsymbol{ G}^\intercal\boldsymbol{ A}\boldsymbol{ G}) - \text{tr}(\boldsymbol{ R}^\intercal(\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}\boldsymbol{ R}) \right)$

where $\boldsymbol{ G}$ has random sign-bit entries and $\boldsymbol{ R}\;{\vcentcolon=}\;\boldsymbol{ Y}^\intercal\boldsymbol{ G}$. Formally, we can write this out as the following algorithm:

Algorithm 1: NYS-Hutch++

---

input: Matrix-Vector Oracle access to $\boldsymbol{ A}\in\mathbb R^{d \times d}$. Number $m$ of queries.

output: Approximation to $\text{tr}(\boldsymbol{ A})$

1. Sample $\boldsymbol{ S}\in\mathbb R^{d \times \frac{m}{4}}$ and $\boldsymbol{ G}\in\mathbb R^{d \times \frac{m}{2}}$ with i.i.d. $\{+1,-1\}$ entries.

2. Compute an orthonormal basis $\boldsymbol{ Q}\in\mathbb R^{d \times \frac{m}{4}}$ for the span of $\boldsymbol{ A}\boldsymbol{ S}$

3. Compute $\boldsymbol{ Y}=\boldsymbol{ A}\boldsymbol{ Q}$ and $\boldsymbol{ R}=\boldsymbol{ Y}^\intercal\boldsymbol{ G}$

4. return $\text{tr}((\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}(\boldsymbol{ Y}^\intercal\boldsymbol{ Y})) + \frac{2}{m}\left(\text{tr}(\boldsymbol{ G}^\intercal\boldsymbol{ A}\boldsymbol{ G}) - \text{tr}(\boldsymbol{ R}^\intercal(\boldsymbol{ Q}^\intercal\boldsymbol{ Y})^{-1}\boldsymbol{ R}) \right)$

Or, in matlab:

function trace_est=simple_nystrom_hutchplusplus(A, num_queries)
S = 2*randi(2,size(A,1),ceil(num_queries/4))-3;
G = 2*randi(2,size(A,1),floor(num_queries/2))-3;
[Q,~] = qr(A*S,0);
Y = A*Q;
R = Y' * G;
QYinv = pinv(Q' * Y);
trace_est = trace(QYinv * (Y'*Y)) + 1/size(G,2)*[trace(G'*A*G) - trace(R'*QYinv*R)];
end

# Epilogue

## Quiz Questions

If you have read Hutch++, then quiz yourself by trying to answer some intuitive high-level questions:

1. Why is the constant-factor approximation on our low-rank approximation $\|A-\tilde{A}_k\|_F \leq 2 \|A-A_k\|_F$ sufficient? Why don't we need a relative error approximation like $\|A-\tilde{A}_k\|_F \leq (1+\varepsilon) \|A-A_k\|_F$?

Give an intuitive answer, not a mathematical one.

## Footnotes

 We analyze variance in this page for convenience, but all of these confidence intervals hold with high probability (i.e. with logarithmic dependence on the failure probability), and this is widely analyzed in the formal publications.