where $\mathrm{Var}_{{\boldsymbol{\xi}}}$ denotes the conditional variance
with respect to ${\boldsymbol{\xi}}$. Recall that
$\xi_{i}=\sqrt{\frac{b_{i}}{dp_{i}}}$, where
$b_{i}=\sum_{t=1}^{d}1_{[s_{t}=i]}$, with $s_{t}$ sampled i.i.d. from
$(p_{1},...,p_{n})$ and $p_{i}\approx_{O(1)}\|\mathbf{u}_{i}\|^{2}/d$ (here,
$\mathbf{u}_{i}^{\scriptscriptstyle{\top}}$ denotes the $i$th row of
$\mathbf{U}$). Thus,
$\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}}\mathbf{U}_{{\boldsymbol{\xi}}}=\sum_{t=1}^{d}\frac{\mathbf{u}_{s_{t}}\mathbf{u}_{s_{t}}^{\scriptscriptstyle{\top}}}{dp_{s_{t}}}$
and it follows that:

$\displaystyle\mathrm{Var}\big{[}\mathrm{tr}(\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}})\big{]}$
$\displaystyle=\mathrm{Var}\bigg{[}\sum_{t=1}^{d}\frac{\mathbf{u}_{s_{t}}^{\scriptscriptstyle{\top}}\mathbf{C}\mathbf{u}_{s_{t}}}{dp_{s_{t}}}\bigg{]}=d\,\mathrm{Var}\bigg{[}\frac{\mathbf{u}_{s_{1}}^{\scriptscriptstyle{\top}}\mathbf{C}\mathbf{u}_{s_{1}}}{dp_{s_{1}}}\bigg{]}$
$\displaystyle\leq
d\,\mathbb{E}\bigg{[}\frac{\mathrm{tr}(\mathbf{C}\mathbf{u}_{s_{1}}\mathbf{u}_{s_{1}}^{\scriptscriptstyle{\top}}\mathbf{C}\mathbf{u}_{s_{1}}\mathbf{u}_{s_{1}}^{\scriptscriptstyle{\top}})}{d^{2}p_{s_{1}}^{2}}\bigg{]}\leq\mathbb{E}\bigg{[}\frac{\|\mathbf{u}_{s_{1}}\|^{2}}{dp_{s_{1}}}\,\frac{\mathbf{u}_{s_{1}}^{\scriptscriptstyle{\top}}\mathbf{C}^{2}\mathbf{u}_{s_{1}}}{p_{s_{1}}}\bigg{]}$
$\displaystyle\leq
O(1)\,\mathbb{E}\bigg{[}\frac{\mathbf{u}_{s_{1}}^{\scriptscriptstyle{\top}}\mathbf{C}^{2}\mathbf{u}_{s_{1}}}{p_{s_{1}}}\bigg{]}=O(1)\,\mathrm{tr}(\mathbf{U}\mathbf{C}^{2}\mathbf{U}^{\scriptscriptstyle{\top}})=O(1)\,\mathrm{tr}(\mathbf{B}^{2}).$

The Bai-Silverstein inequality (Lemma 13) implies that
$\mathrm{Var}_{{\boldsymbol{\xi}}}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}}\mathbf{x}]\leq
O(1)\cdot\mathrm{tr}\big{(}(\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}})^{2}\big{)}$,
so we have:

$\displaystyle\mathbb{E}\big{[}\mathrm{Var}_{{\boldsymbol{\xi}}}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}}\mathbf{x}]\big{]}$
$\displaystyle\leq
O(1)\cdot\mathbb{E}\Big{[}\mathrm{tr}\big{(}(\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}})^{2}\big{)}\Big{]}=O(1)\cdot\mathbb{E}\bigg{[}\mathrm{tr}\bigg{(}\Big{(}\sum_{t=1}^{d}\frac{\mathbf{C}\mathbf{u}_{s_{t}}\mathbf{u}_{s_{t}}^{\scriptscriptstyle{\top}}}{dp_{s_{t}}}\Big{)}^{2}\bigg{)}\bigg{]}$
$\displaystyle\leq
O(1)\sum_{t=1}^{d}\mathbb{E}\bigg{[}\frac{\mathrm{tr}(\mathbf{C}\mathbf{u}_{s_{t}}\mathbf{u}_{s_{t}}^{\scriptscriptstyle{\top}}\mathbf{C}\mathbf{u}_{s_{t}}\mathbf{u}_{s_{t}}^{\scriptscriptstyle{\top}})}{d^{2}p_{s_{t}}^{2}}\bigg{]}+O(1)\sum_{t\neq
r}\mathbb{E}\bigg{[}\frac{\mathrm{tr}(\mathbf{C}\mathbf{u}_{s_{t}}\mathbf{u}_{s_{t}}^{\scriptscriptstyle{\top}}\mathbf{C}\mathbf{u}_{s_{r}}\mathbf{u}_{s_{r}}^{\scriptscriptstyle{\top}})}{dp_{s_{t}}\cdot
dp_{s_{r}}}\bigg{]}$ $\displaystyle\leq
O(1)\,\mathrm{tr}(\mathbf{B}^{2})+O(1)\,\mathrm{tr}\Big{(}\mathbf{C}\,\mathbb{E}\Big{[}\frac{\mathbf{u}_{s_{1}}\mathbf{u}_{s_{1}}^{\scriptscriptstyle{\top}}}{p_{s_{1}}}\Big{]}\mathbf{C}\,\mathbb{E}\Big{[}\frac{\mathbf{u}_{s_{2}}\mathbf{u}_{s_{2}}^{\scriptscriptstyle{\top}}}{p_{s_{2}}}\Big{]}\Big{)}\leq
O(1)\cdot\mathrm{tr}(\mathbf{B}^{2}).$

Thus, we obtain the desired bound:

$\displaystyle\mathrm{Var}\big{[}(\mathbf{x}\circ{\boldsymbol{\xi}})^{\scriptscriptstyle{\top}}\mathbf{B}(\mathbf{x}\circ{\boldsymbol{\xi}})\big{]}\leq
O(1)\cdot\mathrm{tr}(\mathbf{B}^{2}),$

which completes the proof. ∎

###  C.3 Lower bounds for other sketching methods

In this section, we show lower bounds for Condition 2 (Restricted Bai-
Silverstein) in the context of existing fast sketching techniques. To do that,
we first discuss the basic requirement of the framework defined by Theorem 11,
namely that the sketching matrix $\mathbf{S}$ must have i.i.d. rows.

Fast sketches with i.i.d. rows. In our discussion, we will focus on three
types of sketches: approximate leverage score sampling [DMM06], Subsampled
Randomized Hadamard Transform [AC09], and sparse embedding matrices
(extensions of the CountSketch [CW17], see [NN13, Coh16]), all of which can be
implemented in time nearly linear in the input size. The i.i.d. row assumption
can be easily satisfied by any row sampling sketch, including approximate
leverage score sampling. The SRHT technically does not satisfy this
assumption, however if we treat the Randomized Hadamard Transform as a
preprocessing step (given that it does not distort the covariance matrix),
then the subsampling part can be analyzed analogously as leverage score
sampling. In the case of sparse embedding matrices, the most commonly studied
variant has a fixed number of non-zeros per column of $\mathbf{S}$ and so it
does not have independent rows, however, it is known that a variant with
independently sparsified entries (which fits into the setup of Theorem 11)
achieves nearly matching approximation guarantees [Coh16].

Leverage score sampling. Let $\mathbf{S}$ be a row sampling sketch of size
$m$, i.e., each row is distributed independently as
$\frac{1}{\sqrt{m}}\mathbf{x}^{\scriptscriptstyle{\top}}$, where
$\mathbf{x}=\frac{1}{\sqrt{p_{s}}}\mathbf{e}_{s}$ and $s$ is an index drawn
from distribution $(p_{1},...,p_{n})$. Given a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$ of rank $d$, we call this an approximate
leverage score sampling sketch if $p_{i}\approx_{O(1)}l_{i}/d$ for all $i$,
where
$l_{i}=\mathbf{a}_{i}^{\scriptscriptstyle{\top}}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{a}_{i}$
is the $i$th leverage score of $\mathbf{A}$. We will present two lower bound
constructions.

1\. Approximate sampling and arbitrary $\mathbf{A}$. Now, suppose that $n$ is
even and consider the following specific example:

$\displaystyle p_{j}=\begin{cases}l_{j}/2d,&\text{for }j\leq
n/2,\\\\[2.84526pt] 3l_{j}/2d,&\text{otherwise.}\end{cases}$

Further, consider the matrix
$\mathbf{B}=\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}=\mathbf{P}$,
which is the projection onto the column-span of $\mathbf{A}$, and therefore
satisfies the restriction requirement in the Restricted Bai-Silverstein
condition. Then, since
$\mathrm{tr}(\mathbf{B}^{2})=\mathrm{tr}(\mathbf{P}^{2})=\mathrm{tr}(\mathbf{P})=d$,
we have:

$\displaystyle\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}]$
$\displaystyle=\mathbb{E}\Big{[}\big{(}\mathbf{e}_{s}^{\scriptscriptstyle{\top}}\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{e}_{s}/p_{s}-d\big{)}^{2}\Big{]}=\mathbb{E}\Big{[}\big{(}l_{s}/p_{s}-d\big{)}^{2}\Big{]}\geq(d/3)^{2}$

2\. Exact sampling and a specific $\mathbf{A}$. Suppose that
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}=\mathbf{I}$, each
$\mathbf{a}_{i}$ is a standard basis vector scaled by $\sqrt{d/n}$ and we are
sampling index $s$ according to exact leverage scores, i.e., uniformly at
random. Then, letting $\mathbf{x}=\frac{1}{\sqrt{p_{s}}}\mathbf{e}_{s}$ and
$\mathbf{B}=\mathbf{A}\mathbf{C}\mathbf{A}^{\scriptscriptstyle{\top}}$, we
have:

$\displaystyle\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}]$
$\displaystyle=\mathbb{E}\Big{[}\big{(}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}-\mathrm{tr}(\mathbf{B})\big{)}^{2}\Big{]}=\mathbb{E}\Big{[}\big{(}d\cdot\frac{\mathbf{a}_{s}^{\scriptscriptstyle{\top}}\mathbf{C}\mathbf{a}_{s}}{\mathbf{a}_{s}^{\scriptscriptstyle{\top}}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{a}_{s}}-\mathrm{tr}(\mathbf{C})\big{)}^{2}\Big{]}$
$\displaystyle=d^{2}\cdot\frac{1}{d}\sum_{j=1}^{d}\Big{(}C_{jj}-\frac{1}{d}\sum_{i=1}^{d}C_{ii}\Big{)}^{2}=d^{2}\cdot\Omega(1),\quad\text{if}\quad
C_{ii}=\begin{cases}1/2,&\text{for even $i$,}\\\ 3/2,&\text{for odd
$i$.}\end{cases}$

In both constructions, we have $\mathrm{tr}(\mathbf{B}^{2})=\Theta(d)$, so
this implies that for leverage score sampling, Condition 2 can only be shown
with factor $\alpha=\Omega(d)$, as opposed to $O(1)$ for sub-gaussian sketches
and LESS embeddings.

Data-oblivious sparse embeddings. Let $\mathbf{S}$ be a sketch of size $m$,
where each row is distributed independently as
$\frac{1}{\sqrt{m}}\mathbf{x}^{\scriptscriptstyle{\top}}$ and
$\mathbf{x}=(\sqrt{\frac{m}{s}}b_{1}r_{1},...,\sqrt{\frac{m}{s}}b_{n}r_{n})$,
with $b_{i}\sim\mathrm{Bernoulli}(\frac{s}{m})$ and $r_{i}$ distributed as a
uniformly random sign. While this is not the most commonly studied variant of
a sparse embedding, it is known to satisfy the subspace embedding property for
sketch size $m=O(d\log d)$ with sparsity level $s=O(\log d)$ [Coh16], which
matches the state-of-the-art for sparse embeddings. Other sparse embeddings
have non-i.i.d. row distributions [CW17, NN13, MM13], and so they do not fit
into the framework laid out by Theorem 11. The key difference between the
sparsification of $\mathbf{x}$ relative to our LESS embeddings is that it is
data-oblivious. We can exploit that in our lower bound example by choosing an
extremely skewed leverage score distribution of matrix $\mathbf{A}$. In
particular, suppose that
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}=\mathbf{I}$ and moreover,
$\mathbf{a}_{i}=\mathbf{e}_{i}$ for $i=1,...,k$ (where $1\leq k\leq d$) and
for all $i>k$, the first $k$ coordinates of $a_{i}$ are zero. This
construction ensures that the first $k$ leverage scores of $\mathbf{A}$ are
equal $1$. Once again setting
$\mathbf{B}=\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}$,
we get:

$\displaystyle\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}]\geq\sum_{i=1}^{k}\mathrm{Var}\Big{[}\frac{m}{s}b_{i}r_{i}\Big{]}=k\cdot\frac{m}{s}\Big{(}1-\frac{s}{m}\Big{)}.$

If we let $k=\Omega(d)$, then we get
$\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}]\geq\Omega(m/s)\cdot\mathrm{tr}(\mathbf{B}^{2})$.
Thus, unless we zero-out merely a constant fraction of entries of
$\mathbf{S}$, the sketching matrix will not satisfy Condition 2 with a
constant factor $\alpha=O(1)$. We conjecture that this example can be extended
to show a general lower bound on the inversion bias, as we did for approximate
leverage score sampling.

##  Appendix D Subspace embedding guarantee for LESS embeddings

In this section, we prove Lemma 12 and Theorem 8. In particular, we prove that
LESS embeddings achieve the subspace embedding property for a sketch of size
$O(d\log d)$ (Lemma 12), thus establishing Condition 1. Then, at the end of
the section we briefly discuss how to combine Lemmas 12 and 28, using the
structural conditions via Theorem 11, to obtain Theorem 8.

###  D.1 Proof of Lemma 12

First, note that instead of directly showing the subspace embedding of
$\mathbf{S}\mathbf{A}$ for the span of $\mathbf{A}$, it suffices to show the
guarantee when replacing $\mathbf{A}$ with its orthonormal basis matrix
$\mathbf{U}=\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1/2}$,
since
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}=(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{\frac{1}{2}}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{\frac{1}{2}}$.
Then, a standard technique, e.g., as used for leverage score sampling
sketches, relies on the following decomposition of
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}$
as an average of independent rank-one p.s.d. random matrices:

$\displaystyle\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}=\sum_{i=1}^{m}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{U},$

where $\mathbf{s}_{i}^{\scriptscriptstyle{\top}}$ represents the $i$th row of
$\mathbf{S}$. For standard leverage score sampling sketches it suffices to use
the matrix Chernoff bound [Tro12, Theorem 1.1], which uses an almost sure
bound on each rank-one matrix to ensure concentration around the mean,
$\mathbb{E}[\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}]=\mathbf{I}$.
However, in the case of a leverage score sparsified embedding an almost sure
bound is not sufficient. Instead, we show that the rank-one matrices
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}$
exhibit sub-exponential tails on all of their moments, as required by the
following variant of the matrix Bernstein bound.

######  Lemma 29 ([Tro12, Theorem 6.2]).

For $i=1,2,...$, consider a finite sequence $\mathbf{X}_{i}$ of $d\times d$
independent and symmetric random matrices such that

$\mathbb{E}[\mathbf{X}_{i}]=\mathbf{0},\quad\mathbb{E}[\mathbf{X}_{i}^{p}]\preceq\frac{p!}{2}\cdot
R^{p-2}\mathbf{A}_{i}^{2}\quad\textmd{for}\quad p=2,3,...$

Then, defining the variance parameter
$\sigma^{2}=\|\sum_{i}\mathbf{A}_{i}^{2}\|$, for any $t>0$ we have:

$\displaystyle\mathrm{Pr}\bigg{\\{}\lambda_{\max}\Big{(}\sum\nolimits_{i}\mathbf{X}_{i}\Big{)}\geq
t\bigg{\\}}$ $\displaystyle\leq
d\cdot\exp\left(\frac{-t^{2}/2}{\sigma^{2}+Rt}\right).$

We apply the above result for
$\mathbf{X}_{i}=\pm(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}-\frac{1}{m}\mathbf{I})$,
where
$\mathbf{s}_{i}=\frac{1}{\sqrt{m}}(\mathbf{x}_{i}\circ{\boldsymbol{\xi}})$ is
a leverage score sparsified sub-gaussian random vector. We next establish the
subexponential moment bound needed for the matrix Bernstein bound.

######  Lemma 30.

Fix a matrix $\mathbf{U}\in\mathbb{R}^{n\times d}$ such that
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{I}$. Suppose that
${\boldsymbol{\xi}}$ is a leverage score sparsifier for $\mathbf{U}$ and
$\mathbf{x}$ has i.i.d. $O(1)$-sub-gaussian entries with mean zero and unit
variance. Then, there is $C=O(1)$ such that for all $p=2,3,...$ we have

$\displaystyle\Big{\|}\mathbb{E}\Big{[}\Big{(}\mathbf{U}^{\scriptscriptstyle{\top}}(\mathbf{x}\circ{\boldsymbol{\xi}})(\mathbf{x}\circ{\boldsymbol{\xi}})^{\scriptscriptstyle{\top}}\mathbf{U}-\mathbf{I}\Big{)}^{p}\Big{]}\Big{\|}\leq\frac{p!}{2}\cdot(Cd)^{p-1}.$

Now, the matrix Bernstein bound (Lemma 29) can be invoked with
$\mathbf{A}_{i}^{2}=\frac{Cd}{m^{2}}\cdot\mathbf{I}$ and
$\sigma^{2}=R=\frac{Cd}{m}$, obtaining that for $\eta\in(0,1)$:

$\displaystyle\mathrm{Pr}\Big{\\{}\big{\|}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}-\mathbf{I}\big{\|}\geq\eta\Big{\\}}\leq
2d\cdot\exp\Big{(}-\frac{\eta^{2}m}{4Cd}\Big{)}\leq\delta\quad\text{for}\quad
m\geq 4Cd\log(2d/\delta)/\eta^{2},$

which completes the proof.

###  D.2 Proof of Lemma 30

The key part of our proof of Lemma 30 involves establishing the following
concentration inequality which can be viewed as a form of the Hanson-Wright
inequality [RV13] that takes advantage of the leverage score sparsifier
${\boldsymbol{\xi}}$, similarly as we did for the Restricted Bai-Silverstein
inequality (Lemma 28).

######  Lemma 31.

Fix a matrix $\mathbf{U}\in\mathbb{R}^{n\times d}$ such that
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{I}$. Suppose that
${\boldsymbol{\xi}}$ is a leverage score sparsifier for $\mathbf{U}$ and
$\mathbf{x}$ has independent $O(1)$-sub-gaussian entries with mean zero and
unit variance. Then, there is $c=\Omega(1)$ and $C=O(1)$ such that for any
$t\geq Cd$ we have:

$\displaystyle\mathrm{Pr}\Big{\\{}(\mathbf{x}\circ{\boldsymbol{\xi}})^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}(\mathbf{x}\circ{\boldsymbol{\xi}})\geq
t\Big{\\}}\leq\exp\Big{(}-c\,\big{(}\sqrt{t}+t/d\big{)}\Big{)}.$

###### Proof.

We use the shorthand
$\mathbf{U}_{\boldsymbol{\xi}}=\operatorname*{\mathop{\mathrm{diag}}}({\boldsymbol{\xi}})\mathbf{U}$.
Similarly as for Lemma 28, our strategy is to show that the sparsification
$\mathbf{U}_{\boldsymbol{\xi}}$ preserves enough of the structure of
$\mathbf{U}$ so that we can apply the classical Hanson-Wright inequality,
which is repeated below, following [RV13],

######  Lemma 32 ([RV13, Theorem 1.1]).

Let $\mathbf{x}$ have independent $O(1)$-sub-gaussian entries with mean zero
and unit variance. Then, there is $c=\Omega(1)$ such that for any $n\times n$
matrix $\mathbf{B}$ and $t\geq 0$,

$\displaystyle\mathrm{Pr}\Big{\\{}|\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}-\mathrm{tr}(\mathbf{B})|\geq
t\Big{\\}}\leq
2\exp\bigg{(}-c\min\Big{\\{}\frac{t^{2}}{\|\mathbf{B}\|_{F}^{2}},\frac{t}{\|\mathbf{B}\|}\Big{\\}}\bigg{)}.$

To show that the leverage score sparsification $\mathbf{U}_{\boldsymbol{\xi}}$
is sufficiently accurate, we can rely on the matrix Chernoff bound, repeated
below, and the following decomposition:

$\displaystyle\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}=\sum_{i=1}^{d}\frac{\mathbf{u}_{s_{i}}\mathbf{u}_{s_{i}}^{\scriptscriptstyle{\top}}}{dp_{s_{i}}},$

where $s_{i}$ are the random indices sampled from the approximate leverage
score distribution $(p_{1},...,p_{n})$ (see Definition 6). For simplicity, we
only repeat the large deviation part of the Chernoff bound, which is the one
relevant to our analysis.

######  Lemma 33 ([Tro12, Theorem 1.1 and Remark 5.3]).

For $i=1,2,...$, consider a finite sequence $\mathbf{X}_{i}$ of $d\times d$
independent positive semi-definite random matrices such that
$\mathbb{E}\big{[}\sum_{i}\mathbf{X}_{i}\big{]}=\mathbf{I}$ and
$\|\mathbf{X}_{i}\|\leq R$. Then, for any $t\geq\mathrm{e}$, we have:

$\displaystyle\mathrm{Pr}\Big{\\{}\Big{\|}\sum_{i}\mathbf{X}_{i}\Big{\|}\geq
t\Big{\\}}\leq d\cdot\Big{(}\frac{\mathrm{e}}{t}\Big{)}^{t/R}.$

We apply the matrix Chernoff to
$\mathbf{X}_{i}=\frac{1}{dp_{s_{i}}}\mathbf{u}_{s_{i}}\mathbf{u}_{s_{i}}^{\scriptscriptstyle{\top}}$,
noting that since $p_{i}\geq\|\mathbf{u}_{i}\|^{2}/Rd$ for $R=O(1)$, it
follows that $\|\mathbf{X}_{i}\|\leq R$. Moreover,
$\mathbb{E}[\sum_{i=1}^{d}\mathbf{X}_{i}]=\mathbf{I}$, so for $t\geq O(1)\cdot
d$ we have:

$\displaystyle\mathrm{Pr}\big{\\{}\|\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\|\geq\sqrt{t}\big{\\}}\leq
d\exp\big{(}-\sqrt{t}\ln(\sqrt{t}/\mathrm{e})/R\big{)}\leq\exp(-c\sqrt{t}),$

for some $c=\Omega(1)$. Also, note that
$\|\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\|\leq\mathrm{tr}(\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}})\leq
Rd$ almost surely, which implies that event
$\mathcal{E}:\big{[}\|\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\|\leq\min\\{\sqrt{t},Rd\\}\big{]}$
holds with probability at least $1-\exp(-c(\sqrt{t}+t/d))$. Conditioned on
$\mathcal{E}$, it holds that
$\|\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\|_{F}^{2}\leq\|\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\|\cdot\mathrm{tr}(\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}})\leq\min\\{\sqrt{t},Rd\\}\cdot
Rd$, so applying Lemma 32 for fixed ${\boldsymbol{\xi}}$ we get:

$\displaystyle\mathrm{Pr}\big{\\{}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\geq
Rd+t\ \mid\ {\boldsymbol{\xi}},\mathcal{E}\big{\\}}$ $\displaystyle\leq
2\exp\bigg{(}-c\min\Big{\\{}\frac{t^{2}}{\|\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\|_{F}^{2}},\frac{t}{\|\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}\|}\Big{\\}}\bigg{)}$
$\displaystyle\leq
2\exp\bigg{(}-c\min\Big{\\{}\frac{t^{2}}{\min\\{\sqrt{t},Rd\\}\cdot{Rd}},\frac{t}{\min\\{\sqrt{t},Rd\\}}\Big{\\}}\bigg{)}$
$\displaystyle\leq 2\exp\big{(}-c(\sqrt{t}+t/Rd)\big{)}.$

Appropriately rescaling $t$, we obtain the claim. ∎

We are now ready to present the proof of Lemma 30, obtaining subexponential
moment bounds for the random matrix
$\mathbf{U}^{\scriptscriptstyle{\top}}(\mathbf{x}\circ{\boldsymbol{\xi}})(\mathbf{x}\circ{\boldsymbol{\xi}})^{\scriptscriptstyle{\top}}\mathbf{U}$,
thus completing the proof of the subspace embedding guarantee for leverage
score sparsified sketches.

###### Proof of Lemma 30.

Throughout the proof, we will use the shorthand
$\mathbf{U}_{{\boldsymbol{\xi}}}=\operatorname*{\mathop{\mathrm{diag}}}({\boldsymbol{\xi}})\mathbf{U}$.
It is easy to show by induction over $p$ that:

$\displaystyle\underbrace{\big{(}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{{\boldsymbol{\xi}}}-\mathbf{I}\big{)}^{p}}_{\mathbf{Z}^{p}}=\big{(}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}-1\big{)}^{p-1}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}-\underbrace{\big{(}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{{\boldsymbol{\xi}}}-\mathbf{I}\big{)}^{p-1}}_{\mathbf{Z}^{p-1}}.$

Thus, it follows that for any $p=2,3,...$ (both even and odd) we have the
following upper bound:

$\displaystyle\big{\|}\mathbb{E}[\mathbf{Z}^{p}]\big{\|}\leq\Big{\|}\mathbb{E}\big{[}\underbrace{|\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}-1|^{p-1}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}}_{\mathbf{T}_{p}}\big{]}\Big{\|}+\big{\|}\mathbb{E}[\mathbf{Z}^{p-1}]\big{\|}.$

To bound the quadratic form
$\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}$
in the first term, we can use Lemma 31. In particular, the lemma implies that
the event
$\mathcal{E}:[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\leq
Cpd]$ fails with probability at most $\mathrm{e}^{-\sqrt{pd}}$ for a
sufficiently large $C=O(1)$, so we have:

$\displaystyle\big{\|}\mathbb{E}[\mathbf{T}_{p}]\big{\|}$
$\displaystyle\leq\big{\|}\mathbb{E}[\mathbf{T}_{p}\cdot\mathbf{1}_{\mathcal{E}}]\big{\|}+\big{\|}\mathbb{E}[\mathbf{T}_{p}\cdot\mathbf{1}_{\neg\mathcal{E}}]\big{\|}$
$\displaystyle\leq(pd)^{p-1}\big{\|}\mathbb{E}[\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}]\big{\|}+\mathbb{E}\big{[}(\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\cdot\mathbf{1}_{\neg\mathcal{E}})^{p}\big{]}$
$\displaystyle=(Cpd)^{p-1}+\int_{0}^{\infty}pt^{p-1}\mathrm{Pr}\big{\\{}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{\boldsymbol{\xi}}\mathbf{U}_{\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}\mathbf{x}\cdot\mathbf{1}_{\neg\mathcal{E}}>t\big{\\}}dt$
$\displaystyle\leq(Cpd)^{p-1}+p(Cpd)^{p}\mathrm{e}^{-\sqrt{pd}}+\int_{Cpd}^{\infty}pt^{p-1}\mathrm{e}^{-c(\sqrt{t}+t/d)}dt.$

Note that $(O(1)\,p)^{p+O(1)}d^{p-1}\leq
p^{p}(O(1)\,d)^{p-1}\leq(p!/2)(O(1)\,d)^{p-1}$, and also
$\mathrm{e}^{-\sqrt{pd}}\leq O(1/d)$, so the first two terms can be easily
bounded as desired. To bound the last term, we use the following integral
formula:

$\displaystyle\int t^{p-1}\mathrm{e}^{-\alpha
t^{\theta}}dt=-\frac{\Gamma(p/\theta,\alpha
t^{\theta})}{\theta\alpha^{p/\theta}}+\mathrm{const},$

which follows from the definition of the upper incomplete Gamma function
$\Gamma$. Note that for $p=2,3,...$ this function also satisfies:

$\displaystyle\Gamma(p,\lambda)$
$\displaystyle=(p-1)!\cdot\mathrm{Pr}\\{x<p\\}\qquad\text{for}\quad
x\sim\mathrm{Poisson}(\lambda),$
$\displaystyle\leq(p-1)!\cdot\mathrm{e}^{-c\lambda}\qquad\text{for}\quad\lambda\geq
2p,\ c=\Omega(1),$

where the last inequality is a standard tail bound for a Poisson random
variable. With a slight abuse of notation, we let $c$ denote the minimum of
the above constant $c$ and the constant $c$ from Lemma 31. We apply the
integral formula in two different ways, depending on $p$. First, if $p<d$ then
we have:

$\displaystyle\int_{Cpd}^{\infty}pt^{p-1}\mathrm{e}^{-c(\sqrt{t}+t/d)}dt$
$\displaystyle\leq\int_{Cpd}^{\infty}pt^{p-1}\mathrm{e}^{-c\sqrt{t}}dt=2pc^{-2p}\Gamma(2p,c\sqrt{Cpd})\leq
2c^{-2p}(2p)!\mathrm{e}^{-c^{2}\sqrt{Cpd}}.$

By using the fact that $\exp(-c^{2}\sqrt{Cpd})\leq\exp(-c^{2}p)=O(1/p)$, this
expression can be bounded by $(p!/2)(O(1)\,p)^{p-1}\leq(p!/2)(O(1)\,d)^{p-1}$.
Next, we consider the case when $p\geq d$. We have:

$\displaystyle\int_{Cpd}^{\infty}pt^{p-1}\mathrm{e}^{-c(\sqrt{t}+t/d)}dt$
$\displaystyle\leq\int_{Cpd}^{\infty}pt^{p-1}\mathrm{e}^{-ct/d}dt=p(d/c)^{p}\Gamma(p,cCp)\leq
p!d^{p}\mathrm{e}^{-c^{2}Cp},$

where the last inequality holds as long as $C\geq 2/c$. Here, we note that
$\mathrm{e}^{-c^{2}Cp}\leq O(1/d)$ since $p\geq d$, thus again obtaining a
bound of the form $(p!/2)(O(1)\,d)^{p-1}$. Putting everything together, we
conclude that:

$\displaystyle\|\mathbb{E}[\mathbf{Z}^{p}]\|\leq\frac{p!}{2}\,(O(1)\,d)^{p-1}+\|\mathbb{E}[\mathbf{Z}^{p-1}]\|.$

Recursively summing up this bound concludes the proof. ∎

###  D.3 Proof of Theorem 8

In Lemma 12, we showed that a LESS embedding of size $m\geq
4Cd\log(3d/\delta)$ satisfies Condition 1 (subspace embedding) for $\eta=1/2$
with probability $1-\delta/3$, as required by Theorem 11. Also, in Lemma 28 we
showed Condition 2 (Restricted Bai-Silverstein) with $\alpha=O(1)$ for a
leverage score sparsified sub-gaussian vector. Thus, as long as $\delta\leq
1/m^{3}$ and $m/3\geq 4Cd\log(3d/\delta)$, it follows that
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
is an $(\epsilon,\delta)$-unbiased estimator of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ for
$\epsilon=O(\sqrt{d}/m)$, and we obtain the desired guarantee. Note that the
condition for invoking Theorem 11 can be written as $m\geq
C^{\prime}d\log(m)$. This is satisfied for $m=C^{\prime}d\log(C^{\prime
2}d^{2})$, and since $m$ grows faster than $\log(m)$, it will also be
satisfied for all $m\geq C^{\prime}d\log(C^{\prime 2}d^{2})=O(d\log(d))$. This
completes the proof of Theorem 8.

##  Appendix E Averaging nearly-unbiased estimators

In this section, we show that averaging improves spectral approximation for
matrix estimators with small inversion bias, and as a consequence we prove
Corollaries 5 and 9 for averaging sketched inverse covariance matrix
estimators based on sub-gaussian sketches and LESS embeddings respectively.

###  E.1 Conditions for effective averaging of random matrices

We start with a more general result, which should be of interest to averaging
nearly-unbiased matrix estimators in settings other than inverse covariance
matrix estimation.

######  Lemma 34 (Conditions for effective averaging).

Suppose that $\delta\leq{\epsilon}\leq{\eta}\leq 1$ and
$\tilde{\mathbf{C}}_{1},...,\tilde{\mathbf{C}}_{q}$ are i.i.d. positive semi-
definite $d$-dimensional random matrices such that:

  1. 1.

$\tilde{\mathbf{C}}_{i}$ is an $({\epsilon},\delta/2q)$-unbiased estimator of
$\mathbf{C}$;

  2. 2.

$\tilde{\mathbf{C}}_{i}$ is an $({\eta},\delta/2q)$-approximation of
$\mathbf{C}$.

Then, $\frac{1}{q}\sum_{i=1}^{q}\tilde{\mathbf{C}}_{i}$ is an
$(\epsilon^{\prime},2\delta)$-approximation of $\mathbf{C}$ for
$\epsilon^{\prime}={\epsilon}+{\eta}\cdot
O\big{(}\sqrt{\frac{\ln(d/\delta)}{q}}\,\big{)}$.

###### Proof.

For this, we use a variant of the matrix Bernstein inequality given below.

######  Lemma 35 ([Tro12, Theorem 1.4]).

For $i=1,2,...$, consider a finite sequence $\mathbf{X}_{i}$ of $d\times d$
independent and symmetric random matrices such that

$\mathbb{E}[\mathbf{X}_{i}]=\mathbf{0},\quad\lambda_{\max}(\mathbf{X}_{i})\leq
R\quad\text{almost surely.}$

Then, defining the variance parameter
$\sigma^{2}=\|\sum_{i}\mathbb{E}[\mathbf{X}_{i}^{2}]\|$, for any $t>0$ we
have:

$\displaystyle\mathrm{Pr}\bigg{\\{}\lambda_{\max}\Big{(}\sum\nolimits_{i}\mathbf{X}_{i}\Big{)}\geq
t\bigg{\\}}$ $\displaystyle\leq
d\cdot\exp\left(\frac{-t^{2}/2}{\sigma^{2}+Rt/3}\right).$

Suppose that $\tilde{\mathbf{C}}$ is an $({\epsilon},\delta/2q)$-unbiased
estimator and an $({\eta},\delta/2q)$-spectral approximation for $\mathbf{C}$,
with $\mathcal{E}_{inv}$ and $\mathcal{E}_{sub}$ the associated high
probability events. For concreteness, let the $O(1)$ constant factor in
Definition 3 be denoted as $M$. Further, let
$\tilde{\mathbf{C}}_{i},\mathcal{E}_{inv}^{i},\mathcal{E}_{sub}^{i}$ be the
i.i.d. copies of $\tilde{\mathbf{C}}$ with their associated events. Finally,
let $\tilde{\mathbf{C}}_{i}^{\prime}$ be a random matrix obtained from
conditioning $\tilde{\mathbf{C}}_{i}$ on
$\mathcal{E}_{inv}^{i}\wedge\mathcal{E}_{sub}^{i}$, and coupled with
$\tilde{\mathbf{C}}_{i}$ so that
$\mathrm{Pr}(\tilde{\mathbf{C}}_{i}^{\prime}=\tilde{\mathbf{C}}_{i})\geq\mathrm{Pr}(\mathcal{E}_{inv}^{i}\wedge\mathcal{E}_{sub}^{i})\geq
1-\delta/q$ (this coupling can be obtained by considering a construction of
$\tilde{\mathbf{C}}_{i}^{\prime}$ via rejection sampling from
$\tilde{\mathbf{C}}_{i}$). We can bound the bias of
$\tilde{\mathbf{C}}_{i}^{\prime}$ (for any $i$) by observing that:

$\displaystyle-\delta/q\cdot\mathbb{E}[\tilde{\mathbf{C}}_{i}\mid\mathcal{E}_{inv}^{i},\neg\mathcal{E}_{sub}^{i}]\preceq\mathbb{E}[\tilde{\mathbf{C}}_{i}^{\prime}]-\mathbb{E}[\tilde{\mathbf{C}}_{i}\mid\mathcal{E}_{inv}^{i}]\preceq\frac{\delta/q}{1-\delta/q}\cdot\mathbb{E}[\tilde{\mathbf{C}}_{i}\mid\mathcal{E}_{inv}^{i}].$

Since we have
$\mathbb{E}[\tilde{\mathbf{C}}_{i}\mid\mathcal{E}_{inv}^{i}]\approx_{{\epsilon}}\mathbf{C}$
and
$\mathbb{E}[\tilde{\mathbf{C}}_{i}\mid\mathcal{E}_{inv}^{i},\neg\mathcal{E}_{sub}^{i}]\preceq
M\cdot\mathbf{C}$, it follows that
$\mathbb{E}[\tilde{\mathbf{C}}_{i}^{\prime}]$ is an
$\epsilon^{\prime}$-spectral approximation of $\mathbf{C}$ for
$\epsilon^{\prime}={\epsilon}+\frac{2\delta}{q}(1+{\epsilon}+M)$.

We will now apply the matrix Bernstein inequality (Lemma 35) to the sequence
of matrices:

$\mathbf{X}_{i}=\frac{1}{q}\Big{(}\mathbf{C}^{-\frac{1}{2}}\tilde{\mathbf{C}}_{i}^{\prime}\mathbf{C}^{-\frac{1}{2}}-\mathbb{E}\big{[}\mathbf{C}^{-\frac{1}{2}}\tilde{\mathbf{C}}_{i}^{\prime}\mathbf{C}^{-\frac{1}{2}}\big{]}\Big{)},\quad
i=1,...,q.$

Note that we have
$\tilde{\mathbf{C}}_{i}^{\prime}\approx_{{\eta}}\frac{1}{q}\mathbf{C}$, so it
follows that $\|\mathbf{X}_{i}\|\leq({\eta}+\epsilon^{\prime})/q$ and
$\sum_{i}\|\mathbf{X}_{i}^{2}\|\leq({\eta}+\epsilon^{\prime})^{2}/q$. Thus, we
conclude that for $t\in(0,1)$:

$\displaystyle\mathrm{Pr}\bigg{\\{}\Big{\|}\sum_{i=1}^{q}\mathbf{X}_{i}\Big{\|}\geq
t\,({\eta}+\epsilon^{\prime})\bigg{\\}}\leq 2d\exp\big{(}-t^{2}q/4\big{)}.$

Setting $t=\sqrt{4\ln(2d/\delta)/q}$ (without loss of generality, assume that
$t\leq 1$), we obtain that with probability $1-\delta$,

$\displaystyle\Big{\|}\frac{1}{q}\sum_{i=1}^{q}\mathbf{C}^{-\frac{1}{2}}\tilde{\mathbf{C}}_{i}^{\prime}\mathbf{C}^{-\frac{1}{2}}-\mathbf{I}\Big{\|}$
$\displaystyle\leq\Big{\|}\sum_{i=1}^{q}\mathbf{X}_{i}\Big{\|}+\Big{\|}\frac{1}{q}\sum_{i=1}^{q}\mathbf{C}^{-\frac{1}{2}}\mathbb{E}[\tilde{\mathbf{C}}_{i}^{\prime}]\mathbf{C}^{-\frac{1}{2}}-\mathbf{I}\Big{\|}$
$\displaystyle\leq
t\cdot({\eta}+\epsilon^{\prime})+\epsilon^{\prime}\leq{\epsilon}+{\eta}\cdot
O\Big{(}\sqrt{\tfrac{\log(d/\delta)}{q}}\,\Big{)}+O\Big{(}\frac{\delta
M}{q}\Big{)}.$

Note that under the assumptions that $M=O(1)$ and $\delta\leq{\eta}$, we can
absorb the last term into the middle term. Finally, observe that thanks to the
coupling and a union bound, the above bound holds with probability $1-2\delta$
after we replace $\tilde{\mathbf{C}}_{i}^{\prime}$ with
$\tilde{\mathbf{C}}_{i}$, completing the proof of Lemma 34. ∎

###  E.2 Proof of Corollary 5

Consider a sub-gaussian sketching matrix $\mathbf{S}$ of size $m\geq
C(d+\sqrt{d}/\epsilon+\log(2q/\delta))$. From Proposition 4, it follows that
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
is an $(\epsilon,\delta/2q)$-unbiased estimator of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$. Further, it is an
$({\eta},\delta/2q)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$, where
${\eta}=O(\sqrt{d/m})=O(\epsilon\cdot\sqrt{m})$. Thus, using Lemma 34, it
follows that for $q$ i.i.d. copies $\mathbf{S}_{1},...,\mathbf{S}_{q}$, the
averaged estimator
$\frac{1}{q}\sum_{i=1}^{q}(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}\mathbf{A})^{-1}$
is an $(\epsilon^{\prime\prime},2\delta)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ for

$\epsilon^{\prime\prime}=\epsilon+O\Big{(}\epsilon\cdot\sqrt{m\log(d/\delta)/q}\,\Big{)}.$

Setting $q=O(m\log(d/\delta))$ and adjusting the constants appropriately, we
obtain the claim.

###  E.3 Proof of Corollary 9

Consider a LESS embedding matrix $\mathbf{S}$ of size $m\geq
C(d\log(2dq/\delta)+\sqrt{d}/\epsilon)$. From Theorem 8, it follows that
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
is an $(\epsilon,\delta/2q)$-unbiased estimator of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$. Furthermore, the
theorem also implies that this matrix is an $({\eta},\delta/2q)$-approximation
of $(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ for
${\eta}=O(\sqrt{d\log(2dq/\delta)/m})=O(\epsilon\cdot\sqrt{m\log(d/\delta)})$.
Using Lemma 34, it follows that for $q$ i.i.d. copies
$\mathbf{S}_{1},...,\mathbf{S}_{q}$, the averaged estimator
$\frac{1}{q}\sum_{i=1}^{q}(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}\mathbf{A})^{-1}$
is an $(\epsilon^{\prime\prime},2\delta)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ for

$\epsilon^{\prime\prime}=\epsilon+O\bigg{(}\epsilon\cdot\sqrt{m\log^{2}(2dq/\delta)/q}\,\bigg{)}.$

Setting $q=O(m\log^{2}(d/\delta))$ and adjusting the constants appropriately,
we obtain the claim.

Note that in both Corollaries there is a slight interdependence in the
conditions for $m$ and $q$. This is in general unavoidable, since as $q$ grows
large with fixed $m$, the average has to eventually converge to the true
expectation of
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$,
which may be unbounded.

##  Appendix F Inversion bias lower bound for leverage score sampling

In this section, we show a lower bound on the inversion bias of approximate
leverage score sampling, proving Theorem 10. In the proof, we show a lower
bound for the inverse moment of a shifted Binomial random variable (Lemma 36),
which should be of independent interest.

###  F.1 Proof of Theorem 10

Without loss of generality, suppose that $n=2d$ (otherwise the matrix
$\mathbf{A}$ can be padded by zeros). We can also assume that $m\geq d$, since
the other cases follow easily. Our construction is designed so that uniform
row sampling is a $1/2$-approximation of leverage score sampling. Let
$\mathbf{S}$ be a uniform row sampling sketch of size $m$, i.e., its $i$th row
is $\sqrt{\\!\frac{n}{m}}\,\mathbf{e}_{s_{i}}^{\scriptscriptstyle{\top}}$,
where $s_{1},...,s_{m}$ are independent uniformly random indices from
$1,...,n$. Our matrix $\mathbf{A}$ consists of $n=2d$ scaled standard basis
vectors such that pairs of consecutive rows are given by
$\mathbf{a}_{2(i-1)+1}^{\scriptscriptstyle{\top}}=\mathbf{a}_{2(i-1)+2}^{\scriptscriptstyle{\top}}=\frac{1}{\sqrt{2}}\
\mathbf{e}_{i}^{\scriptscriptstyle{\top}}$ for $i\geq 2$, whereas the first
two rows are
$\mathbf{a}_{1}^{\scriptscriptstyle{\top}}=\frac{1}{\sqrt{4}}\mathbf{e}_{1}^{\scriptscriptstyle{\top}}$
and
$\mathbf{a}_{2}^{\scriptscriptstyle{\top}}=\sqrt{\frac{3}{4}}\mathbf{e}_{1}^{\scriptscriptstyle{\top}}$:

$\displaystyle\mathbf{A}=\begin{bmatrix}\frac{1}{\sqrt{4}}&&&\\\
\sqrt{\frac{3}{4}}&&0&\\\ &\frac{1}{\sqrt{2}}&&\\\ &\frac{1}{\sqrt{2}}&&\\\
&&\ddots&\\\ &0&&\frac{1}{\sqrt{2}}\\\ &&&\frac{1}{\sqrt{2}}\end{bmatrix}.$

First, note that $\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}=\mathbf{I}$,
and all of the squared row norms are within
$[\frac{1}{2}\frac{d}{n},\frac{3}{2}\frac{d}{n}]$, so uniform sampling is
indeed a $1/2$-approximate leverage score sampling scheme. Further, for any
$\gamma>0$, the matrix
$\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}$
is diagonal, and its diagonal entries are given by:

$\displaystyle\big{[}\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}\big{]}_{ii}=\begin{cases}\frac{\gamma
n}{m}\sum_{j=1}^{m}\big{(}\frac{1}{4}\mathbf{1}_{[s_{j}=1]}+\frac{3}{4}\mathbf{1}_{[s_{j}=2]}\big{)}=\frac{\gamma
n}{m}\cdot\frac{x+b_{1}/2}{2}&\text{for }i=1,\\\ \frac{\gamma
n}{m}\sum_{j=1}^{m}\big{(}\frac{1}{2}\mathbf{1}_{[s_{j}=2(i-1)+1]}+\frac{1}{2}\mathbf{1}_{[s_{j}=2(i-1)+2]}\big{)}=\frac{\gamma
n}{m}\cdot\frac{b_{i}}{2}&\text{otherwise,}\end{cases}$

where $b_{i}$’s are all identically (but not independently) distributed as
${\mathrm{Binomial}}(m,1/d)$ and $x$ is distributed, conditionally on $b_{1}$,
as ${\mathrm{Binomial}}(b_{1},1/2)$. Here $b_{i}$ denote the number of times
$s_{j}\in\\{2i-1,2i\\}$, while $x$ denotes the number of times $s_{j}=2$. Due
to the symmetry of the problem, conditionally on a given value of $b_{1}$
(i.e., a given value of counts $s_{j}$ that are equal to either unity or two),
each $s_{j}\in\\{1,2\\}$ is distributed uniformly over $\\{1,2\\}$, hence the
value $x$ of counts $s_{j}$ that are equal to two is distributed as
${\mathrm{Binomial}}(b_{1},1/2)$. This leads to the claimed distributional
representation.

The key idea in the construction is that the first diagonal entry of the
sketch has more variance than the others, and thus it will also have more
inversion bias. As a result, there is no scaling $\gamma$ that will
simultaneously correct the inversion bias of the first entry and of all the
other entries. To that end, we lower bound a shifted inverse moment of the
Binomial distribution in the following lemma, potentially of independent
interest, proven at the end of this section.

######  Lemma 36.

There is a universal constant $C>0$ such that for any positive integer $b$, if
$x\sim{\mathrm{Binomial}}(b,1/2)$ then:

$\displaystyle\mathbb{E}\bigg{[}\frac{1}{x+b/2}\bigg{]}\geq\Big{(}1+\frac{1}{Cb}\Big{)}\cdot\frac{1}{b}.$

Note that the expected inverse of
$\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}$
is undefined since the matrix may not be invertible. Thus, as in the
definition of an $(\epsilon,\delta)$-unbiased estimator, we must condition on
a high probability event which ensures invertibility. We start by considering
the largest such event, $\mathcal{E}^{*}:[\forall_{i}b_{i}>0]$. Using the fact
that, conditioned on $b_{1}$, the variable $x$ is independent of
$\mathcal{E}^{*}$, we have:

$\displaystyle\mathbb{E}\Big{[}\big{[}(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\big{]}_{11}\mid\mathcal{E}^{*}\Big{]}$
$\displaystyle=\left(\frac{\gamma
n}{m}\right)^{-1}\sum_{b>0}\mathbb{E}\Big{[}\frac{2}{x+b_{1}/2}\mid
b_{1}=b\Big{]}\,\mathrm{Pr}(b_{1}=b\mid\mathcal{E}^{*})$
$\displaystyle\overset{(a)}{\geq}\left(\frac{\gamma
n}{m}\right)^{-1}\sum_{b>0}\Big{(}1+\frac{1}{Cb}\Big{)}\,\frac{2}{b}\,\mathrm{Pr}(b_{1}=b\mid\mathcal{E}^{*})$
$\displaystyle=\sum_{b>0}\Big{(}1+\frac{1}{Cb}\Big{)}\,\mathbb{E}\Big{[}\big{[}(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\big{]}_{22}\mid
b_{2}=b\Big{]}\,\mathrm{Pr}(b_{2}=b\mid\mathcal{E}^{*})$
$\displaystyle\geq\mathbb{E}\Big{[}\big{[}(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\big{]}_{22}\mid\mathcal{E}^{*}\Big{]}$
$\displaystyle\quad+\frac{1}{2C}\frac{d}{m}\sum_{b=1}^{2m/d}\mathbb{E}\Big{[}\big{[}(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\big{]}_{22}\mid
b_{2}=b\Big{]}\,\mathrm{Pr}(b_{2}=b\mid\mathcal{E}^{*})$
$\displaystyle\overset{(b)}{\geq}\Big{(}1+\frac{d}{4Cm}\Big{)}\cdot\mathbb{E}\Big{[}\big{[}(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\big{]}_{22}\mid\mathcal{E}^{*}\Big{]},$

where in $(a)$ we used Lemma 36 and in $(b)$ we observed that
$\big{[}(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\big{]}_{22}$
decreases with $b_{2}$ and moreover, since $\mathbb{E}[b_{2}]=m/d\geq 1$, it
is easy to verify that the range $[1,2m/d]$ contains more than half of the
probability mass of ${\mathrm{Binomial}}(m,1/d)$.

The above derivation shows that when conditioned on $\mathcal{E}^{*}$, for any
scaling $\gamma>0$ the inversion bias will be at least $\Omega(d/m)$, since
the estimated matrix
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}=\mathbf{I}$ has the
same entries on the diagonal, whereas the expectation of the first two
diagonal entries of the estimator
$(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
differs by a factor of $1+\Omega(d/m)$. To complete the proof of Theorem 10,
it remains to show that the same is true not just for $\mathcal{E}^{*}$, but
for any event $\mathcal{E}\subseteq\mathcal{E}^{*}$ with sufficiently high
probability. Suppose that $\mathcal{E}$ is such an event, with
$\delta=\mathrm{Pr}(\mathcal{E}\mid\mathcal{E}^{*})\leq\mathrm{Pr}(\neg\mathcal{E})\leq\frac{1}{4C\cdot
16}(\frac{d}{m})^{2}$. Then, using $\tau_{i}=\frac{m}{\gamma
n}[(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}]_{ii}$
as a shorthand, we have:

$\displaystyle\mathbb{E}[\tau_{1}\mid\mathcal{E}]$
$\displaystyle=\mathbb{E}[\tau_{1}\mid\mathcal{E}^{*}]+\frac{\delta}{1-\delta}\Big{(}\mathbb{E}[\tau_{1}\mid\mathcal{E}^{*}]-\mathbb{E}[\tau_{1}\mid\mathcal{E}^{*},\neg\mathcal{E}]\Big{)}\geq\mathbb{E}[\tau_{1}\mid\mathcal{E}^{*}]-8\delta,$

where we used that $\delta\leq 1/2$ and, conditioned on $\mathcal{E}^{*}$, we
have $\tau_{1}\leq 4$. On the other hand,

$\displaystyle\mathbb{E}[\tau_{2}\mid\mathcal{E}]$
$\displaystyle=\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*}]+\frac{\delta}{1-\delta}\Big{(}\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*}]-\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*},\neg\mathcal{E}]\Big{)}\leq(1+2\delta)\,\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*}].$

Combining the two inequalities and using that
$\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*}]\geq d/m$ and
$\delta\leq\frac{1}{4C\cdot 16}(\frac{d}{m})^{2}$, we get:

$\displaystyle\frac{\mathbb{E}\big{[}[(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}]_{11}\mid\mathcal{E}\big{]}}{\mathbb{E}\big{[}[(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}]_{22}\mid\mathcal{E}\big{]}}$
$\displaystyle=\frac{\mathbb{E}[\tau_{1}\mid\mathcal{E}]}{\mathbb{E}[\tau_{2}\mid\mathcal{E}]}\geq\frac{(1+\frac{d}{4Cm})\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*}]-8\delta}{(1+2\delta)\mathbb{E}[\tau_{2}\mid\mathcal{E}^{*}]}$
$\displaystyle\geq\frac{1+\frac{d}{4Cm}-8\delta\frac{m}{d}}{1+2\delta}\geq\frac{1+\frac{d}{8Cm}}{1+\frac{d}{32Cm}}\geq
1+\frac{d}{64Cm}.$

Thus, as discussed above, we conclude that for any scaling $\gamma>0$ and any
event $\mathcal{E}$ with probability $\mathrm{Pr}(\mathcal{E})\geq
1-\frac{1}{4C\cdot 16}(\frac{d}{m})^{2}$, we have
$\|\mathbb{E}[(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\mid\mathcal{E}]-\mathbf{I}\|=\Omega(\frac{d}{m})$,
which concludes the proof.

###  F.2 Proof of Lemma 36

We conclude this section with a proof of the Binomial inverse moment bound
from Lemma 36. While existing work has focused on asymptotic expansions of
inverse moments of the Binomial [Zni09], those precise characterizations
either break down or appear to be impractical to work with when the variable
is significantly shifted, as in our case. Thus, we use a different strategy:
reducing the inverse moment bound to showing an anti-concentration inequality
for the Binomial distribution. For this, we use the classical Paley-Zygmund
inequality, stated below.

######  Lemma 37.

For any non-negative variable $\mathbf{Z}$ with finite variance and
$\theta\in(0,1)$, we have:

$\displaystyle\mathrm{Pr}\big{(}Z\geq\theta\,\mathbb{E}[Z]\big{)}\geq(1-\theta)^{2}\frac{\mathbb{E}[Z]^{2}}{\mathbb{E}[Z^{2}]}.$

Let $x\sim{\mathrm{Binomial}}(b,1/2)$ for a positive integer $b$. It follows
that:

$\displaystyle\mathbb{E}\Big{[}\frac{1}{x+b/2}-\frac{1}{b}\Big{]}$
$\displaystyle=\sum_{i=0}^{b}\mathrm{Pr}(x=i)\Big{(}\frac{1}{i+b/2}-\frac{1}{b}\Big{)}=\frac{1}{b}\sum_{i=0}^{b}\mathrm{Pr}(x=i)\frac{b/2-i}{b/2+i}$
$\displaystyle=\frac{1}{b}\sum_{i=0}^{\lfloor
b/2\rfloor}\mathrm{Pr}(x=i)(b/2-i)\Big{(}\frac{1}{b/2+i}-\frac{1}{3b/2-i}\Big{)},$

where the last equality is obtained by symmetrically pairing up the terms $i$
and $b-i$ in the first sum. Next, observe that for $0\leq i\leq
b/2-\sqrt{b}/4$, we have:

$\displaystyle(b/2-i)\Big{(}\frac{1}{b/2+i}-\frac{1}{3b/2-i}\Big{)}$
$\displaystyle\geq\frac{\sqrt{b}}{4}\Big{(}\frac{1}{b-\sqrt{b}/4}-\frac{1}{b+\sqrt{b}/4}\Big{)}=\frac{\sqrt{b}}{4}\cdot\frac{\sqrt{b}/2}{b^{2}-b/16}\geq\frac{1}{8b}.$

Putting this together, we conclude that:

$\displaystyle\mathbb{E}\Big{[}\frac{1}{x+b/2}\Big{]}\geq\Big{(}1+\frac{1}{8b}\mathrm{Pr}\big{\\{}x-b/2\leq-\sqrt{b}/4\big{\\}}\Big{)}\cdot\frac{1}{b}.$
(8)

Thus, it suffices to show that, with constant probability, $x$ is smaller than
its mean, $b/2$, by at least $\sqrt{b}/4$. This follows from the Paley-Zygmund
inequality (Lemma 37) by setting $\mathbf{Z}=(x-b/2)^{2}$. Using standard
formulas for the second and fourth centered moment of the Binomial
distribution, we have $\mathbb{E}[Z]=b/4$ and
$\mathbb{E}[Z^{2}]=\frac{b}{4}(1+\frac{3b-6}{4})\leq 3b^{2}/16$. Therefore,
setting $\theta=1/4$ in Lemma 37, we obtain:

$\displaystyle\mathrm{Pr}\big{(}x-b/2\,\leq\,-\sqrt{b}/4\big{)}$
$\displaystyle=\frac{1}{2}\,\mathrm{Pr}\big{(}|x-b/2|\,\geq\,\sqrt{b}/4\big{)}=\frac{1}{2}\,\mathrm{Pr}\big{(}Z\geq\theta\,\mathbb{E}[Z]\big{)}$
$\displaystyle\geq\frac{1}{2}\,\Big{(}1-\frac{1}{4}\Big{)}^{2}\,\frac{b^{2}/16}{3b^{2}/16}=\frac{3}{32}.$

Combining this with (8), we obtain the desired claim for $C=8\cdot 32/3$.

##  Appendix G Exact bias-correction for orthogonally invariant embeddings

In this section we prove that orthogonal invariance implies no inversion bias.
This claim has been mentioned in the main text, in Section 2. Here we give a
formal statement.

######  Proposition 38 (Orthogonal invariance implies no inversion bias).

Let $\mathbf{S}$ be a random and right-orthogonally invariant matrix;
specifically an $m\times n$ matrix (with $m\leq n$) such that for any
orthogonal $n\times n$ matrix $\mathbf{O}$, we have
$\mathbf{S}\stackrel{{\scriptstyle d}}{{=}}\mathbf{S}\mathbf{O}$. Assume that
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
exists with probability one. Then the inversion bias is exactly correctable,
i.e., there exists a constant $c=c_{m,n,d}$ such that
$\mathbb{E}\mathbf{\hat{\Sigma}}^{-1}=c\cdot\mathbf{\Sigma}^{-1}$; where
$\mathbf{\Sigma}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$ and
$\mathbf{\hat{\Sigma}}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}$.

Examples of orthogonal ensembles can be constructed in the following way:

  1. 1.

Let $\mathbf{S}$ have i.i.d. normal entries with variance $m^{-1}$. Due to the
properties of the Wishart ensemble, the constant $c_{m,n,d}$ is
$c_{m,n,d}=m/(m-d-1)$.

  2. 2.

Let $\mathbf{S}_{u}$ be a uniformly random $m\times n$ partial orthogonal
matrix (with $m\leq n$) such that
$\mathbf{S}_{u}\mathbf{S}_{u}^{\scriptscriptstyle{\top}}=\mathbf{I}_{m}$.
Equivalently, these are the first few rows of a Haar matrix. Then define
$\mathbf{S}=\sqrt{n/m}\cdot\mathbf{S}_{u}$, scaled such that
$\mathbb{E}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}=\mathbf{I}_{m}$. We
will call this the Haar sketch.

  3. 3.

The class of orthogonally invariant matrices has several closure properties.
Specifically, it is closed with respect to left-multiplication by any
matrices, right-multiplication by orthogonal matrices, and with respect to
vector space operations (addition and multiplication by scalars). Several
examples can be obtained this way. For instance, matrices $\mathbf{S}$ of the
form $\mathbf{S}=\mathbf{M}\mathbf{Z}$, where $\mathbf{Z}$ has i.i.d. normal
entries with variance $m^{-1}$, and $\mathbf{M}$ is an arbitrary matrix fixed
or random and independent of $\mathbf{Z}$ are orthogonally invariant.

###### Proof.

We start with a reduction to orthogonal matrices: Let
$\mathbf{A}=\mathbf{U}\mathbf{\Lambda}\mathbf{V}^{\scriptscriptstyle{\top}}$
be the SVD of $\mathbf{A}$. Here recall that $\mathbf{A}$ is an $n\times d$
matrix, with $n\geq d$ and with full column rank, and thus $\mathbf{U}$ is an
$n\times d$ partial orthogonal matrix with $n\geq d$, $\mathbf{\Lambda}$ is
$d\times d$ diagonal, and $\mathbf{V}$ is $d\times d$ orthogonal. Our goal is
to show that
$\mathbb{E}\mathbf{\hat{\Sigma}}^{-1}=c\cdot\mathbf{\Sigma}^{-1}$, or
equivalently that

$\mathbb{E}(\mathbf{V}\mathbf{\Lambda}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}\mathbf{\Lambda}\mathbf{V}^{\scriptscriptstyle{\top}})^{-1}=c\cdot(\mathbf{V}\mathbf{\Lambda}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{\Lambda}\mathbf{V}^{\scriptscriptstyle{\top}})^{-1}.$

Then, by cancelling $\mathbf{\Lambda}$ and $\mathbf{V}$ above (using that they
are deterministic square invertible matrices), and using that
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{I}$, we see that the
above inequality is equivalent to

$\mathbb{E}(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}=c\cdot\mathbf{I}.$

Thus, the problem is reduced to studying orthogonal matrices
$\mathbf{A}=\mathbf{U}$, such that
$\mathbf{\Sigma}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}=\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{I}$.
We claim that the right-orthogonal invariance implies that
$\mathbf{S}\mathbf{U}\stackrel{{\scriptstyle
d}}{{=}}\mathbf{S}\mathbf{U}\mathbf{O}$ for any $d\times d$ orthogonal matrix
$\mathbf{O}$. Here is a geometric argument. We have that
$\mathbf{S}\mathbf{U}$ are the angles that the random orthogonal rows of
$\mathbf{S}$ form with the fixed set of basis vectors formed by the columns of
$\mathbf{U}$. Also, $\mathbf{S}\mathbf{U}\mathbf{O}$ corresponds to the same
quantity, but with respect to the basis formed by $\mathbf{U}\mathbf{O}$.
Since $\mathbf{S}$ is right-rotationally invariant, these angles have the same
distribution.

Another, more algebraic proof is as follows. Since $\mathbf{S}$ is right-
rotationally invariant, for any orthogonal $n\times n$ matrix $\mathbf{R}$, we
have $\mathbf{S}\stackrel{{\scriptstyle d}}{{=}}\mathbf{S}\mathbf{R}$. Thus,
for any fixed matrix $\mathbf{U}$, we have
$\mathbf{S}\mathbf{U}\stackrel{{\scriptstyle
d}}{{=}}\mathbf{S}\mathbf{R}\mathbf{U}$. Choose a rotation matrix $\mathbf{R}$
such that
$\mathbf{R}\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}=\mathbf{U}\mathbf{O}\mathbf{U}^{\scriptscriptstyle{\top}}$,
while $\mathbf{R}\mathbf{U}^{\perp}$ is arbitrary, where $\mathbf{U}^{\perp}$
is an orthogonal complement of $\mathbf{U}$. Then, multiplying the above with
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}$ we have

$\mathbf{S}\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}\stackrel{{\scriptstyle
d}}{{=}}\mathbf{S}\mathbf{R}\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{S}\mathbf{U}\mathbf{O}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{S}\mathbf{U}\mathbf{O}.$

We get that $\mathbf{S}\mathbf{U}\stackrel{{\scriptstyle
d}}{{=}}\mathbf{S}\mathbf{U}\mathbf{O}$. Next,
$\mathbf{S}\mathbf{U}\stackrel{{\scriptstyle
d}}{{=}}\mathbf{S}\mathbf{U}\mathbf{O}$ implies that:

$\displaystyle\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\mathbf{O}^{\scriptscriptstyle{\top}}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}\mathbf{O}$
$\displaystyle(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\mathbf{O}^{\scriptscriptstyle{\top}}(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}\mathbf{O}$
$\displaystyle\mathbb{E}(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}$
$\displaystyle=\mathbf{O}^{\scriptscriptstyle{\top}}\mathbb{E}(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}\mathbf{O}.$

Thus, since
$\mathbf{J}:=\mathbb{E}(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}$
is preserved under conjugation by any orthogonal matrix, $\mathbf{J}$ must be
a multiple of the identity matrix, so $\mathbf{J}=c\mathbf{I}_{d}$, for some
$c=c_{m,n,d}$. This finishes the proof. ∎