# Sparse sketches with small inversion bias

Michał Dereziński Department of Statistics, University of California, Berkeley
<EMAIL_ADDRESS>Zhenyu Liao ICSI and Department of Statistics,
University of California, Berkeley<EMAIL_ADDRESS>Edgar Dobriban
Department of Statistics and Data Science, University of Pennsylvania
<EMAIL_ADDRESS>Michael W. Mahoney ICSI and Department of
Statistics, University of California, Berkeley<EMAIL_ADDRESS>
###### Abstract

For a tall $n\times d$ matrix $\mathbf{A}$ and a random $m\times n$ sketching
matrix $\mathbf{S}$, the sketched estimate of the inverse covariance matrix
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ is typically biased:
$\mathbb{E}[(\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}]\neq(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$,
where $\tilde{\mathbf{A}}=\mathbf{S}\mathbf{A}$. This phenomenon, which we
call _inversion bias_ , arises, e.g., in statistics and distributed
optimization, when averaging multiple independently constructed estimates of
quantities that depend on the inverse covariance matrix. We develop a
framework for analyzing inversion bias, based on our proposed concept of an
$(\epsilon,\delta)$-unbiased estimator for random matrices. We show that when
the sketching matrix $\mathbf{S}$ is dense and has i.i.d. sub-gaussian
entries, then after simple rescaling, the estimator
$(\frac{m}{m-d}\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}$
is $(\epsilon,\delta)$-unbiased for
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ with a sketch of size
$m=O(d+\sqrt{d}/\epsilon)$. In particular, this implies that for $m=O(d)$, the
inversion bias of this estimator is $O(1/\sqrt{d})$, which is much smaller
than the $\Theta(1)$ approximation error obtained as a consequence of the
subspace embedding guarantee for sub-gaussian sketches. We then propose a new
sketching technique, called LEverage Score Sparsified (LESS) embeddings, which
uses ideas from both data-oblivious sparse embeddings as well as data-aware
leverage-based row sampling methods, to get $\epsilon$ inversion bias for
sketch size $m=O(d\log d+\sqrt{d}/\epsilon)$ in time
$O({\mathrm{nnz}}(\mathbf{A})\log n+md^{2})$, where ${\mathrm{nnz}}$ is the
number of non-zeros. The key techniques enabling our analysis include an
extension of a classical inequality of Bai and Silverstein for random
quadratic forms, which we call the _Restricted Bai-Silverstein inequality_ ;
and anti-concentration of the Binomial distribution via the Paley-Zygmund
inequality, which we use to prove a lower bound showing that leverage score
sampling sketches generally do not achieve small inversion bias.

##  1 Introduction

Sketching has been widely used in the design of scalable algorithms, perhaps
most prominently in Randomized Numerical Linear Algebra (RandNLA) due to
applications in machine learning and data analysis. In this approach, one
randomly samples or computes a random projection of the data matrix to
construct a smaller matrix, the sketch. One then uses the sketch as a
surrogate to approximate quantities of interest. The analysis of these methods
typically proceeds via a Johnson-Lindenstrauss-type argument to establish that
the geometry of the matrix is not perturbed too much under the sketching
operation. These methods have yielded state-of-the-art in worst-case analysis,
high-quality numerical implementations, and numerous applications in machine
learning [Mah11, HMT11, Woo14, DM16, DM18].

In many cases, either to preserve the structure of the data or for algorithmic
reasons, one is interested in sparse sketches, i.e., random transformations
that are represented by matrices with most entries exactly equal to zero. One
class of sparse sketches includes row sampling techniques, such as leverage
score sampling [DMM08, DMIMW12, MMY15], which are typically _data-aware_ , in
the sense that the sampling distribution depends on the given data matrix.
Another important class of methods uses _data-oblivious_ sparse embeddings,
such as the CountSketch [CCFC02, CW17, MM13, NN13], to construct sketches in
time depending on the number of non-zeros (nnz) in the input.

In all these cases, one can show that the sketch will be an approximation of
the solution with high probability. However, comparatively little is known
about the average performance of these sketches. In particular, there may be a
systematic bias away from the solution, which is problematic in many
situations in statistics, machine learning, and data analysis. Perhaps the
most ubiquitous example of this phenomenon is the systematic bias caused by
matrix inversion, a key component of algorithms in the aforementioned domains.
In this paper, we introduce the fundamental notion of inversion bias, which
provides a finer control over the sketched estimates involving matrix
inversion. We show that one can conveniently make the inversion bias small
with dense Gaussian and sub-gaussian sketches. We also show that some sparse
sketches do not have this desired property. Then, we provide a non-trivial new
construction and algorithm, using ideas from both data-oblivious projections
and data-aware sampling, to get small inversion bias even for very sparse
sketches.

###  1.1 Overview

Consider an $n\times d$ data matrix $\mathbf{A}$ of rank $d$, where $n\geq d$.
In many applications, we wish to approximate quantities of the form
$F((\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1})$, where
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ is the $d\times d$
inverse data covariance and $F(\cdot)$ is a linear functional. Our goal is to
provide a finer control over the effect of matrix inversion on the quality of
such estimates. Here are some of the motivating examples:

  * •

The vector $(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{b}$
is the solution of ordinary least squares (OLS) when
$\mathbf{b}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{y}$ for a vector
$\mathbf{y}$, arguably the most widely used multivariate statistical method
[And03, Rao73, HTF09], and it is also crucial for the Newton’s method in
numerical optimization [BV04, NW06]. In particular, accurate approximations of
this vector lead directly to improved convergence guarantees for many
optimization algorithms [PW16, WRXM18].

  * •

The scalar
$\mathbf{x}^{\scriptscriptstyle{\top}}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{x}$
for a vector $\mathbf{x}$, has numerous use-cases: When
$\mathbf{x}=\mathbf{a}_{i}$ is one of the rows of $\mathbf{A}$, then it
represents the statistical leverage scores [DMM06]; If
$\mathbf{x}=\mathbf{e}_{i}$ is a standard basis vector, then this is the
squared length of the confidence interval for the $i$-th coefficient in OLS
[And03, HTF09].

  * •

The scalar
$\mathrm{tr}\,\mathbf{C}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$
for a matrix $\mathbf{C}$, is used to quantify uncertainty in statistical
results, e.g., via the mean squared error (MSE) of estimating the regression
coefficients in OLS [And03, HTF09], and to formulate widely used criteria from
experimental design, e.g., A-designs and V-designs [Puk06, CR00].

More generally, our work is also motivated by the important problem of inverse
covariance estimation in statistics, machine learning, finance, signal
processing, and related areas [Dem72, MB06, FHT08, LF09, CLL11, LW12, MH15].
In this area, we wish to estimate statistically the inverse covariance matrix
of a population, or some of its functionals, based on a finite number of
samples. Furthermore, inverting covariance matrices occurs in Bayesian
statistics [Har69, GCS+13], Gaussian processes [Ras03], as well as time series
analysis and control, e.g., via the Kalman filter [WB95, BD09].

When $n$ and $d$ are large, and particularly when $n\gg d$, then the costs of
storing the matrix $\mathbf{A}$ and of computing
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ are prohibitively
large. Matrix sketching has proven successful at drastically reducing these
costs by approximating the inverse covariance with a sketched estimate
$(\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}$ based
on a smaller matrix $\tilde{\mathbf{A}}=\mathbf{S}\mathbf{A}$, where
$\mathbf{S}$ is a random $m\times n$ matrix and $m\ll n$ [Mah11, HMT11, Woo14,
DM16, DM18]. As a concrete algorithmic motivation for our work, consider the
following popular strategy for boosting the quality of such estimates:
Construct multiple copies in parallel, based on independent sketches, and then
average the estimates. This strategy is especially useful in distributed
architectures, where storage and computing resources are spread out across
many machines, and has commonly appeared in the literature [KBRR16, KBY+16,
WRXM18, DM19]. While promising in practice, this averaging technique is
fundamentally limited by the _inversion bias_ : even though the sketched
covariance estimate is unbiased,
$\mathbb{E}[\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}]=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$,
its inverse in general is not unbiased, i.e.,
$\mathbb{E}[(\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}]\neq(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$.
When the sketch size $m$ is not much larger than the dimension $d$, the size
of this bias can be very significant, even as large as the approximation
error, in which case averaging becomes ineffective. Motivated by this, we ask:

> When is the inversion bias small, relative to the approximation error?

In this paper, we develop a framework for analyzing the inversion bias of
sketching, via the notion of an $(\epsilon,\delta)$-unbiased estimator
(Definition 3), and we show how it can be used to provide improved
approximation guarantees for averaging. Through this framework, we provide
several contributions towards addressing the above question.

Sub-gaussian sketches have small inversion bias. Arguably the most classical
family of sketches consists of dense random matrices $\mathbf{S}$ with i.i.d.
sub-gaussian entries. These sketches offer strong relative error approximation
guarantees via the so-called _subspace embedding_ property, at the expense of
high computational cost of the matrix product
$\tilde{\mathbf{A}}=\mathbf{S}\mathbf{A}$. We show that, upon a simple
correction, sub-gaussian sketches are nearly-unbiased, i.e., their inversion
bias is much smaller than the approximation error, which means that averaging
can be used to significantly improve the approximation quality. In particular,
we show that, after a simple scalar rescaling, the inverse covariance
estimator of the form
$(\frac{m}{m-d}\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}$
achieves $\epsilon$ inversion bias relative to
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ with a sketch of size
only $m=O(d+\sqrt{d}/\epsilon)$ (Proposition 4). In contrast, to ensure that
$(\frac{m}{m-d}\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}$
is an ${\eta}$ relative error approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ via the subspace
embedding property, we need a sub-gaussian sketch of size
$m=\Theta(d/{\eta}^{2})$, which is comparatively larger if we let
$\eta=\Theta(\epsilon)$. This implies that an aggregate estimator obtained via
averaging can with high probability produce a relative error approximation
that is by a factor of $O(1/\sqrt{m})$ better than the approximation error
offered by any one of the estimators being averaged.

LEverage Score Sparsified (LESS) embeddings. We show that existing
algorithmically efficient sketching techniques may not provide guarantees for
the inversion bias that match those satisfied by dense sub-gaussian sketches
(see Theorem 10 for a lower bound on leverage score sampling, and a discussion
of other methods in Appendix C.3). To address this, we propose a new family of
sketching methods, called LEverage Score Sparsified (LESS) embeddings, which
combines a data-oblivious sparsification strategy reminiscent of the
CountSketch with the data-aware approach of approximate leverage score
sampling. LESS embeddings have time complexity
$O({\mathrm{nnz}}(\mathbf{A})\log n+md^{2})$ and achieve $\epsilon$ inversion
bias with the sketch of size $m=O(d\log d+\sqrt{d}/\epsilon)$, nearly matching
our guarantee for sub-gaussian sketches (Theorem 8). Thus, our new algorithm
provides a promising way to address the fundamental problem of inversion bias,
and it may have many other applications in the future. Finally, our analysis
reveals two structural conditions for small inversion bias (Theorem 11), one
of which (Condition 2, called the Restricted Bai-Silverstein condition) leads
to a generalization of a classical inequality used in random matrix theory,
and should be of independent interest.

###  1.2 Related work

#### Distributed averaging.

Averaging strategies have been studied extensively in the literature,
particularly in the context of machine learning and numerical optimization.
This line of work has proven particularly effective for _federated learning_
[KBRR16, KBY+16], where local storage and communication bandwidth are
particularly constrained. The performance of averaged estimates was analyzed
in numerous statistical learning settings [MMS+09, MHM10, ZDW13, DS18, DS20]
and in stochastic first-order optimization [ZWLS10, AD11]. Of particular
relevance to our results is a recent line of works on distributed second-order
optimization [SSZ14, ZL15, RKR+16, WRXM18], as well as large-scale second-
order optimization [YGKM19, YGS+20], since sketching is used there to estimate
(implicitly) the inverse Hessian matrix which arises in Newton-type methods.
In particular, [DM19, DBPM20] pointed to Hessian inversion bias as a key
challenge in these approaches. To address it, their algorithms use non-i.i.d.
sampling sketches based on Determinantal Point Processes (DPPs) [DM21]. DPP-
based sketches are known to correct inversion bias exactly [DW18, DWH19,
DCMW19]. However, state-of-the-art DPP sampling algorithms [DWH18, DCV19,
CDV20] have time complexity $O({\mathrm{nnz}}(\mathbf{A})\log n+d^{4}\log d)$,
which is considerably more expensive than fast sketching techniques when
dimension $d$ is large.

#### Random matrix theory.

When considering $\mathbf{S}\in\mathbb{R}^{m\times n}$ having i.i.d. zero-mean
rows,
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}$
can be viewed as the popular _sample covariance estimator_ of the “population
covariance matrix”
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}\in\mathbb{R}^{d\times d}$. In
this area, one often considers the matrix
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}-z\mathbf{I})^{-1}$
for $z\in\mathbb{C}\setminus\mathbb{R}_{+}$, the so-called _resolvent_ matrix,
which plays a fundamental role in the literature of random matrix theory (RMT)
[MP67, BS10, ER05, AGZ10, CD11, Tao12, BBP17] and which is directly connected
to the popular Marchenko-Pastur law [MP67]. The RMT literature focuses on the
Stieltjes transform (that is, the normalized trace of the resolvent) to
investigate the limiting eigenvalue distribution of large random matrices of
the form
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}$
as $m,n,d\to\infty$ at the same rate. Here, we provide _precise_ and _finite-
dimensional_ results on the inverse sketched matrix. This addresses the
important case of $z=0$, which is typically avoided in RMT analyses, due to
the difficulty of dealing with the possible singularity. More generally, the
resolvent also appears as the key object of study in the spectrum analysis of
linear operators in general Hilbert space [AG13], as well as in modern convex
optimization theory [BC11], thereby showing a much broader interest of the
proposed analysis.

#### Sketching.

For overviews of sketching and random projection methods, we refer to [Vem05,
HMT11, Mah11, Woo14, DM16, DM17, DM18, DM21]. A key result in this area is the
Johnson-Lindenstrauss lemma, which states that norms, and thus also relative
distances between points, are _approximately_ preserved after sketching, i.e.,
$(1-\eta)\|\mathbf{x}_{i}\|^{2}\leq\|\mathbf{S}\mathbf{x}_{i}\|^{2}\leq(1+\eta)\|\mathbf{x}_{i}\|^{2}$
for $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in\mathbb{R}^{p}$. This is further
extended to the _subspace embedding property_ : for all $\mathbf{x}$, the norm
of $\mathbf{x}$ is preserved up to an $\eta$ factor. Subspace embeddings were
first used in RandNLA by [DMM06], where they were used in a data-aware context
to obtain relative-error approximations for $\ell_{2}$ regression and low-rank
matrix approximation [DMM08]. Subsequently, data-oblivious subspace embeddings
were used by [Sar06] and popularized by [Woo14]. Both data-aware and data-
oblivious subspace embeddings can be used to derive bounds for the accuracy of
various algorithms [DM16, DM18].

The most popular sketching methods include random projections with i.i.d.
entries, random sampling of the datapoints, uniform orthogonal projections,
Subsampled Randomized Hadamard Transform (SRHT) [Sar06, AC06], leverage score
sampling [DMM08, DMIMW12, MMY15], and CountSketch [CCFC02, CW17, NN13, MM13].
Random projection based approaches have been developed for a wide variety of
problems in data science, statistics, machine learning etc., including linear
regression [Sar06, DMMS11, RM16, DL18], ridge regression [LDFU13, CLL+15,
WGM18, LD19], two sample testing [LJW11, SLR16], classification [CS17], PCA
[FKV04, DKM06, Sar06, LWM+07, HMST11, HMT11, WLRT08, MM15, TYUC17, DSBN15,
YLDW20, GWS20], convex optimization [PW15, PW16, PW17], etc.; see [Woo14,
DM16, DM18] for a more comprehensive list. Our new LESS embeddings have the
potential to be relevant for all those important applications.

##  2 Dense Gaussian and sub-gaussian sketches have small inversion bias

Consider first the classical Gaussian sketch, i.e., where the entries of
$\mathbf{S}$ are i.i.d. standard normal scaled by $1/\sqrt{m}$. In this
special case, the sketched covariance matrix
$\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$ is a
Wishart-distributed random matrix, and we have:

$\displaystyle\mathbb{E}\big{[}(\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}\big{]}=\tfrac{m}{m-d-1}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\quad\text{for}\quad
m\geq d+2.$ (1)

In other words, even though the sketched inverse covariance is not an unbiased
estimate, the bias can be corrected by simply scaling the matrix, after which
averaging can be used effectively without encountering any inversion bias.

The key property which enables exact bias-correction for the Gaussian sketch
is _orthogonal invariance_. This property requires that for any orthonormal
matrix $\mathbf{O}$, the distributions of the random matrices $\mathbf{S}$ and
$\mathbf{S}\mathbf{O}$ are identical. An example beyond Gaussians are Haar
sketches, which are uniform over all partial orthogonal matrices. If a sketch
$\mathbf{S}$ is orthogonally invariant and
$\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$ is
invertible with probability one, then we can show that the inversion bias can
be corrected exactly, in that, (1) holds with some constant factor $c$
(replacing the factor $\frac{m}{m-d-1}$) that depends on the distribution of
the sketch (see Proposition 38 in Appendix G).

Exact bias-correction, achieved by the Gaussian sketch and other orthogonally
invariant sketches, is no longer possible for general sub-gaussian sketches.
Here, we consider sketching matrices with i.i.d. entries that (after scaling
by $\sqrt{m}$) have $O(1)$ sub-gaussian Orlicz norm. Consider for example the
so-called Rademacher sketch, with $\mathbf{S}$ consisting of scaled i.i.d.
random sign entries (which is useful for reducing the cost of randomness
relative to the Gaussian sketch). In this case, an exact bias-correction
analogous to (1) is clearly infeasible for any $d>1$, simply because, with
some positive (but exponentially small) probability, the matrix
$\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$ will be non-
invertible, making the expectation undefined. Yet, any task where we observe
at most polynomially many independent estimates (such as averaging) should not
be affected by such low-probability events, so we need a notion of near-
unbiasedness that is robust to this. To that end, we first recall a standard
definition of a relative error approximation for a positive semi-definite
matrix.

######  Definition 1 (Relative error approximation).

A positive semi-definite (p.s.d.) matrix $\tilde{\mathbf{C}}$ (or a non-
negative scalar) is an ${\eta}$-approximation of $\mathbf{C}$, denoted as
$\tilde{\mathbf{C}}\approx_{\eta}\mathbf{C}$, if

$\mathbf{C}/(1+{\eta})\preceq\tilde{\mathbf{C}}\preceq(1+{\eta})\cdot\mathbf{C}.$

If $\tilde{\mathbf{C}}$ is random and the above holds with probability
$1-\delta$, then we call it an $({\eta},\delta)$-approximation.

######  Remark 2 (Subspace embedding).

If
$\tilde{\mathbf{C}}=\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$
where $\tilde{\mathbf{A}}\in\mathbb{R}^{m\times d}$ is a sketch of
$\mathbf{A}\in\mathbb{R}^{n\times d}$, then the condition
$\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}\approx_{\eta}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$
is called the subspace embedding property with error $\eta$.

For instance, any sketching matrix $\mathbf{S}$ with i.i.d. $O(1)$ sub-
gaussian random entries, of size $m=O((d+\ln(1/\delta))/{\eta}^{2})$, where
${\eta}\in(0,1)$, ensures that $\tilde{\mathbf{A}}=\mathbf{S}\mathbf{A}$ with
probability $1-\delta$ satisfies the subspace embedding property with error
${\eta}$. In other words,
$\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$ is an
$({\eta},\delta)$-approximation of
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$ (This is known to be tight;
see, e.g., [NN14]). As a consequence, the same guarantee applies to the
inverse
$(\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}$,
relative to $(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$. The
$\delta$ failure probability makes this definition robust to the rare events
where $\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$ is not
invertible. It is natural to desire a similar robustness in the definition of
near-unbiasedness. We achieve this as follows.

######  Definition 3 ($(\epsilon,\delta)$-unbiased estimator).

A random p.s.d. matrix $\tilde{\mathbf{C}}$ is an $(\epsilon,\delta)$-unbiased
estimator of $\mathbf{C}$ if there is an event $\mathcal{E}$ that holds with
probability $1-\delta$ such that

$\displaystyle\mathbb{E}\big{[}\tilde{\mathbf{C}}\mid\mathcal{E}\big{]}\approx_{\epsilon}\mathbf{C},\quad\text{and}\quad\tilde{\mathbf{C}}\preceq
O(1)\cdot\mathbf{C}\quad\text{when conditioned on $\mathcal{E}$.}$

Note that this definition only becomes meaningful if we use it with an
$\epsilon$ that is much smaller than the approximation error ${\eta}$ in
Definition 1 (for instance, we will often have ${\eta}=\Omega(1)$ and
$\epsilon\ll 1$). Further, note the following two important aspects of
Definition 3. First, instead of a simple expectation, we condition on some
high probability event $\mathcal{E}$, which, similarly as in Definition 1,
allows robustness against such corner cases as when the sketch
$\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}$ is not
invertible. Second, conditioned on the event $\mathcal{E}$, in addition to an
$\epsilon$-approximation holding in expectation, we require a weaker upper
bound to hold almost surely, in terms of the target matrix $\mathbf{C}$ scaled
by some constant factor. This condition is important to guard against certain
corner cases where the probability mass is extremely skewed. For instance,
suppose that $\tilde{C}$ is a scalar random variable which is uniform over
$[0,1]$ and has an additional probability mass of $10^{-10}$ at the value
$10^{100}$. Here, averaging will not prove effective at converging to the true
expectation of $\tilde{C}$, but we could still use the notion of
$(\epsilon,\delta)$-unbiasedness to show that the average of an appropriately
chosen number (much smaller than $10^{10}$) of i.i.d. copies will converge
very close to $0.5$, by choosing an event $\mathcal{E}$ that avoids the
$10^{100}$ (see Appendix E).

We are now ready to state our main result for sub-gaussian sketches (this is
in fact a corollary of our more general result, Theorem 11, discussed in
Section 4), which asserts that after proper rescaling, not only the Gaussian
sketch, but in fact all sub-gaussian sketches (including the Rademacher
sketch) enjoy small inversion bias.

######  Proposition 4 (Near-unbiasedness of sub-gaussian sketches).

Let $\mathbf{S}$ be an $m\times n$ random matrix such that
$\sqrt{m}\,\mathbf{S}$ has i.i.d. $O(1)$-sub-gaussian entries with mean zero
and unit variance. There is $C=O(1)$ such that for any
$\epsilon,\delta\in(0,1)$ if $m\geq C(d+\sqrt{d}/\epsilon+\log(1/\delta))$,
then for all $\mathbf{A}\in\mathbb{R}^{n\times d}$ of rank $d$,
$(\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}$.

Observe that the scaling $\frac{m}{m-d}$ essentially matches the exact bias-
correction for Gaussian sketches, which is $\frac{m}{m-d-1}$. In fact, the
same statement of the theorem holds with either scaling, and we merely chose
the simplest form of the scaling.

As a corollary of the near-unbiasedness of sub-gaussian sketches, we can show
the following approximation guarantee for averaging the inverse covariance
matrix estimates. Recall that our primary motivation is parallel and
distributed averaging, where the computational cost does not grow with the
number of independent estimates.

######  Corollary 5.

Let $\mathbf{S}$ be a sub-gaussian sketching matrix of size $m$, and let
$\mathbf{S}_{1},...,\mathbf{S}_{q}$ be i.i.d. copies of $\mathbf{S}$. There is
$C=O(1)$ such that if $m\geq C(d+\sqrt{d}/\epsilon+\log(q/\delta))$ and $q\geq
Cm\log(d/\delta)$, then for any $\mathbf{A}\in\mathbb{R}^{n\times d}$ of rank
$d$,
$\frac{1}{q}\sum_{i=1}^{q}(\tfrac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}\mathbf{A})^{-1}$
is an $(\epsilon,\delta)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$.

Proposition 4 shows that for a sub-gaussian sketch
$\tilde{\mathbf{A}}=\mathbf{S}\mathbf{A}$ of size $m\geq Cd$, the sketched
inverse covariance
$(\frac{m}{m-d}\tilde{\mathbf{A}}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}})^{-1}$
has inversion bias $O(\sqrt{d}/m)$. This means that the inversion bias of this
estimator is smaller than the approximation error, which is
$\Theta(\sqrt{d/m})$, by a factor of $O(1/\sqrt{m})$. Thus, using Corollary 5,
we can reduce the approximation error by averaging $q=O(m\log(d/\delta))$
copies of this estimator, obtaining that
$\frac{1}{q}\sum_{i=1}^{q}(\tfrac{m}{m-d}\tilde{\mathbf{A}}_{i}^{\scriptscriptstyle{\top}}\tilde{\mathbf{A}}_{i})^{-1}$
is with high probability an $O(\sqrt{d}/m)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$. In particular, when
$m=\Theta(d)$, then the approximation error of a single estimate (without
averaging) is $\Theta(1)$, whereas the approximation error of the averaged
estimate is only $O(1/\sqrt{d})$.

##  3 Main results: Less inversion bias with LESS embeddings

To address the high computational cost of sub-gaussian sketches, while
preserving their good near-unbiasedness properties, we propose a new family of
sketches, which we call LEverage Score Sparsified (LESS) embeddings. A LESS
embedding is defined simply as a sparsified sub-gaussian sketch, where the
sparsification is designed so as to ensure small inversion bias for a
particular matrix $\mathbf{A}$. Our approach combines ideas from approximate
leverage score sampling (which is data-aware) with ideas from sparse embedding
matrices (which are normally data-oblivious). Importantly, neither strategy by
itself is sufficient to ensure small inversion bias (see our lower bound in
Theorem 10 and discussion in Section 4). Each row of a LESS embedding is
sparsified independently using a sparsification pattern defined as follows.
Recall that for a tall full rank matrix $\mathbf{A}$, we use
$\mathbf{a}_{i}^{\scriptscriptstyle{\top}}$ to denote the $i$th row of
$\mathbf{A}$, and the $i$th leverage score of $\mathbf{A}$ is defined as
$l_{i}=\mathbf{a}_{i}^{\scriptscriptstyle{\top}}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{a}_{i}$.

######  Definition 6 (LESS: LEverage Score Sparsified embedding).

Fix a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ of rank $d$ with leverage
scores $l_{1},...,l_{n}$, and let $s_{1},...,s_{d}$ be sampled i.i.d. from a
probability distribution $(p_{1},...,p_{n})$ such that
$p_{i}\approx_{O(1)}l_{i}/d$ for all $i$. Then, the random vector
${\boldsymbol{\xi}}^{\scriptscriptstyle{\top}}=\big{(}\sqrt{\\!\frac{b_{1}}{dp_{1}}},...,\sqrt{\\!\frac{b_{n}}{dp_{n}}}\big{)}$,
where $b_{i}=\sum_{t=1}^{d}1_{[s_{t}=i]}$, is called a leverage score
sparsifier for $\mathbf{A}$.

Sketching matrix $\mathbf{S}$ is a LESS embedding of size $m$ for a matrix
$\mathbf{A}$, if it consists of $m$ i.i.d. row vectors distributed as
$\frac{1}{\sqrt{m}}(\mathbf{x}\circ{\boldsymbol{\xi}})^{\scriptscriptstyle{\top}}$,
where $\circ$ denotes an entry-wise product and $\mathbf{x}$ is a random
vector with i.i.d. mean zero, unit variance, $O(1)$-sub-gaussian entries.

######  Remark 7 (Time complexity of LESS).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ of rank $d$, there is an
algorithm with an $O({\mathrm{nnz}}(\mathbf{A})\log n+d^{3}\log d)$ time
preprocessing step, that can then construct a LESS embedding sketch
$\mathbf{S}\mathbf{A}$ of size $m$ in time $O(md^{2})$. In the following
results we always use $m\geq d\log d$, in which case the total cost of
constructing a LESS embedding is $O({\mathrm{nnz}}(\mathbf{A})\log n+md^{2})$.

The matrix product $\mathbf{S}\mathbf{A}$ costs only $O(md^{2})$ because, by
definition, the number of non-zeros per row of $\mathbf{S}$ is bounded almost
surely by $d$. It is not essential for our analysis that we sample exactly $d$
indices in each row of a LESS embedding, but we fix it here for the sake of
simplicity. We could also have approximately $d$ non-zeros per row, and
similar results would still hold. To construct the distribution
$(p_{1},...,p_{n})$, the sparsifier requires a constant relative error
approximation of all the leverage scores of $\mathbf{A}$, which can be
computed in $O({\mathrm{nnz}}(\mathbf{A})\log n+d^{3}\log d)$ time [DMIMW12,
CW17]. Alternatively we can use our approach in a data-oblivious way, by
combining LEverage Score Sparsification with the Randomized Hadamard Transform
[AC09, DMMS11], which we may abbreviate as LESSRHT. Here, the matrix
$\mathbf{A}$ is first transformed so that it has approximately uniform
leverage scores [DMIMW12], and then we can sparsify it using a uniform
distribution, i.e., $p_{i}=1/n$ for all $i$, with total cost $O(nd\log
n+md^{2})$. Finally, computing the sketched inverse covariance matrix
estimator
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
only adds an $O(md^{2})$ cost. These costs can be further optimized using fast
matrix multiplication [Wil12].111The cost of computing the matrix product
$\mathbf{S}\mathbf{A}$ can be optimized beyond $O(md^{2})$ by adapting the
fast matrix multiplication routines to take advantage of the sparsity pattern;
see, e.g., [YZ05].

In our main result, we show that LESS embeddings enjoy small inversion bias,
nearly matching our guarantee for sub-gaussian sketches (Proposition 4).

######  Theorem 8 (Near-unbiasedness for LESS).

Suppose that $\mathbf{S}$ is a LESS embedding of size $m$ for a rank $d$
matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$. There is $C=O(1)$ such that if
$m\geq C(d\log(d/\delta)+\sqrt{d}/\epsilon)$ then the sketch
$(\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}$.

Thus, we show that the inversion bias guarantee for LESS embeddings matches
our result for sub-gaussian sketches up to a logarithmic factor. This
additional factor is standard in the analysis of fast sketching methods. It
comes from the fact that, as an artifact of the matrix concentration bounds
[Tro12] we use in our analysis of LESS embeddings, a sketch of size $m=O(d\log
d)$ is needed to satisfy the subspace embedding property, which is one of our
two structural conditions for small inversion bias (see Section 4). As a
corollary, we obtain an improved guarantee for parallel and distributed
averaging of i.i.d. sketched inverse covariance estimates which also matches
the corresponding statement for sub-gaussian sketches (Corollary 5) up to
logarithmic factors.

######  Corollary 9.

Let $\mathbf{S}$ be a LESS embedding matrix of size $m$ for a rank $d$ matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$, and let
$\mathbf{S}_{1},...,\mathbf{S}_{q}$ be i.i.d. copies of $\mathbf{S}$. There is
$C=O(1)$ such that if $m\geq C(d\log(q/\delta)+\sqrt{d}/\epsilon)$ and $q\geq
Cm\log^{2}(d/\delta)$, then
$\frac{1}{q}\sum_{i=1}^{q}(\tfrac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}\mathbf{A})^{-1}$
is an $(\epsilon,\delta)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$.

To motivate and place our new algorithm into context, we demonstrate that
existing fast sketching techniques may not achieve an inversion bias bound
comparable to that of sub-gaussian sketches, even if they achieve a nearly
matching subspace embedding guarantee. This lower bound demonstrates the
hardness of constructing an $(\epsilon,\delta)$-unbiased estimator of the
inverse covariance matrix from its sketch. We show this here for leverage
score sampling [DMM06, DMM08, DMIMW12, MMY15]. However, based on evidence from
our analysis, we conjecture that similar lower bounds hold for other methods
such as Subsampled Randomized Hadamard Transform [AC09, DMMS11] and data-
oblivious sparse embedding matrices [CW17, NN13, MM13].222An alternative
approach to achieving small inversion bias is to chain together a fast sketch
having a larger size, say, $t=\widetilde{O}(d/\epsilon^{2})$, with a sub-
gaussian sketch having a smaller size $m=O(d+\sqrt{d}/\epsilon)$. However,
this leads to a sub-optimal time complexity in terms of the polynomial
dependence on $d$ due to the cost $O(tmd)$ of the sub-gaussian sketch. For
example, with $\epsilon=1/\sqrt{d}$, the overall cost is
$\widetilde{O}({\mathrm{nnz}}(\mathbf{A})+d^{4})$ compared to
$\widetilde{O}({\mathrm{nnz}}(\mathbf{A})+d^{3})$ with LESS.

######  Theorem 10 (Lower bound for leverage score sampling).

For any $n\geq 2d\geq 4$, there is an $n\times d$ matrix $\mathbf{A}$ and a
row sampling $(p_{1},...,p_{n})$, with a corresponding $m\times n$ sketching
matrix $\mathbf{S}$, s.t.:

  1. 1.

The row sampling $(p_{1},...,p_{n})$ is a $1/2$-approximation of leverage
score sampling; and

  2. 2.

For any sketch size $m$ and scaling $\gamma$,
$(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
is not an $(\epsilon,\delta)$-unbiased estimator of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ with any
$\epsilon\leq c\frac{d}{m}$ and $\delta\leq c(\frac{d}{m})^{2}$, where $c>0$
is an absolute constant.

In the proof of Theorem 10, we develop a new lower bound for the inverse
moment of the Binomial distribution (Lemma 36), by using anti-concentration of
measure via the Paley-Zygmund inequality, which should be of independent
interest. To illustrate Theorem 10, consider a sketch of size $m=O(d\log d)$.
This is sufficient to ensure that approximate leverage score sampling achieves
the subspace embedding property with relative error $O(1)$. In particular, it
implies that for any $\gamma=\Theta(1)$, the inverse covariance matrix
estimator
$(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}$
is with high probability an $O(1)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$. Our lower bound
implies that the inversion bias of any such estimator is $\Omega(1/\log d)$,
which is up to logarithmic factors the same as the approximation achieved by a
single estimator.

Thus, Theorem 10 shows that when $m=O(d\log d)$, averaging i.i.d. copies of
the sketched inverse covariance estimator obtained from approximate leverage
score sampling may lead to only $\Omega(1/\log d)$ factor improvement in the
approximation, which is merely inverse-logarithmic in $d$. In contrast,
Theorem 8 shows that, when using our new LESS embeddings with the same sketch
size and time complexity, averaging i.i.d. copies of the sketched inverse
covariance reduces the approximation error by a factor of $O(1/\sqrt{d})$,
which is inverse-polynomial in $d$ and thus far superior to what is achievable
by approximate leverage score sampling.

##  4 Our techniques: Structural conditions for near-unbiasedness

In order for our analysis of inversion bias to apply to a wide range of
sketching techniques, we give two key structural conditions for a random
sketching matrix $\mathbf{S}$ that are sufficient to achieve provably small
inversion bias. The first is the subspace embedding property discussed in
Remark 2, which we now use as one of the key conditions needed in our
analysis.

######  Condition 1 (Subspace embedding).

The (sketching) matrix $\mathbf{S}\in\mathbb{R}^{m\times n}$ satisfies the
subspace embedding condition with $\eta\geq 0$ for a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$, if
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}\approx_{\eta}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$.

The second structural condition for small inversion bias is a property of each
individual row of $\mathbf{S}$. We use an $n$-dimensional random row vector
$\mathbf{x}^{\scriptscriptstyle{\top}}$ to denote the marginal distribution of
a row of $\mathbf{S}$ (after scaling by $\sqrt{m}$). This condition represents
a key novelty in our analysis.

######  Condition 2 (Restricted Bai-Silverstein).

The random vector $\mathbf{x}\in\mathbb{R}^{n}$ satisfies the Restricted Bai-
Silverstein condition with $\alpha>0$ for a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$, if
$\mathrm{Var}\\!\big{[}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}\big{]}\leq\alpha\cdot\mathrm{tr}(\mathbf{B}^{2})$
for all p.s.d. matrices $\mathbf{B}$ such that
$\mathbf{B}=\mathbf{P}\mathbf{B}\mathbf{P}$, where $\mathbf{P}$ is the
projection onto the column span of $\mathbf{A}$.

Based on these two structural conditions, we show the following result, which
we use to prove both Proposition 4 and Theorem 8. In this result, we will
refer to an $m\times n$ sketching matrix $\mathbf{S}_{m}$, indexed by the
number of rows $m$.

######  Theorem 11 (Structural conditions for near-unbiasedness).

Fix $\mathbf{A}\in\mathbb{R}^{n\times d}$ with rank $d$ and let
$\mathbf{S}_{m}$ consist of $m\geq 8d$ i.i.d. rows distributed as
$\frac{1}{\sqrt{m}}\mathbf{x}^{\scriptscriptstyle{\top}}$, where
$\mathbb{E}[\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}]=\mathbf{I}_{n}$.
Suppose that $\mathbf{S}_{m/3}$ satisfies Condition 1 (subspace embedding) for
$\eta=1/2$, with probability $1-\delta/3$, where $\delta\leq 1/m^{3}$. Suppose
also that $\mathbf{x}$ satisfies Condition 2 (Restricted Bai-Silverstein) with
some $\alpha\geq 1$. Then
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A})^{-1}$
is an $(\epsilon,\delta)$-unbiased estimator of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ for
$\epsilon=O(\alpha\sqrt{d}/m)$.

The proof of Theorem 11 adapts and extends techniques for analyzing the
limiting Stieltjes transform for high-dimensional random matrices in the so-
called Marchenko-Pastur regime (also called the proportional or mean-field
limit). This regime arises if we let $n$, $m$ and $d$ all go to infinity and
let the ratio $m/d$ converge to a fixed constant larger than unity. Crucially,
our analysis is non-asymptotic, and it is not restricted to the constant
aspect ratio between the sketch size and the dimension. Further, while
classical random matrix theory analysis considers matrix resolvents, which
take the form
$(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A}+z\mathbf{I})^{-1}$
for $z,\gamma\neq 0$, and are well-defined with full probability, we consider
the case of $z=0$ where the matrix in question may be undefined with positive
probability. We address this by defining a high probability event which
ensures that the sketch
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A})^{-1}$
is well-defined and bounded, while preserving enough of the independence
structure in the conditional distribution for the expectation analysis to go
through. Specifically, we split the sketch into three parts, and we condition
on the event that each part satisfies the subspace embedding property. This
way, for any pair of rows, there is a part of the sketch that ensures
invertibility while being independent from the two rows, which is important
for the analysis.

#### Subspace embedding condition.

Our first structural condition for small inversion bias (Condition 1) is a
variant of the subspace embedding property, which is standard in the sketching
literature. In particular, this condition immediately implies that
$(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A})^{-1}$
is with probability $1-\delta$ an $O(1)$-approximation of
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$. For sub-gaussian
sketches this is known to hold with sketch size $O(d+\log(1/\delta))$ [NN14].
We prove this for LESS of size $O(d\log(d/\delta))$.

######  Lemma 12 (Subspace embedding for LESS).

Suppose that $\mathbf{S}$ is a LESS embedding of size $m$ for a rank $d$
matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$. There is $C=O(1)$ such that if
$m\geq Cd\log(d/\delta)/{\eta}^{2}$ for ${\eta}\in(0,1)$, then the sketch
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A}$
is an $({\eta},\delta)$-approximation of
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$.

The subspace embedding guarantee for LESS embeddings is as good as that for
existing fast sketching methods. However, the analysis differs from the ones
used for either data-aware leverage score sampling or for data-oblivious
sparse sketches. We show the result by deriving a subexponential bound on the
matrix moments of a LESS embedding (Lemma 30), relying on a novel variant of
the Hanson-Wright concentration inequality for quadratic forms based on
orthogonal projection matrices (Lemma 31). We then use this to invoke a matrix
Bernstein inequality for random matrices with subexponential moments [Tro12,
Theorem 6.2].

#### Restricted Bai-Silverstein condition.

Our second structural condition for small inversion bias (Condition 2) is not
commonly seen in sketching, but we expect that it will be of broader interest
in adapting high-dimensional random matrix theory to RandNLA [DLLM20, DLM20,
DM21]. It is based on the classical inequality of Bai and Silverstein [BS10]
which bounds the deviation of a random quadratic form
$\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}$ from its mean. We
call it the _Restricted Bai-Silverstein condition_ because, unlike in the
classical version, we only require the inequality to hold for matrices
$\mathbf{B}$ that are restricted to the subspace spanned by the columns of
$\mathbf{A}$. By contrast, in classical random matrix theory it is often
assumed that the the following (unrestricted) condition holds.

######  Condition 3 (Bai-Silverstein).

Random vector $\mathbf{x}\in\mathbb{R}^{n}$ satisfies the (unrestricted) Bai-
Silverstein condition with $\alpha>0$, if
$\mathrm{Var}\\!\big{[}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}\big{]}\leq\alpha\cdot\mathrm{tr}(\mathbf{B}^{2})$
for all $n\times n$ p.s.d. matrices $\mathbf{B}$.

When the random vector $\mathbf{x}$ is $O(1)$-sub-gaussian, then Condition 3
is satisfied with $\alpha=O(1)$, as a consequence of the original inequality
of [BS10].333The original lemma applies more broadly to higher moments; we
cite only the case relevant to our analysis.

######  Lemma 13 (Bai-Silverstein inequality).

Let $\mathbf{x}$ have $n$ independent entries with mean zero and unit variance
such that $\mathbb{E}[x_{i}^{4}]=O(1)$. Then, Condition 3 is satisfied with
$\alpha=O(1)$.

##  5 Restricted Bai-Silverstein inequality

The Bai-Silverstein inequality from Lemma 13 does not directly apply to any of
the fast sketching methods discussed above (see Appendix C.3 for lower
bounds). However, we state and prove a generalization of this lemma, which
allows us to show the Restricted Bai-Silverstein condition (Condition 2) for
our new LESS embeddings.

To provide some intuition behind this result, consider the variance term
$\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}]$
which appears in the Restricted Bai-Silverstein condition, where
$\frac{1}{\sqrt{m}}\mathbf{x}^{\scriptscriptstyle{\top}}$ represents a random
row vector of the sketching matrix $\mathbf{S}$. The condition requires that
just this one row vector carries enough randomness to produce an accurate
sketch of the trace of a quadratic form $\mathbf{B}$. This is in contrast to
the subspace embedding condition, which uses the joint randomness of all the
rows of $\mathbf{S}$. Lemma 13 achieves this by enforcing a fourth-moment
bound on all of the entries of $\mathbf{x}$. Suppose that we sparsify this
vector, following the strategy of sparse embedding matrices, by multiplying
each entry of $\mathbf{x}$ with an independent scaled Bernoulli variable,
obtaining $\sqrt{\frac{m}{s}}\,b_{i}x_{i}$ for
$b_{i}\sim\mathrm{Bernoulli}(\frac{s}{m})$, where $s\ll m$ is the sparsity
level and $i$ is the entry index.444Most commonly studied sparse embedding
matrices have non-independent entries. However, the i.i.d. variant we consider
offers an equivalent guarantee for the subspace embedding property. See
[Coh16] and Appendix C.3. This preserves the mean and variance assumptions
from Lemma 13, but as long as $s=o(m)$, it violates the fourth-moment
assumption. Thus, it is natural to ask whether we can relax this fourth-moment
assumption. It turns out that, if we do the sparsification in a data-oblivious
manner, then the answer is no, since the random vector may not capture most of
the relevant directions in the matrix $\mathbf{B}$ (see Appendix C.3).
Importantly, this can occur even when the rows of the sketch _together_
capture all of the directions, ensuring the subspace embedding property, which
is already the case when we set the sparsity level to be as small as $s=O(\log
d)$. In other words, there is a wide gap between the sparsity needed to
preserve the Bai-Silverstein inequality, $s=\Omega(m)$, and sparsity needed to
ensure the subspace embedding.

Crucially, Theorem 11 does not require the Bai-Silverstein inequality to hold
for all $n\times n$ p.s.d. quadratic forms $\mathbf{B}$. Rather, it restricts
the family of quadratic forms to those that lie within the column-span of the
$n\times d$ data matrix $\mathbf{A}$. In particular, this restriction implies
that the matrix $\mathbf{B}$ is low-rank (it has at most rank $d$) and its
important directions are captured by the leverage scores of $\mathbf{A}$. We
take advantage of this additional information to relax the fourth-moment
assumptions, obtaining the following generalization of Lemma 13, which should
be of independent interest.

######  Theorem 14 (Restricted Bai-Silverstein inequality).

Fix a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ with rank $d$ and leverage
scores $l_{i}$, and let $\mathbf{x}$ have $n$ independent entries with mean
zero and unit variance such that $\mathbb{E}x_{i}^{4}\leq C/l_{i}$. Then,
$\mathbf{x}$ satisfies Condition 2 with $\alpha=C+2$ for matrix $\mathbf{A}$.

By setting $\mathbf{A}=\mathbf{I}_{n}$, where all leverage scores are 1 and
the restriction on $\mathbf{B}$ is vacuous, we not only recover the statement
of Lemma 13, but also our new analysis uses the Perron-Frobenius theorem to
obtain a tight constant factor in the bound (see Appendix C). However, when
$\mathbf{A}$ is a tall matrix, then the fourth-moment assumption becomes
potentially much more broadly applicable (for example, when the leverage
scores are uniform, we only need $\mathbb{E}x_{i}^{4}\leq C\cdot n/d$). In
particular, consider an i.i.d. sub-gaussian random vector $\mathbf{x}$
sparsified as follows: $\mathbf{x}\circ{\boldsymbol{\xi}}$, where we let
$\xi_{i}=b_{i}/\sqrt{l_{i}}$ and $b_{i}\sim\mathrm{Bernoulli}(l_{i})$. Then,
the entries satisfy the assumptions of Theorem 14, with expected number of
non-zeros equal to $d$. Note that this is different than the data-oblivious
sparsification discussed above, since the entries of the vector corresponding
to large leverage scores are less likely to be zeroed-out than others. This
form of sparsification is nearly equivalent to the one we use for our LESS
embeddings (see Definition 6; our analysis can be applied to either variant),
except that it leads to a non-deterministic level of sparsification. In
Appendix C.2 we prove the Restricted Bai-Silverstein condition with
$\alpha=O(1)$ for a leverage score sparsified vector constructed as in
Definition 6, which has non-independent entries.

##  6 Conclusions

We analyzed the phenomenon of inversion bias in sketching-based estimation
tasks involving the inverse covariance matrix. Inversion bias is a significant
bottleneck in methods that use parallel and distributed averaging. We showed
that certain classical sketching methods (such as sub-gaussian sketches) have
small inversion bias, while many algorithmically efficient sketches (such as
leverage score sampling) may not provide such a guarantee. Finally, we
developed a new efficient sketching method, called LEverage Score Sparsified
(LESS) embeddings, which has small inversion bias and its computational cost
is nearly-linear in the input size.

Estimation of the inverse covariance matrix and its various linear functionals
is motivated by a rich body of literature in statistics, data science,
numerical optimization, machine learning, signal processing, etc., which we
summarized in detail in Section 1.1. Here, we additionally remark that the
$(\epsilon,\delta)$-approximation guarantee we provide for the averaged
estimates of the inverse covariance (see Corollaries 5 and 9) immediately
implies corresponding approximation guarantees for linear functionals of the
inverse covariance in numerous tasks. In a distributed environment, one can
use this to build a system for querying such functionals, by aggregating
coarse estimates computed locally from $q$ sketches to produce an improved
global estimate with minimal communication cost. We illustrate this here for a
family of linear functionals of the form
$\mathrm{tr}\,\mathbf{C}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$,
parameterized by any p.s.d. matrix $\mathbf{C}$, as motivated by applications
in statistical inference (see Section 1.1). The claim follows from Corollary 9
by letting
$\mathbf{Q}_{i}=(\frac{m}{m-d}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}^{\scriptscriptstyle{\top}}\mathbf{S}_{i}\mathbf{A})^{-1}$.

######  Corollary 15 (Querying linear functionals).

For any matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and
$\epsilon,\delta\in(0,1)$, we can use LESS embeddings of size
$m=O(d\log(d/\epsilon\delta)+\sqrt{d}/\epsilon)$ to construct
$\mathbf{Q}_{1},...,\mathbf{Q}_{q}\in\mathbb{R}^{d\times d}$ in parallel time
$O({\mathrm{nnz}}(\mathbf{A})\log n+md^{2})$, where
$q=O(m\log^{2}(d/\delta))$, so that with probability $1-\delta$:

$\displaystyle\text{For all p.s.d. matrices $\mathbf{C}\in\mathbb{R}^{d\times
d}$,}\qquad\frac{1}{q}\sum_{i=1}^{q}\mathrm{tr}\,\mathbf{C}\mathbf{Q}_{i}\,\approx_{\epsilon}\,\mathrm{tr}\,\mathbf{C}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}.$

In the context of distributed optimization, our results can be directly
applied to show improved convergence guarantees, for instance, in the case of
the Distributed Iterative Hessian Sketch algorithm [PW16, DBPM20] and
Distributed Newton Sketch method [WRXM18, DM19]. Here, the quantity of
interest is of the form
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{b}$ for some
vector $\mathbf{b}$ (where $\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$
corresponds to the Hessian and $\mathbf{b}$ corresponds to the gradient). For
those methods, an $\epsilon$-approximation guarantee for the average of the
sketched inverse covariance matrices, as in Corollaries 5 and 9, directly
implies that the iterates $\mathbf{x}_{t}$ produced by the algorithms achieve
a convergence rate of the form $\Delta_{t}\leq
O(\epsilon^{t})\cdot\Delta_{0}$, where $\Delta_{t}$ represents distance from
the optimum in the $t$-th iteration. We illustrate this by applying Corollary
9 to the existing analysis of Distributed Newton Sketch, as outlined in
Section 4 of [DBPM20], obtaining an improved linear-quadratic convergence rate
for distributed empirical risk minimization.

######  Corollary 16 (Distributed Newton Sketch).

Consider a twice differentiable convex function of the form
$f(\mathbf{x})=\frac{1}{n}\sum_{i=1}^{n}\ell_{i}(\mathbf{x}^{\scriptscriptstyle{\top}}\phi_{i})+\frac{\lambda}{2}\|\mathbf{x}\|^{2}$,
where $\mathbf{x}\in\mathbb{R}^{d}$ and $\phi_{i}^{\scriptscriptstyle{\top}}$
is the $i$th row of an $n\times d$ data matrix $\Phi$. Given $\mathbf{x}_{t}$,
we can use LESS embeddings to construct $q$ independent randomized estimates
$\widehat{\mathbf{H}}_{1}(\mathbf{x}_{t}),...,\widehat{\mathbf{H}}_{q}(\mathbf{x}_{t})$
of the Hessian $\nabla^{2}f(\mathbf{x}_{t})$ in parallel time
$O({\mathrm{nnz}}(\Phi)\log n+md^{2})$, where
$m=O(d\log(d/\epsilon\delta)+\sqrt{d}/\epsilon)$ and
$q=O(m\log^{2}(d/\delta))$, so that

$\displaystyle\mathbf{x}_{t+1}$
$\displaystyle=\mathbf{x}_{t}-\frac{1}{q}\sum_{i=1}^{q}\widehat{\mathbf{H}}_{i}(\mathbf{x}_{t})^{-1}\nabla
f(\mathbf{x}_{t})\quad\text{with probability $1-\delta$ satisfies:}$
$\displaystyle\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|$
$\displaystyle\leq\max\big{\\{}\epsilon\cdot\sqrt{\kappa}\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|,\
\tfrac{2L}{\lambda_{\min}}\|\mathbf{x}_{t+1}-\mathbf{x}^{*}\|^{2}\big{\\}}\quad\text{for}\quad\mathbf{x}^{*}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{x}}f(\mathbf{x}),$

where $\kappa$, $L$, $\lambda_{\min}$ are the condition number, Lipschitz
constant and smallest eigenvalue of $\nabla^{2}f(\mathbf{x})$.

This result provides an improvement over the recently proposed DPP-based
sketching methods of [DBPM20], which suffer no inversion bias but are more
expensive, as well as over other fast sketching methods like row sampling
[WRXM18], which, as shown in this work, may indeed suffer from large inversion
bias.

### Acknowledgments

We would like to acknowledge DARPA, IARPA, NSF, and ONR via its BRC on RandNLA
for providing partial support of this work. Our conclusions do not necessarily
reflect the position or the policy of our sponsors, and no official
endorsement should be inferred.

## References

  * [AC06] Nir Ailon and Bernard Chazelle.  Approximate nearest neighbors and the fast johnson-lindenstrauss transform.  In Proceedings of the thirty-eighth annual ACM symposium on Theory of computing, pages 557–563. ACM, 2006.
  * [AC09] Nir Ailon and Bernard Chazelle.  The fast johnson–lindenstrauss transform and approximate nearest neighbors.  SIAM Journal on computing, 39(1):302–322, 2009.
  * [AD11] Alekh Agarwal and John C Duchi.  Distributed delayed stochastic optimization.  In J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 24, pages 873–881. Curran Associates, Inc., 2011.
  * [AG13] Naum Ilich Akhiezer and Izrail Markovich Glazman.  Theory of linear operators in Hilbert space.  Courier Corporation, 2013.
  * [AGZ10] Greg W Anderson, Alice Guionnet, and Ofer Zeitouni.  An Introduction to Random Matrices.  Number 118. Cambridge University Press, 2010.
  * [And03] Theodore W Anderson.  An Introduction to Multivariate Statistical Analysis.  Wiley New York, 2003.
  * [BBP17] Joël Bun, Jean-Philippe Bouchaud, and Marc Potters.  Cleaning large correlation matrices: tools from random matrix theory.  Physics Reports, 666:1–109, 2017.
  * [BC11] Heinz H Bauschke and Patrick L Combettes.  Convex analysis and monotone operator theory in Hilbert spaces, volume 408.  Springer, 2011.
  * [BD09] Peter J Brockwell and Richard A Davis.  Time series: theory and methods.  Springer, 2009.
  * [BS10] Zhidong Bai and Jack W Silverstein.  Spectral analysis of large dimensional random matrices, volume 20.  Springer, 2010.
  * [Bur73] Donald L Burkholder.  Distribution function inequalities for martingales.  the Annals of Probability, pages 19–42, 1973.
  * [BV04] Stephen Boyd and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [CCFC02] Moses Charikar, Kevin Chen, and Martin Farach-Colton.  Finding frequent items in data streams.  In International Colloquium on Automata, Languages, and Programming, pages 693–703. Springer, 2002.
  * [CD11] Romain Couillet and Merouane Debbah.  Random matrix methods for wireless communications.  Cambridge University Press, 2011.
  * [CDV20] Daniele Calandriello, Michał Dereziński, and Michal Valko.  Sampling from a k-dpp without looking at all items.  In Advances in Neural Information Processing Systems, volume 33, pages 6889–6899, 2020.
  * [CLL11] Tony Cai, Weidong Liu, and Xi Luo.  A constrained l1 minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association, 106(494):594–607, 2011.
  * [CLL+15] Shouyuan Chen, Yang Liu, Michael R Lyu, Irwin King, and Shengyu Zhang.  Fast relative-error approximation algorithm for ridge regression.  In UAI, pages 201–210, 2015.
  * [Coh16] Michael B Cohen.  Nearly tight oblivious subspace embeddings by trace inequalities.  In Proceedings of the twenty-seventh annual ACM-SIAM symposium on Discrete algorithms, pages 278–287. SIAM, 2016.
  * [CR00] David Roxbee Cox and Nancy Reid.  The theory of the design of experiments.  CRC Press, 2000.
  * [CS17] Timothy I Cannings and Richard J Samworth.  Random-projection ensemble classification.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(4):959–1035, 2017.
  * [CW17] Kenneth L. Clarkson and David P. Woodruff.  Low-rank approximation and regression in input sparsity time.  J. ACM, 63(6):54:1–54:45, January 2017.
  * [DBPM20] Michał Dereziński, Burak Bartan, Mert Pilanci, and Michael W Mahoney.  Debiasing distributed second order optimization with surrogate sketching and scaled regularization.  In Advances in Neural Information Processing Systems, volume 33, pages 6889–6899, 2020.
  * [DCMW19] Michał Dereziński, Kenneth L. Clarkson, Michael W. Mahoney, and Manfred K. Warmuth.  Minimax experimental design: Bridging the gap between statistical and worst-case approaches to least squares regression.  In Alina Beygelzimer and Daniel Hsu, editors, Proceedings of the Thirty-Second Conference on Learning Theory, volume 99 of Proceedings of Machine Learning Research, pages 1050–1069, Phoenix, USA, 25–28 Jun 2019\.
  * [DCV19] Michał Dereziński, Daniele Calandriello, and Michal Valko.  Exact sampling of determinantal point processes with sublinear time preprocessing.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 11542–11554. Curran Associates, Inc., 2019.
  * [Dem72] Arthur P Dempster.  Covariance selection.  Biometrics, pages 157–175, 1972.
  * [DKM06] P. Drineas, R. Kannan, and M. W. Mahoney.  Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix.  SIAM Journal on Computing, 36:158–183, 2006.
  * [DL18] Edgar Dobriban and Sifan Liu.  Asymptotic for sketching in least squares regression.  arXiv preprint arXiv:1810.06089, NeurIPS 2019, 2018.
  * [DLLM20] Michał Dereziński, Feynman Liang, Zhenyu Liao, and Michael W Mahoney.  Precise expressions for random projections: Low-rank approximation and randomized Newton.  In Advances in Neural Information Processing Systems, volume 33, pages 18272–18283, 2020.
  * [DLM20] Michał Dereziński, Feynman Liang, and Michael W Mahoney.  Exact expressions for double descent and implicit regularization via surrogate random design.  In Advances in Neural Information Processing Systems, volume 33, pages 5152–5164, 2020.
  * [DM16] Petros Drineas and Michael W. Mahoney.  RandNLA: Randomized numerical linear algebra.  Communications of the ACM, 59:80–90, 2016.
  * [DM17] Petros Drineas and Michael W Mahoney.  Lectures on randomized numerical linear algebra.  arXiv preprint arXiv:1712.08880, 2017.
  * [DM18] P. Drineas and M. W. Mahoney.  Lectures on randomized numerical linear algebra.  In M. W. Mahoney, J. C. Duchi, and A. C. Gilbert, editors, The Mathematics of Data, IAS/Park City Mathematics Series, pages 1–48. AMS/IAS/SIAM, 2018.
  * [DM19] Michał Dereziński and Michael W Mahoney.  Distributed estimation of the inverse Hessian by determinantal averaging.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 11401–11411. Curran Associates, Inc., 2019.
  * [DM21] Michał Dereziński and Michael W Mahoney.  Determinantal point processes in randomized numerical linear algebra.  Notices of the American Mathematical Society, 68(1):34–45, 2021\.
  * [DMIMW12] Petros Drineas, Malik Magdon-Ismail, Michael W. Mahoney, and David P. Woodruff.  Fast approximation of matrix coherence and statistical leverage.  J. Mach. Learn. Res., 13(1):3475–3506, December 2012.
  * [DMM06] Petros Drineas, Michael W Mahoney, and S Muthukrishnan.  Sampling algorithms for l 2 regression and applications.  In Proceedings of the seventeenth annual ACM-SIAM symposium on Discrete algorithm, pages 1127–1136. Society for Industrial and Applied Mathematics, 2006.
  * [DMM08] Petros Drineas, Michael W Mahoney, and Shan Muthukrishnan.  Relative-error CUR matrix decompositions.  SIAM Journal on Matrix Analysis and Applications, 30(2):844–881, 2008.
  * [DMMS11] Petros Drineas, Michael W Mahoney, S Muthukrishnan, and Tamás Sarlós.  Faster least squares approximation.  Numerische mathematik, 117(2):219–249, 2011.
  * [DS18] Edgar Dobriban and Yue Sheng.  Distributed linear regression by averaging.  arXiv preprint arxiv:1810.00412, to appear in The Annals of Statistics, 2018.
  * [DS20] Edgar Dobriban and Yue Sheng.  Wonder: Weighted one-shot distributed ridge regression in high dimensions.  Journal of Machine Learning Research, 21(66):1–52, 2020.
  * [DSBN15] G. Dasarathy, P. Shah, B. N. Bhaskar, and R. D. Nowak.  Sketching sparse matrices, covariances, and graphs via tensor products.  IEEE Transactions on Information Theory, 61(3):1373–1388, 2015\.
  * [DW18] Michał Dereziński and Manfred K. Warmuth.  Reverse iterative volume sampling for linear regression.  Journal of Machine Learning Research, 19(23):1–39, 2018.
  * [DWH18] Michał Dereziński, Manfred K. Warmuth, and Daniel Hsu.  Leveraged volume sampling for linear regression.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 2510–2519. Curran Associates, Inc., 2018.
  * [DWH19] Michał Dereziński, Manfred K. Warmuth, and Daniel Hsu.  Unbiased estimators for random design regression.  arXiv e-prints, page arXiv:1907.03411, Jul 2019.
  * [ER05] Alan Edelman and N Raj Rao.  Random matrix theory.  Acta numerica, 14:233, 2005.
  * [FHT08] Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9(3):432–441, 2008.
  * [FKV04] Alan Frieze, Ravi Kannan, and Santosh Vempala.  Fast monte-carlo algorithms for finding low-rank approximations.  Journal of the ACM (JACM), 51(6):1025–1041, 2004.
  * [GCS+13] AB Gelman, JB Carlin, HS Stern, DB Dunson, A Vehtari, and D Rubin.  Bayesian data analysis third edition. boca raton.  FL: CRC Press, 2013.
  * [GWS20] Milana Gataric, Tengyao Wang, and Richard J Samworth.  Sparse principal component analysis via axis-aligned random projections.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2020.
  * [Har69] JA Hartigan.  Linear bayesian methods.  Journal of the Royal Statistical Society. Series B (Methodological), pages 446–454, 1969.
  * [HMST11] Nathan Halko, Per-Gunnar Martinsson, Yoel Shkolnisky, and Mark Tygert.  An algorithm for the principal component analysis of large data sets.  SIAM Journal on Scientific computing, 33(5):2580–2594, 2011.
  * [HMT11] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM review, 53(2):217–288, 2011.
  * [HTF09] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The elements of statistical learning.  Springer series in statistics, 2009.
  * [KBRR16] Jakub Konecný, H. Brendan McMahan, Daniel Ramage, and Peter Richtárik.  Federated Optimization: Distributed Machine Learning for On-Device Intelligence.  arXiv e-prints, page arXiv:1610.02527, Oct 2016.
  * [KBY+16] Jakub Konecný, H. Brendan McMahan, Felix X. Yu, Peter Richtárik, Ananda Theertha Suresh, and Dave Bacon.  Federated Learning: Strategies for Improving Communication Efficiency.  arXiv e-prints, page arXiv:1610.05492, Oct 2016.
  * [LD19] Sifan Liu and Edgar Dobriban.  Ridge regression: Structure, cross-validation, and sketching.  arXiv preprint arXiv:1910.02373, International Conference on Learning Representations (ICLR) 2020, 2019.
  * [LDFU13] Yichao Lu, Paramveer Dhillon, Dean P Foster, and Lyle Ungar.  Faster ridge regression via the subsampled randomized hadamard transform.  In Advances in neural information processing systems, pages 369–377, 2013.
  * [LF09] Clifford Lam and Jianqing Fan.  Sparsistency and rates of convergence in large covariance matrix estimation.  Annals of statistics, 37(6B):4254, 2009.
  * [LJW11] Miles Lopes, Laurent Jacob, and Martin J Wainwright.  A more powerful two-sample test in high dimensions using random projection.  In Advances in Neural Information Processing Systems, pages 1206–1214, 2011.
  * [LW12] Olivier Ledoit and Michael Wolf.  Nonlinear shrinkage estimation of large-dimensional covariance matrices.  The Annals of Statistics, 40(2):1024–1060, 2012.
  * [LWM+07] Edo Liberty, Franco Woolfe, Per-Gunnar Martinsson, Vladimir Rokhlin, and Mark Tygert.  Randomized algorithms for the low-rank approximation of matrices.  Proceedings of the National Academy of Sciences, 104(51):20167–20172, 2007.
  * [Mah11] Michael W. Mahoney.  Randomized algorithms for matrices and data.  Foundations and Trends in Machine Learning, 3(2):123–224, 2011\.  Also available at: arXiv:1104.5557.
  * [MB06] Nicolai Meinshausen and Peter Bühlmann.  High-dimensional graphs and variable selection with the lasso.  The annals of statistics, 34(3):1436–1462, 2006.
  * [Mey00] Carl D Meyer.  Matrix analysis and applied linear algebra, volume 71.  Siam, 2000.
  * [MH15] Goran Marjanovic and Alfred O Hero.  l0 sparse inverse covariance estimation.  IEEE Transactions on Signal Processing, 63(12):3218–3231, 2015\.
  * [MHM10] Ryan McDonald, Keith Hall, and Gideon Mann.  Distributed training strategies for the structured perceptron.  In Human Language Technologies: The 2010 Annual Conference of the North American Chapter of the Association for Computational Linguistics, HLT ’10, pages 456–464, Stroudsburg, PA, USA, 2010. Association for Computational Linguistics.
  * [MM13] Xiangrui Meng and Michael W. Mahoney.  Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression.  In Proceedings of the Forty-fifth Annual ACM Symposium on Theory of Computing, STOC ’13, pages 91–100, New York, NY, USA, 2013. ACM.
  * [MM15] Cameron Musco and Christopher Musco.  Randomized block krylov methods for stronger and faster approximate singular value decomposition.  arXiv preprint arXiv:1504.05477, 2015.
  * [MM19] Song Mei and Andrea Montanari.  The generalization error of random features regression: Precise asymptotics and double descent curve.  arXiv preprint arXiv:1908.05355, 2019.
  * [MMS+09] Ryan McDonald, Mehryar Mohri, Nathan Silberman, Dan Walker, and Gideon S. Mann.  Efficient large-scale distributed training of conditional maximum entropy models.  In Y. Bengio, D. Schuurmans, J. D. Lafferty, C. K. I. Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems 22, pages 1231–1239. 2009.
  * [MMY15] P. Ma, M. W. Mahoney, and B. Yu.  A statistical perspective on algorithmic leveraging.  Journal of Machine Learning Research, 16:861–911, 2015.
  * [MP67] Vladimir A Marchenko and Leonid A Pastur.  Distribution of eigenvalues for some sets of random matrices.  Mat. Sb., 114(4):507–536, 1967.
  * [NN13] Jelani Nelson and Huy L. Nguyên.  OSNAP: Faster numerical linear algebra algorithms via sparser subspace embeddings.  In Proceedings of the 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, FOCS ’13, pages 117–126, Washington, DC, USA, 2013. IEEE Computer Society.
  * [NN14] Jelani Nelson and Huy L Nguyen.  Lower bounds for oblivious subspace embeddings.  In International Colloquium on Automata, Languages, and Programming, pages 883–894. Springer, 2014.
  * [NW06] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [Puk06] Friedrich Pukelsheim.  Optimal design of experiments.  SIAM, 2006.
  * [PW15] Mert Pilanci and Martin J Wainwright.  Randomized sketches of convex programs with sharp guarantees.  IEEE Transactions on Information Theory, 61(9):5096–5115, 2015\.
  * [PW16] Mert Pilanci and Martin J Wainwright.  Iterative hessian sketch: Fast and accurate solution approximation for constrained least-squares.  The Journal of Machine Learning Research, 17(1):1842–1879, 2016\.
  * [PW17] Mert Pilanci and Martin J Wainwright.  Newton sketch: A near linear-time optimization algorithm with linear-quadratic convergence.  SIAM Journal on Optimization, 27(1):205–245, 2017.
  * [PZ32] REAC Paley and Antoni Zygmund.  A note on analytic functions in the unit circle.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 28, pages 266–272. Cambridge University Press, 1932.
  * [Rao73] Calyampudi Radhakrishna Rao.  Linear statistical inference and its applications, volume 2.  Wiley New York, 1973.
  * [Ras03] Carl Edward Rasmussen.  Gaussian processes in machine learning.  In Summer School on Machine Learning, pages 63–71. Springer, 2003\.
  * [RKR+16] Sashank J. Reddi, Jakub Konecný, Peter Richtárik, Barnabás Póczós, and Alex Smola.  AIDE: Fast and Communication Efficient Distributed Optimization.  arXiv e-prints, page arXiv:1608.06879, Aug 2016.
  * [RM16] Garvesh Raskutti and Michael W Mahoney.  A statistical perspective on randomized sketching for ordinary least-squares.  The Journal of Machine Learning Research, 17(1):7508–7538, 2016\.
  * [RV13] Mark Rudelson and Roman Vershynin.  Hanson-wright inequality and sub-gaussian concentration.  Electronic Communications in Probability, 18, 2013.
  * [Sar06] Tamas Sarlos.  Improved approximation algorithms for large matrices via random projections.  In Foundations of Computer Science, 2006. FOCS’06. 47th Annual IEEE Symposium on, pages 143–152. IEEE, 2006.
  * [SLR16] Radhendushka Srivastava, Ping Li, and David Ruppert.  Raptt: An exact two-sample test in high dimensions using random projections.  Journal of Computational and Graphical Statistics, 25(3):954–970, 2016.
  * [SSZ14] Ohad Shamir, Nati Srebro, and Tong Zhang.  Communication-efficient distributed optimization using an approximate Newton-type method.  In Eric P. Xing and Tony Jebara, editors, Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 1000–1008, Bejing, China, 22–24 Jun 2014. PMLR.
  * [Tao12] Terence Tao.  Topics in random matrix theory, volume 132.  American Mathematical Soc., 2012.
  * [Tro12] Joel A Tropp.  User-friendly tail bounds for sums of random matrices.  Foundations of computational mathematics, 12(4):389–434, 2012.
  * [TYUC17] Joel A Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher.  Practical sketching algorithms for low-rank matrix approximation.  SIAM Journal on Matrix Analysis and Applications, 38(4):1454–1485, 2017.
  * [Vem05] Santosh S Vempala.  The random projection method, volume 65.  American Mathematical Soc., 2005.
  * [WB95] Greg Welch and Gary Bishop.  An introduction to the kalman filter, 1995.
  * [WGM18] Shusen Wang, Alex Gittens, and Michael W Mahoney.  Sketched ridge regression: Optimization perspective, statistical perspective, and model averaging.  Journal of Machine Learning Research, 18:1–50, 2018.
  * [Wil12] Virginia Vassilevska Williams.  Multiplying matrices faster than coppersmith-winograd.  In Proceedings of the forty-fourth annual ACM symposium on Theory of computing, pages 887–898, 2012.
  * [WLRT08] Franco Woolfe, Edo Liberty, Vladimir Rokhlin, and Mark Tygert.  A fast randomized algorithm for the approximation of matrices.  Applied and Computational Harmonic Analysis, 25(3):335–366, 2008\.
  * [Woo14] David P Woodruff.  Sketching as a tool for numerical linear algebra.  Foundations and Trends® in Theoretical Computer Science, 10(1–2):1–157, 2014.
  * [WRXM18] Shusen Wang, Fred Roosta, Peng Xu, and Michael W Mahoney.  GIANT: Globally improved approximate Newton method for distributed optimization.  In Advances in Neural Information Processing Systems, pages 2332–2342, 2018.
  * [YGKM19] Z. Yao, A. Gholami, K. Keutzer, and M. W. Mahoney.  PyHessian: Neural networks through the lens of the Hessian.  Technical Report Preprint: arXiv:1912.07145, 2019.
  * [YGS+20] Z. Yao, A. Gholami, S. Shen, K. Keutzer, and M. W. Mahoney.  ADAHESSIAN: An adaptive second order optimizer for machine learning.  Technical Report Preprint: arXiv:2006.00719, 2020.
  * [YLDW20] Fan Yang, Sifan Liu, Edgar Dobriban, and David P Woodruff.  How to reduce dimension with pca and random projections?  arXiv preprint arXiv:2005.00511, 2020.
  * [YZ05] Raphael Yuster and Uri Zwick.  Fast sparse matrix multiplication.  ACM Transactions On Algorithms (TALG), 1(1):2–13, 2005.
  * [ZDW13] Yuchen Zhang, John C. Duchi, and Martin J. Wainwright.  Communication-efficient algorithms for statistical optimization.  J. Mach. Learn. Res., 14(1):3321–3363, January 2013.
  * [ZL15] Yuchen Zhang and Xiao Lin.  Disco: Distributed optimization for self-concordant empirical loss.  In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 362–370, Lille, France, 07–09 Jul 2015\. PMLR.
  * [Zni09] Marko Znidaric.  Asymptotic expansion for inverse moments of binomial and poisson distributions.  The Open Mathematics, Statistics and Probability Journal, 1(1), 2009\.
  * [ZWLS10] Martin Zinkevich, Markus Weimer, Lihong Li, and Alex J. Smola.  Parallelized stochastic gradient descent.  In J. D. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, and A. Culotta, editors, Advances in Neural Information Processing Systems 23, pages 2595–2603. Curran Associates, Inc., 2010.

##  Appendix A Preliminaries

#### Notations.

In the remainder of the article, we follow the convention of denoting scalars
by lowercase, vectors by lowercase boldface, and matrices by uppercase
boldface letters. The norm $\|\cdot\|$ is the Euclidean norm for vectors and
the spectral or operator norm for matrices, and $\|\cdot\|_{F}$ is the
Frobenius norm for matrices. For vector $\mathbf{v}\in\mathbb{R}^{d}$, we let
$\|\mathbf{v}\|_{1}:=\sum_{i=1}^{d}|v_{i}|$ denote the $\ell_{1}$ norm and
$\|\mathbf{v}\|_{\infty}:=\max_{i}|v_{i}|$ denote the $\ell_{\infty}$ norm of
$\mathbf{v}$. We use $\lambda_{\max}(\mathbf{A})$ to denote the maximum
eigenvalue of a symmetric matrix $\mathbf{A}$. We say
$\mathbf{A}\preceq\mathbf{B}$ if and only if $\mathbf{B}-\mathbf{A}$ is
positive semi-definite. We use $\mathbf{A}\circ\mathbf{B}$ to denote the
entry-wise Hadamard product of matrices or vectors. For random vectors or
matrices, we say $\mathbf{A}\stackrel{{\scriptstyle d}}{{=}}\mathbf{B}$ if
$\mathbf{A}$ follows the same distribution as $\mathbf{B}$. For positive semi-
definite (p.s.d.) matrices $\mathbf{A}$ and $\mathbf{B}$, or non-negative
scalars $a$ and $b$, we use $\mathbf{A}\approx_{\eta}\mathbf{B}$ and
$a\approx_{\eta}b$ to denote the relative error approximation (Definition 1).
The big-O notation is used to absorb constant factors in upper bounds, where
the constant only depends on other big-O constants appearing in a given
statement (thus, all constants can be made absolute).

An important linear algebraic result that will be used in proving the
restricted Bai-Silverstein inequality (Theorem 14) is the following Perron-
Frobenius theorem on non-negative matrices. While the most well known version
of the Perron-Frobenius theorem concerns matrices with strictly positive
entries, there is also a version for matrices with only non-negative entries.

######  Lemma 17 (Perron-Frobenius theorem, [Mey00, claims 8.3.1 and 8.3.2]).

For a non-negative symmetric matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ such
that $[\mathbf{A}]_{ij}\geq 0$ for all $i,j\in\\{1,\ldots,n\\}$, then the
largest eigenvalue of $\mathbf{A}$ is non-negative, i.e.,
$r=\lambda_{\max}(\mathbf{A})\geq 0$. Moreover, there is a corresponding
eigenvector $\mathbf{z}$, i.e., $\mathbf{A}\mathbf{z}=r\mathbf{z}$ with non-
negative entries $\mathbf{z}_{i}\geq 0$ for all $i$.

Our analysis of the inversion bias (proof of Theorem 11) crucially relies on a
standard rank-one update formula for the matrix inverse, which is given below.

######  Lemma 18 (Sherman-Morrison formula).

For an invertible matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ and
$\mathbf{u},\mathbf{v}\in\mathbb{R}^{n}$,
$\mathbf{A}+\mathbf{u}\mathbf{v}^{\scriptscriptstyle{\top}}$ is invertible if
and only if
$1+\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{A}^{-1}\mathbf{u}\neq 0$. If
this holds, then

$(\mathbf{A}+\mathbf{u}\mathbf{v}^{\scriptscriptstyle{\top}})^{-1}=\mathbf{A}^{-1}-\frac{\mathbf{A}^{-1}\mathbf{u}\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{A}^{-1}}{1+\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{A}^{-1}\mathbf{u}}.$

In particular, it follows that:

$(\mathbf{A}+\mathbf{u}\mathbf{v}^{\scriptscriptstyle{\top}})^{-1}\mathbf{u}=\frac{\mathbf{A}^{-1}\mathbf{u}}{1+\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{A}^{-1}\mathbf{u}}.$

Our proofs rely on different types of concentration and anti-concentration
inequalities, from scalars to quadratic forms of the type
$\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}$, and eventually to
matrix concentration bounds. These technical lemmas are collected in this
section and will be repeatedly used in the proofs of our main results.

###  A.1 Scalar concentration and anti-concentration inequalities

The Burkholder inequality [Bur73] provides moment bounds on the sum of a
martingale difference sequence. It is used to show Lemma 25 as part of the
proof of Theorem 11.

######  Lemma 19 (Burkholder inequality, [Bur73]).

For $\\{x_{j}\\}_{j=1}^{m}$ a real martingale difference sequence with respect
to the increasing $\sigma$ field $\mathcal{F}_{j}$, we have, for $L>1$, there
exists $C_{L}>0$ such that

$\mathbb{E}\bigg{[}\Big{|}\sum_{j=1}^{m}x_{j}\Big{|}^{L}\bigg{]}\leq
C_{L}\cdot\mathbb{E}\bigg{[}\Big{(}\sum_{j=1}^{m}|x_{j}|^{2}\Big{)}^{L/2}\bigg{]}.$

The Paley-Zygmund inequality is used to establish an anti-concentration
inequality for the Binomial distribution (Lemma 36), which is the key in
deriving a lower bound for the inversion bias of leverage score sampling in
Appendix F.

######  Lemma 20 (Paley-Zygmund inequality, [PZ32]).

For any non-negative variable $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}]}.$

###  A.2 Quadratic form concentration

Being the key object of (one of) the structural conditions in Theorem 11, the
(random) quadratic form of the type
$\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}$ will consistently
appear in our analysis, for instance in the form of the Bai-Silverstein
inequality in Lemma 13 on quadratic form variance, as well as the following
Hanson-Wright inequality on the tail probability.

######  Lemma 21 (Hanson-Wright inequality, [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{)}.$

###  A.3 Matrix concentration inequalities

When random matrices are considered, different variants of Matrix
Chernoff/Bernstein inequalities will be needed to handle the case where the
random matrix under study is known to have (almost surely) bounded operator
norm, or only to admit a subexponential decay for its higher order moments.

######  Lemma 22 (Matrix Bernstein: Bounded Case, [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).$

######  Lemma 23 (Matrix Bernstein: Subexponential Case, [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).$

######  Lemma 24 (Matrix Chernoff, [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}.$

##  Appendix B Structural conditions for small inversion bias

In this section, we prove Theorem 11, which gives two structural conditions
for a random sketch of a rank $d$ matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$
to have small inversion bias. We assume that the sketching matrix
$\mathbf{S}_{m}\in\mathbb{R}^{m\times n}$ consists of $m\geq 8d$ i.i.d. rows
$\frac{1}{\sqrt{m}}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}$, where
$\mathbb{E}[\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}]=\mathbf{I}$.
To simplify the analysis, we assume that $m$ is divisible by $3$.

###  B.1 Proof of Theorem 11

Note that the subspace embedding assumption (based on Condition 1) immediately
implies the result with $\epsilon=O(1)$, so without loss of generality we can
assume that $\alpha\sqrt{d}/m\leq 1$. Let
$\mathbf{H}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$ and
$\mathbf{Q}=(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A})^{-1}$
for $\gamma=\frac{m}{m-d}$. Moreover, let $\mathbf{S}_{-i}$ denote
$\mathbf{S}_{m}$ without the $i$th row, with
$\mathbf{Q}_{-i}=(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{A})^{-1}$.
Finally, for $t=m/3$, we define the following events:

$\displaystyle\mathcal{E}_{j}:\
\frac{1}{t}\mathbf{A}^{\scriptscriptstyle{\top}}\Big{(}\sum_{i=t(j-1)+1}^{tj}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\Big{)}\mathbf{A}\succeq\frac{1}{2}\cdot\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A},\quad
j=1,2,3,\qquad\mathcal{E}=\bigwedge_{j=1}^{3}\mathcal{E}_{j}.$ (2)

For each $j$, the meaning of the event $\mathcal{E}_{j}$ is that the average
of the rank one matrices
$\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}$ over the
corresponding $j$-th third of indices $1,\ldots,m$ forms a sketch for
$\mathbf{A}$ that is a “lower” spectral approximation of
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$.

Note that events $\mathcal{E}_{1}$, $\mathcal{E}_{2}$ and $\mathcal{E}_{3}$
are independent, and for each $i\in\\{1,...,m\\}$ there is a
$j=j(i)\in\\{1,2,3\\}$ such that:

  1. 1.

$\mathcal{E}_{j}$ is independent of $\mathbf{x}_{i}$; and

  2. 2.

$\mathcal{E}_{j}$ implies that
$\mathbf{Q}_{-i}\preceq\gamma\mathbf{Q}_{-i}=(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{A})^{-1}\preceq
6\cdot(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}=6\cdot\mathbf{H}^{-1}$.

Here we use that
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A}$
is the average of the three matrices to which the conditions in
$\mathcal{E}_{j}$ refer to, and also that $m\geq 2d$.

From the subspace embedding assumption and the union bound we conclude that
$\mathrm{Pr}(\mathcal{E})\geq 1-\delta$. Letting $\mathbb{E}_{\mathcal{E}}$
denote the expectation conditioned on $\mathcal{E}$ and
$\gamma_{i}=1+\frac{\gamma}{m}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}$,
we have:

$\displaystyle\mathbf{I}-\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H}$
$\displaystyle=-\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H}+\gamma\,\mathbb{E}_{\mathcal{E}}[\mathbf{Q}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}^{\scriptscriptstyle{\top}}\mathbf{S}_{m}\mathbf{A}]=-\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H}+\gamma\,\mathbb{E}_{\mathcal{E}}[\mathbf{Q}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}]$
$\displaystyle\overset{(*)}{=}-\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H}+\gamma\,\mathbb{E}_{\mathcal{E}}\Big{[}\tfrac{\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}}{1+\frac{\gamma}{m}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}}\Big{]}$
$\displaystyle=-\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H}+\mathbb{E}_{\mathcal{E}}[\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}]+\mathbb{E}_{\mathcal{E}}\big{[}(\tfrac{\gamma}{\gamma_{i}}-1)\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\big{]}$
$\displaystyle=\underbrace{\mathbb{E}_{\mathcal{E}}[\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}]}_{\mathbf{Z}_{0}}+\underbrace{\mathbb{E}_{\mathcal{E}}[\mathbf{Q}_{-i}-\mathbf{Q}]\mathbf{H}}_{\mathbf{Z}_{1}}+\underbrace{\mathbb{E}_{\mathcal{E}}\big{[}(\tfrac{\gamma}{\gamma_{i}}-1)\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\big{]}}_{\mathbf{Z}_{2}},$

for a fixed $i$, where $(*)$ uses the Sherman-Morrison rank-one update formula
(Lemma 18). We also used the fact that due to symmetry in the definition of
event $\mathcal{E}$, the marginal distributions of the random vectors
$\mathbf{x}_{i}$ are identical after conditioning (even though they are no
longer independent and identically distributed). To obtain the result, it
suffices to bound:

$\displaystyle\|\mathbf{I}-\mathbf{H}^{\frac{1}{2}}\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H}^{\frac{1}{2}}\|$
$\displaystyle=\|\mathbf{H}^{\frac{1}{2}}(\mathbf{Z}_{0}+\mathbf{Z}_{1}+\mathbf{Z}_{2})\mathbf{H}^{-\frac{1}{2}}\|$
$\displaystyle\leq\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{0}\mathbf{H}^{-\frac{1}{2}}\|+\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{1}\mathbf{H}^{-\frac{1}{2}}\|+\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{2}\mathbf{H}^{-\frac{1}{2}}\|.$
(3)

We start by bounding the first term. Without loss of generality, assume that
events $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ are both independent of
$\mathbf{x}_{i}$, and let
$\mathcal{E}^{\prime}=\mathcal{E}_{1}\wedge\mathcal{E}_{2}$ as well as
$\delta_{3}=\mathrm{Pr}(\neg\mathcal{E}_{3})$. We have:

$\displaystyle\mathbf{Z}_{0}$
$\displaystyle=\frac{1}{1-\delta_{3}}\cdot\Big{(}\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}]-\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}\cdot\mathbf{1}_{\neg\mathcal{E}_{3}}]\Big{)}$
$\displaystyle=-\frac{1}{1-\delta_{3}}\cdot\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}\cdot\mathbf{1}_{\neg\mathcal{E}_{3}}\big{]}.$

Above, we evaluated the expectation
$\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}]$
by first conditioning on all randomness except $\mathbf{x}_{i}$, and using the
independence of $\mathbf{x}_{i}$ and $\mathcal{E}^{\prime}$, as well as
$\mathbb{E}[\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}]=\mathbf{I}$.

Thus, since $\delta_{3}\leq\frac{1}{2}$, we obtain that:

$\displaystyle\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{0}\mathbf{H}^{-\frac{1}{2}}\|$
$\displaystyle\leq
2\,\Big{\|}\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}\mathbf{H}^{-\frac{1}{2}}\cdot\mathbf{1}_{\neg\mathcal{E}_{3}}\big{]}\Big{\|}$
$\displaystyle\leq
2\,\mathbb{E}_{\mathcal{E}^{\prime}}\Big{[}\big{\|}\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}\mathbf{H}^{-\frac{1}{2}}\big{\|}\cdot\mathbf{1}_{\neg\mathcal{E}_{3}}\Big{]}$
$\displaystyle\leq
2\,\mathbb{E}_{\mathcal{E}^{\prime}}\Big{[}\|\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{H}^{\frac{1}{2}}\|\cdot\big{\|}\mathbf{H}^{-\frac{1}{2}}\mathbf{A}^{\scriptscriptstyle{\top}}(\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I})\mathbf{A}\mathbf{H}^{-\frac{1}{2}}\big{\|}\cdot\mathbf{1}_{\neg\mathcal{E}_{3}}\Big{]}$
$\displaystyle\leq
12\,\mathbb{E}_{\mathcal{E}^{\prime}}\Big{[}\big{(}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}+1\big{)}\cdot\mathbf{1}_{\neg\mathcal{E}_{3}}\Big{]}.$

Note that
$\mathbb{E}[\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}]=d$,
and using Condition 2 (Restricted Bai-Silverstein), we have
$\mathrm{Var}[\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}]$
$\leq\alpha\cdot d$ (and both are still true after conditioning on
$\mathcal{E}^{\prime}$, because it is independent of $\mathbf{x}_{i}$).
Chebyshev’s inequality thus implies that for $x\geq 2d$ we have
$\mathrm{Pr}(\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\geq
x\mid\mathcal{E}^{\prime})\leq C\alpha d/x^{2}$. Combining this with the
assumption that $\delta_{3}\leq\delta\leq 1/m^{3}$, we have:

$\displaystyle\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\cdot\mathbf{1}_{\neg\mathcal{E}}\big{]}$
$\displaystyle=\int_{0}^{\infty}\mathrm{Pr}(\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\cdot\mathbf{1}_{\neg\mathcal{E}}\geq
x\mid\mathcal{E}^{\prime})\,dx$ $\displaystyle\leq
2m^{2}\delta_{3}+\int_{2m^{2}}^{\infty}\mathrm{Pr}(\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\geq
x)\,dx$ $\displaystyle\leq\frac{2}{m}+C\alpha
d\int_{2m^{2}}^{\infty}\frac{1}{x^{2}}\,dx\leq\frac{2}{m}+C\frac{\alpha
d}{m^{2}},$

which implies that
$\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{0}\mathbf{H}^{-\frac{1}{2}}\|=O(1/m+\alpha
d/m^{2})=O(\alpha\sqrt{d}/m)$. We now move on to bounding the second term in
(3). In the following, we will use the observation that for a p.s.d. random
matrix $\mathbf{C}$ (or non-negative random variable) in the probability space
of $\mathbf{S}_{m}$, we have:

$\displaystyle\mathbb{E}_{\mathcal{E}}[\mathbf{C}]=\frac{\mathbb{E}[(\prod_{j=1}^{3}\mathbf{1}_{\mathcal{E}_{j}})\cdot\mathbf{C}]}{\mathrm{Pr}(\mathcal{E})}\preceq\frac{1}{1-\delta}\,\mathbb{E}[\mathbf{1}_{\mathcal{E}^{\prime}}\cdot\mathbf{C}]\preceq
2\cdot\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{C}].$ (4)

Using the above, and the fact that event $\mathcal{E}^{\prime}$ is independent
of $\mathbf{x}_{i}$, we have:

$\displaystyle\mathbb{E}_{\mathcal{E}}[\mathbf{Q}_{-i}-\mathbf{Q}]$
$\displaystyle\preceq
2\cdot\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}-\mathbf{Q}]=\frac{2\gamma}{m}\cdot\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}\gamma_{i}^{-1}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\big{]}\preceq\frac{2\gamma}{m}\cdot\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}\mathbf{H}\mathbf{Q}_{-i}].$

We now bound the second term in (3) by using the fact that
$\mathcal{E}^{\prime}$ implies
$\mathbf{H}^{\frac{1}{2}}\gamma\mathbf{Q}_{-i}\mathbf{H}^{\frac{1}{2}}\preceq
6\mathbf{I}$:

$\displaystyle\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{1}\mathbf{H}^{-\frac{1}{2}}\|=\|\mathbf{H}^{\frac{1}{2}}\mathbb{E}_{\mathcal{E}}[\mathbf{Q}_{-i}-\mathbf{Q}]\mathbf{H}^{\frac{1}{2}}\|\leq\frac{2\gamma}{m}\cdot\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}\|\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{H}^{\frac{1}{2}}\cdot\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{H}^{\frac{1}{2}}\|\big{]}\leq\frac{2}{m}\cdot
36.$

We next bound the last term in (3), applying the Cauchy-Schwarz inequality
twice:

$\displaystyle\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{2}\mathbf{H}^{-\frac{1}{2}}\|$
$\displaystyle=\sup_{\|\mathbf{v}\|=1,~{}\|\mathbf{u}\|=1}\mathbb{E}_{\mathcal{E}}\Big{[}\left|\tfrac{\gamma}{\gamma_{i}}-1\right|\cdot\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-\frac{1}{2}}\mathbf{u}\Big{]}$
$\displaystyle\leq\sqrt{\mathbb{E}_{\mathcal{E}}\big{[}(\tfrac{\gamma}{\gamma_{i}}-1)^{2}\big{]}}\cdot\sup_{\|\mathbf{v}\|=1,~{}\|\mathbf{u}\|=1}\sqrt{\mathbb{E}_{\mathcal{E}}\big{[}(\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\cdot\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-\frac{1}{2}}\mathbf{u})^{2}\big{]}}$
$\displaystyle\leq\underbrace{\sqrt{\mathbb{E}_{\mathcal{E}}\big{[}(\gamma_{i}-\gamma)^{2}\big{]}}}_{T_{1}}\cdot\underbrace{\sup_{\|\mathbf{u}\|=1}\sqrt[4]{\mathbb{E}_{\mathcal{E}}\big{[}(\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i})^{4}\big{]}}}_{T_{2}}\cdot\underbrace{\sup_{\|\mathbf{u}\|=1}\sqrt[4]{\mathbb{E}_{\mathcal{E}}\big{[}(\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{-\frac{1}{2}}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i})^{4}\big{]}}}_{T_{3}}.$

To bound $T_{3}$, we rely on Restricted Bai-Silverstein with
$\mathbf{B}=\mathbf{A}\mathbf{H}^{-\frac{1}{2}}\mathbf{u}\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{-\frac{1}{2}}\mathbf{A}^{\scriptscriptstyle{\top}}$,
noting that
$\mathrm{tr}(\mathbf{B}^{2})=\mathrm{tr}(\mathbf{B})=(\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{-\frac{1}{2}}\mathbf{H}\mathbf{H}^{-\frac{1}{2}}\mathbf{u})^{2}=\|\mathbf{u}\|^{4}=1$.
Recall that event $\mathcal{E}^{\prime}$ is independent of $\mathbf{x}_{i}$,
so we have:

$\displaystyle\mathbb{E}_{\mathcal{E}}\big{[}(\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{-\frac{1}{2}}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i})^{4}\big{]}$
$\displaystyle\leq
2\,\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}(\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{-\frac{1}{2}}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i})^{4}\big{]}$
$\displaystyle=2\,\mathbb{E}\big{[}(\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}_{i})^{2}\big{]}$
$\displaystyle=2\mathrm{Var}[\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}_{i}]+2\big{(}\mathbb{E}[\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}_{i}]\big{)}^{2}$
$\displaystyle\leq
2\alpha\cdot\mathrm{tr}(\mathbf{B}^{2})+2\big{(}\mathrm{tr}(\mathbf{B})\big{)}^{2}=2(\alpha+1),$

obtaining that $T_{3}=O(\sqrt[4]{\alpha+1})$. We can similarly bound $T_{2}$
by letting
$\mathbf{B}=\mathbf{A}\mathbf{Q}_{-i}\mathbf{H}^{\frac{1}{2}}\mathbf{u}\mathbf{u}^{\scriptscriptstyle{\top}}\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}$.
Note that, conditioned on $\mathcal{E}^{\prime}$, we again have

$\mathrm{tr}(\mathbf{B}^{2})=\big{(}\mathbf{u}^{\scriptscriptstyle{\top}}(\mathbf{H}^{\frac{1}{2}}\mathbf{Q}_{-i}\mathbf{H}^{\frac{1}{2}})^{2}\mathbf{u}\big{)}^{2}\leq
6^{4},$

so analogously as above we conclude that $T_{2}=O(\sqrt[4]{\alpha+1})$.

It thus remains to bound the term $T_{1}$. First, note that:

$\displaystyle\mathbb{E}_{\mathcal{E}}[(\gamma-\gamma_{i})^{2}]\leq
2\,\mathbb{E}_{\mathcal{E}^{\prime}}[(\gamma-\gamma_{i})^{2}]=2\,(\gamma-\bar{\gamma})^{2}+2\,\mathbb{E}_{\mathcal{E}^{\prime}}[(\gamma_{i}-\bar{\gamma})^{2}],$
(5)

where
$\bar{\gamma}=\mathbb{E}_{\mathcal{E}^{\prime}}[\gamma_{i}]=1+\frac{\gamma}{m}\mathrm{tr}(\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}]\mathbf{H})$.
To bound the second term in (5), we write:

$\displaystyle\mathbb{E}_{\mathcal{E}^{\prime}}[(\gamma_{i}-\bar{\gamma})^{2}]$
$\displaystyle=\frac{\gamma^{2}}{m^{2}}\,\mathbb{E}_{\mathcal{E}^{\prime}}\Big{[}\big{(}\mathrm{tr}(\mathbf{Q}_{-i}-\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}])\mathbf{H}\big{)}^{2}\Big{]}+\frac{\gamma^{2}}{m^{2}}\,\mathbb{E}_{\mathcal{E}^{\prime}}\Big{[}\big{(}\mathrm{tr}(\mathbf{Q}_{-i}\mathbf{H})-\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{)}^{2}\Big{]}.$

The latter term can be bounded by again using Condition 2, with
$\mathbf{B}=\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}$,
obtaining:

$\displaystyle\frac{\gamma^{2}}{m^{2}}\,\mathbb{E}_{\mathcal{E}^{\prime}}\Big{[}\big{(}\mathrm{tr}(\mathbf{Q}_{-i}\mathbf{H})-\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{)}^{2}\Big{]}$
$\displaystyle\leq\frac{1}{m^{2}}\,\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}\alpha\cdot\mathrm{tr}((\gamma\mathbf{Q}_{-i}\mathbf{H})^{2})\big{]}\leq
36\cdot\frac{\alpha d}{m^{2}}.$

The former term can be bounded using the Burkholder inequality for martingale
difference sequences. We state this bound as a lemma, proven separately in
Appendix B.2.

######  Lemma 25.

Let $\mathrm{Var}_{\mathcal{E}^{\prime}}[\cdot]$ be the conditional variance
with respect to event
$\mathcal{E}^{\prime}=\mathcal{E}_{1}\wedge\mathcal{E}_{2}$, see (2), with
$\mathbf{x}_{i}$ independent of $\mathcal{E}^{\prime}$. Then, there is an
absolute constant $C>0$ such that:

$\displaystyle\mathrm{Var}_{\mathcal{E}^{\prime}}\\!\big{[}\mathrm{tr}(\mathbf{Q}_{-i}\mathbf{H})\big{]}\leq
C\cdot d.$

Using Lemma 25, we conclude that
$\mathbb{E}_{\mathcal{E}^{\prime}}[(\gamma_{i}-\bar{\gamma})^{2}]\leq
C^{\prime}\cdot\alpha d/m^{2}$ for some absolute constant $C^{\prime}$. It
remains to bound the term:

$\displaystyle|\gamma-\bar{\gamma}|=\bigg{|}\frac{m}{m-d}-\Big{(}1+\frac{\gamma}{m}\mathrm{tr}(\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}]\mathbf{H})\Big{)}\bigg{|}=\frac{|d-\mathrm{tr}(\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}]\mathbf{H})|}{m-d}.$

Observe that we have:

$\displaystyle\big{|}d-\mathrm{tr}\,\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}]\mathbf{H}\big{|}$
$\displaystyle=\big{|}\mathrm{tr}((\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]-\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}])\mathbf{H})+\mathrm{tr}(\mathbf{I}-\mathbb{E}_{\mathcal{E}}[\mathbf{Q}]\mathbf{H})\,\big{|}$
$\displaystyle=\big{|}\mathrm{tr}((\mathbb{E}_{\mathcal{E}}-\mathbb{E}_{\mathcal{E}^{\prime}})[\mathbf{Q}_{-i}]\mathbf{H})+\mathrm{tr}(-\mathbf{Z}_{1})+\mathrm{tr}(\mathbf{Z}_{0}+\mathbf{Z}_{1}+\mathbf{Z}_{2})\big{|}$
$\displaystyle\leq\big{|}\mathrm{tr}((\mathbb{E}_{\mathcal{E}}-\mathbb{E}_{\mathcal{E}^{\prime}})[\mathbf{Q}_{-i}]\mathbf{H})\big{|}+|\mathrm{tr}(\mathbf{Z}_{0})|+|\mathrm{tr}(\mathbf{Z}_{2})|.$

The first two terms can be bounded similarly as we did
$\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{0}\mathbf{H}^{-\frac{1}{2}}\|$,
obtaining that $|\mathrm{tr}(\mathbf{Z}_{0})|=O(\alpha d/m)$, and also:

$\displaystyle\big{|}\mathrm{tr}((\mathbb{E}_{\mathcal{E}}-\mathbb{E}_{\mathcal{E}^{\prime}})[\mathbf{Q}_{-i}]\mathbf{H})\big{|}=\frac{\delta_{3}}{1-\delta_{3}}\big{|}\mathrm{tr}((\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}]-\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}\mid\neg\mathcal{E}_{3}])\mathbf{H})\big{|}=O(d\delta_{3})=O(d/m^{3}).$

For the last term, we have:

$\displaystyle|\mathrm{tr}(\mathbf{Z}_{2})|$
$\displaystyle=\Big{|}\mathbb{E}_{\mathcal{E}}\big{[}(\tfrac{\gamma}{\gamma_{i}}-1)\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{]}\Big{|}$
$\displaystyle\leq\Big{|}\mathbb{E}_{\mathcal{E}}\big{[}\tfrac{\gamma-\bar{\gamma}}{\gamma_{i}}\,\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{]}\Big{|}+\Big{|}\mathbb{E}_{\mathcal{E}}\big{[}\tfrac{\bar{\gamma}-\gamma_{i}}{\gamma_{i}}\,\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{]}\Big{|}$
$\displaystyle\leq|\gamma-\bar{\gamma}|\cdot\mathbb{E}_{\mathcal{E}}[\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}]+(m-d)\cdot\mathbb{E}\big{[}|\gamma_{i}-\bar{\gamma}|\big{]}$
$\displaystyle\leq|\gamma-\bar{\gamma}|\cdot\frac{6}{1-\delta}\,\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{H}^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}]+(m-d)\cdot\sqrt{\mathbb{E}[(\gamma_{i}-\bar{\gamma})^{2}]}$
$\displaystyle\leq|\gamma-\bar{\gamma}|\cdot\frac{6}{1-\delta}\,d+\sqrt{C^{\prime}\alpha
d}.$

The bound for the second term
$\Big{|}\mathbb{E}_{\mathcal{E}}\big{[}\tfrac{\bar{\gamma}-\gamma_{i}}{\gamma_{i}}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{]}\Big{|}$
comes from the definition of
$\gamma_{i}=1+\frac{\gamma}{m}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}$,
because
$\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{/}\gamma_{i}\leq
m/\gamma=m-d$.

Thus, putting this together we conclude that:

$\displaystyle|\gamma-\bar{\gamma}|\cdot\big{(}1-\tfrac{6d}{(m-d)(1-\delta)}\big{)}\leq
O(\alpha d/m^{2})+O(\sqrt{\alpha d}/m)=O(\sqrt{\alpha d}/m),$

which for $m\geq 8d$ and $\delta\leq 1/m^{3}$ implies that
$(\gamma-\bar{\gamma})^{2}=O(\alpha d/m^{2})$ so we get $T_{1}=O(\sqrt{\alpha
d}/m)$. Finally, we obtain the bound
$\|\mathbf{H}^{\frac{1}{2}}\mathbf{Z}_{2}\mathbf{H}^{-\frac{1}{2}}\|\leq
T_{1}\cdot T_{2}\cdot T_{3}=O(\alpha\sqrt{d}/m)$, which concludes the proof.

###  B.2 Proof of Lemma 25

Let $\mathbf{Q}_{-ij}$ denote the matrix
$(\gamma\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-ij}^{\scriptscriptstyle{\top}}\mathbf{S}_{-ij}\mathbf{A})^{-1}$
where $\mathbf{S}_{-ij}$ is the matrix $\mathbf{S}_{m}$ without the $i$th and
$j$th rows and $\gamma=\frac{m}{m-d}$. Let
$\mathbb{E}_{\mathcal{E}^{\prime},j}[\cdot]$ be the conditional expectation
with respect to $\mathcal{E}^{\prime}$ and the $\sigma$-field
$\mathcal{F}_{j}$ generating the rows
$\frac{1}{\sqrt{m}}\mathbf{x}_{1}^{\scriptscriptstyle{\top}}\ldots,\frac{1}{\sqrt{m}}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}$
of $\mathbf{S}$. First note that

$\displaystyle\mathrm{tr}(\mathbf{Q}_{-i}-\mathbb{E}_{\mathcal{E}^{\prime}}\mathbf{Q}_{-i})\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$
$\displaystyle=\mathbb{E}_{\mathcal{E}^{\prime},m}[\mathrm{tr}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]-\mathbb{E}_{\mathcal{E}^{\prime},0}[\mathrm{tr}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]$
$\displaystyle=\sum_{j=1}^{m}\left(\mathbb{E}_{\mathcal{E}^{\prime},j}[\mathrm{tr}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]-\mathbb{E}_{\mathcal{E}^{\prime},j-1}[\mathrm{tr}\mathbf{Q}_{-i}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]\right)=-\sum_{j=1}^{m}(\psi_{j}+\xi_{j}),$
$\displaystyle\text{where}\quad\psi_{j}$
$\displaystyle=(\mathbb{E}_{\mathcal{E}^{\prime},j}-\mathbb{E}_{\mathcal{E}^{\prime},j-1})[\mathrm{tr}(\mathbf{Q}_{-ij}-\mathbf{Q}_{-i})\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]$
$\displaystyle\text{ and }\ \xi_{j}$
$\displaystyle=-(\mathbb{E}_{\mathcal{E}^{\prime},j}-\mathbb{E}_{\mathcal{E}^{\prime},j-1})[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}].$

This forms a martingale difference sequence and hence falls within the scope
of the Burkholder inequality [Bur73], recalled as follows. We mention that
similar martingale decomposition techniques are common in random matrix
theory, see e.g., [BS10]. Also, for the case $L=2$ that we will use,
Burkholder inequality is nothing but the law of iterated variance.

######  Lemma 26 ([Bur73]).

For $\\{x_{j}\\}_{j=1}^{m}$ a real martingale difference sequence with respect
to the increasing $\sigma$ field $\mathcal{F}_{j}$, we have, for $L>1$, there
exists $C_{L}>0$ such that

$\mathbb{E}\bigg{[}\Big{|}\sum_{j=1}^{m}x_{j}\Big{|}^{L}\bigg{]}\leq
C_{L}\cdot\mathbb{E}\bigg{[}\Big{(}\sum_{j=1}^{m}|x_{j}|^{2}\Big{)}^{L/2}\bigg{]}.$

Note that for each pair $i,j$, one of $\mathcal{E}_{1},\mathcal{E}_{2}$ is
independent of both $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$. Without loss of
generality, suppose that this is $\mathcal{E}_{1}$. Then, in particular,
$\mathcal{E}_{1}$ implies that
$\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\preceq
6\,\mathbf{I}$. Thus, conditioned on $\mathcal{E}_{1}$, it follows that

$\displaystyle\mathrm{tr}(\mathbf{Q}_{-ij}-\mathbf{Q}_{-i})\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$
$\displaystyle=\mathrm{tr}\bigg{(}\frac{\frac{\gamma}{m}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}}{1+\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}\bigg{)}$
$\displaystyle=\frac{\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}(\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}})^{2}\mathbf{x}_{j}}{1+\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}}\leq\frac{6\cdot\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}}{1+\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}}\leq
6,$

which implies that $|\psi_{j}|\leq 6$. We now provide a bound on the second
moment of $\psi_{j}$, bounding the $\mathcal{E}^{\prime}$-conditional
expectation in terms of the $\mathcal{E}_{1}$-conditional expectation
analogously as in (4):

$\displaystyle\mathbb{E}_{\mathcal{E}^{\prime}}[\psi_{j}^{2}]$
$\displaystyle\leq
2\cdot\mathbb{E}_{\mathcal{E}_{1}}\Bigg{[}\bigg{(}\frac{6\cdot\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}}{1+\frac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}}\bigg{)}^{2}\Bigg{]}\leq
72\cdot\mathbb{E}_{\mathcal{E}_{1}}[\tfrac{\gamma}{m}\mathbf{x}_{j}^{\scriptscriptstyle{\top}}\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{x}_{j}]$
$\displaystyle=72\cdot\frac{\mathbb{E}_{\mathcal{E}_{1}}[\mathrm{tr}\,\mathbf{A}\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}]}{m-d}\leq
72\cdot 6\cdot\frac{d}{m-d}.$

We now aim to bound $|\xi_{j}|$. Since $\mathcal{E}_{1}$ is independent of
$\mathbf{x}_{j}$, we have
$\mathbb{E}_{\mathcal{E}_{1},j}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]=\mathbb{E}_{\mathcal{E}_{1},j-1}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]$.
Furthermore, letting $\delta_{2}=\mathrm{Pr}(\neg\mathcal{E}_{2})$, we have:

$\displaystyle\mathbb{E}_{\mathcal{E}_{1},j-1}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]$
$\displaystyle=\mathbb{E}_{\mathcal{E}^{\prime},j-1}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}](1-\delta_{2})+\mathbb{E}_{\mathcal{E}_{1},j-1}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}\mid\neg\mathcal{E}_{2}]\delta_{2},$
$\displaystyle\mathbb{E}_{\mathcal{E}_{1},j}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]$
$\displaystyle=\mathbb{E}_{\mathcal{E}^{\prime},j}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}](1-\delta_{2})+\mathbb{E}_{\mathcal{E}_{1},j}[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}\mid\neg\mathcal{E}_{2}]\delta_{2}.$

Thus, subtracting the two equalities from each other, we conclude that:

$\displaystyle|\xi_{j}|$
$\displaystyle=|(\mathbb{E}_{\mathcal{E}^{\prime},j}-\mathbb{E}_{\mathcal{E}^{\prime},j-1})[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}]|$
$\displaystyle\leq\delta_{2}\cdot\frac{|(\mathbb{E}_{\mathcal{E}_{1},j}-\mathbb{E}_{\mathcal{E}_{1},j-1})[\mathrm{tr}\,\mathbf{Q}_{-ij}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}\mid\neg\mathcal{E}_{2}]|}{1-\delta_{2}}$
$\displaystyle\leq 2\delta_{2}\cdot 6d\leq 12\cdot
d/m,\quad\text{for}\quad\delta_{2}\leq 1/m.$

So, with $x_{j}=\psi_{j}+\xi_{j}$ and
$X=-\mathrm{tr}(\mathbf{Q}_{-i}-\mathbb{E}_{\mathcal{E}^{\prime}}[\mathbf{Q}_{-i}])\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$
in Lemma 26, for $L=2$ we get:

$\displaystyle\mathbb{E}_{\mathcal{E}^{\prime}}[X^{2}]$ $\displaystyle\leq
C_{2}\cdot\sum_{j}\mathbb{E}_{\mathcal{E}^{\prime}}\big{[}(\psi_{j}+\xi_{j})^{2}\big{]}$
$\displaystyle=C_{2}\cdot\sum_{j}\Big{(}\mathbb{E}_{\mathcal{E}^{\prime}}[\psi_{j}^{2}]+2\,\mathbb{E}_{\mathcal{E}^{\prime}}[\psi_{j}\xi_{j}]+\mathbb{E}_{\mathcal{E}^{\prime}}[\xi_{j}^{2}]\Big{)}$
$\displaystyle\leq C_{2}m\cdot\Big{(}72\cdot 6\cdot\frac{d}{m-d}+2\cdot 6\cdot
12\cdot\frac{d}{m}+12^{2}\frac{d^{2}}{m^{2}}\Big{)}\leq Cd,$

where we also used that $m\geq 8d$, thus concluding the proof.

##  Appendix C Restricted Bai-Silverstein inequality

In this section, we prove Theorem 14. Specifically, we study Condition 2
(Restricted Bai-Silverstein), the second structural condition for small
inversion bias in Theorem 11, which describes the deviation of a quadratic
form $\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}$ from its
mean, for a random vector $\mathbf{x}$. We start by showing Theorem 14, a
generalized version of the lemma of Bai and Silverstein (Lemma 13), which
applies when $\mathbf{x}$ has independent entries. Then, in Appendix C.2 we
show a similar result for a leverage score sparsified vector, constructed as
in Definition 6, which has non-independent entries. Finally, in Appendix C.3
we consider random vectors used in other fast sketching methods, and give
lower bounds demonstrating why these methods do not provide satisfactory
guarantees for Condition 2.

###  C.1 Proof of Theorem 14

Since the assumptions on $\mathbf{x}$ only depend on the leverage scores of
$\mathbf{A}$, and the conclusion is about the variance of a quadratic form,
which only depends on the first four moments of the entries of $\mathbf{x}$,
we can assume without loss of generality that the distribution of $\mathbf{x}$
only depends on the leverage scores of $\mathbf{A}$. We will prove the claim
for such random vectors $\mathbf{x}$.

We start by proving the following result:

######  Proposition 27.

Let $\mathbf{A}$ be a fixed $n\times d$ matrix with $n\geq d$, and
$\mathbf{x}$ be a random vector with independent entries with mean zero and
unit variance, whose distribution only depends on the leverage scores of
$\mathbf{A}$. Then, Condition 2 (Restricted Bai-Silverstein) for the matrix
$\mathbf{A}$ is equivalent to

$\displaystyle\lambda_{\max}\left((\mathbf{U}\circ\mathbf{U})^{\scriptscriptstyle{\top}}\mathbf{D}(\mathbf{U}\circ\mathbf{U})\right)\leq\alpha-2,$

where $\mathbf{U}$ is the $n\times d$ matrix of left singular vectors of
$\mathbf{A}$ and $\mathbf{D}$ is the $n\times n$ matrix
$\mathbf{D}=\operatorname*{\mathop{\mathrm{diag}}}(d_{k})$, with
$d_{k}=\mathbb{E}x_{k}^{4}-3$.

###### Proof.

The Restricted Bai-Silverstein condition is equivalent to having, for all
matrices $\mathbf{B}$ of the form
$\mathbf{B}=\mathbf{U}\mathbf{M}\mathbf{U}^{\scriptscriptstyle{\top}}$, where
$\mathbf{U}$ is the matrix of left singular vectors of $\mathbf{A}$,

$\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}]\leq\alpha\cdot\mathrm{tr}(\mathbf{B}^{2}).$

Let $\mathbf{z}=\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{x}$. Then this is
equivalent to

$\displaystyle\mathrm{Var}[\mathbf{z}^{\scriptscriptstyle{\top}}\mathbf{M}\mathbf{z}]\leq\alpha\cdot\mathrm{tr}(\mathbf{M}^{2}).$
(6)

First we claim that it is enough to consider diagonal matrices $\mathbf{M}$.
Suppose that we have a condition
$C(\operatorname*{\mathop{\mathrm{diag}}}(\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}),\alpha)$
that guarantees that (6) holds for diagonal matrices $\mathbf{M}_{d}$, and
that depends only on the leverage scores and $\alpha$. Now, consider a general
matrix $\mathbf{M}$, and suppose it has the eigendecomposition
$\mathbf{M}=\mathbf{O}\mathbf{M}_{d}\mathbf{O}^{\scriptscriptstyle{\top}}$ for
a diagonal matrix $\mathbf{M}_{d}$. We can write the equivalences

$\displaystyle\mathrm{Var}[\mathbf{z}^{\scriptscriptstyle{\top}}\mathbf{M}\mathbf{z}]$
$\displaystyle\leq\alpha\cdot\mathrm{tr}(\mathbf{M}^{2})$
$\displaystyle\mathrm{Var}[\mathbf{z}^{\scriptscriptstyle{\top}}\mathbf{O}\mathbf{M}_{d}\mathbf{O}^{\scriptscriptstyle{\top}}\mathbf{z}]$
$\displaystyle\leq\alpha\cdot\mathrm{tr}([\mathbf{O}\mathbf{M}_{d}\mathbf{O}^{\scriptscriptstyle{\top}}]^{2})$
$\displaystyle\mathrm{Var}[\mathbf{z}_{d}^{\scriptscriptstyle{\top}}\mathbf{M}_{d}\mathbf{z}_{d}]$
$\displaystyle\leq\alpha\cdot\mathrm{tr}(\mathbf{M}_{d}^{2})$

where
$\mathbf{z}_{d}=\mathbf{O}^{\scriptscriptstyle{\top}}\mathbf{z}=(\mathbf{U}\mathbf{O})^{\scriptscriptstyle{\top}}\mathbf{x}$.
Now we apply the condition
$C(\operatorname*{\mathop{\mathrm{diag}}}(\mathbf{U}_{o}\mathbf{U}_{o}^{\scriptscriptstyle{\top}}),\alpha)$
to $\mathbf{U}_{o}=\mathbf{U}\mathbf{O}$ and the diagonal matrix
$\mathbf{M}_{d}$. This condition is applicable, because $\mathbf{M}_{o}$ is a
diagonal matrix, and guarantees (6) for $\mathbf{M}_{d}$. However, we also
have that the row norms of $\mathbf{U}_{o}=\mathbf{U}\mathbf{O}$ are the same
as the row norms of $\mathbf{U}$, because $\mathbf{O}$ simply acts by an
orthogonal rotation of the rows. So
$\operatorname*{\mathop{\mathrm{diag}}}(\mathbf{U}_{o}\mathbf{U}_{o}^{\scriptscriptstyle{\top}})=\operatorname*{\mathop{\mathrm{diag}}}(\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}})$.
Thus, since the distributions of the sketches we consider only depend on the
leverage scores of $\mathbf{A}$, which are the diagonals of the matrix
$\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}=\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}$,
the condition
$C(\operatorname*{\mathop{\mathrm{diag}}}(\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}),\alpha)$
guarantees that (6) holds for the original matrix $\mathbf{M}$. This shows
that it is enough to establish the condition for diagonal matrices
$\mathbf{M}$.

Hence we can rotate $\mathbf{U}$ by the eigenvectors $\mathbf{O}$ of
$\mathbf{M}$ into $\mathbf{U}^{\prime}=\mathbf{U}\mathbf{O}$, and thus assume
without loss of generality that $\mathbf{M}$ is diagonal,
$\mathbf{M}=\operatorname*{\mathop{\mathrm{diag}}}({\mathbf{g}})$, where
${\mathbf{g}}$ is a vector. Then, the condition simplifies to

$\displaystyle\mathrm{Var}[\mathbf{z}^{\scriptscriptstyle{\top}}\mathbf{M}\mathbf{z}]$
$\displaystyle=\mathrm{Var}[\sum_{i=1}^{d}z_{i}^{2}g_{i}]$
$\displaystyle={\mathbf{g}}^{\scriptscriptstyle{\top}}\mathbf{\Gamma}{\mathbf{g}}\leq\alpha\cdot\|{\mathbf{g}}\|^{2},$

where $\mathbf{\Gamma}$ is the covariance matrix of
$\mathbf{z}\circ\mathbf{z}$. Here the symbol $\circ$ means entrywise product.
This condition has to be true for any vector ${\mathbf{g}}$. Thus, this
condition says exactly that the largest eigenvalue of $\mathbf{\Gamma}$ is at
most $\alpha$:

$\displaystyle\lambda_{\max}(\mathbf{\Gamma})\leq\alpha.$

Also we assume that
$\mathbb{E}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}=\mathbf{I}_{m}$,
hence for any symmetric matrix $\mathbf{F}$ (see e.g., [BS10, CD11] and [MM19,
Lemma B.6.]),

$\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{F}\mathbf{x}]=\sum_{k}d_{k}F_{kk}^{2}+2\mathrm{tr}(\mathbf{F}^{2})$
(7)

where $d_{k}=\mathbb{E}x_{k}^{4}-3$. Therefore, applying this for
$\mathbf{F}=\mathbf{U}\operatorname*{\mathop{\mathrm{diag}}}({\mathbf{g}})\mathbf{U}^{\scriptscriptstyle{\top}}$,
and matching terms, one has
$\mathbf{\Gamma}=(\mathbf{U}\circ\mathbf{U})^{\scriptscriptstyle{\top}}\mathbf{D}\mathbf{U}\circ\mathbf{U}+2\mathbf{I}_{n}$,
where $\mathbf{D}=\operatorname*{\mathop{\mathrm{diag}}}(d_{k})$ and with
$d_{k}=\mathbb{E}x_{k}^{4}-3$. This finishes the proof. ∎

We now continue with the proof of the main claim (Theorem 14). Based on the
above results, as long as the random vector $\mathbf{x}$ has independent
entries of zero mean and unit variance, proving Condition 2 boils down to the
control of the fourth moment of the distribution.

Let $\mathbf{R}=\mathbf{U}\circ\mathbf{U}$, and let $\mathbf{r}_{i}$ denote
its rows. Note that $\mathbf{r}_{i}$ have non-negative entries. Let
$\mathbf{L}=\operatorname*{\mathop{\mathrm{diag}}}(1/\|\mathbf{u}_{i}\|^{2})=\operatorname*{\mathop{\mathrm{diag}}}(1/\|\mathbf{r}_{i}\|_{1})=\operatorname*{\mathop{\mathrm{diag}}}(1/l_{i})$
be the matrix of inverse leverage scores of $\mathbf{A}$, which are also the
inverse $\ell_{1}$ norms of the rows $\mathbf{r}_{i}$ of $\mathbf{R}$. We can
simply discard the zero rows to ensure that this is well defined and
$\|\mathbf{r}_{i}\|_{1}>0$ for all indices.

Then if we can bound
$\lambda_{\max}\left(\mathbf{R}^{\scriptscriptstyle{\top}}\mathbf{L}\mathbf{R}\right)\leq\kappa$,
it follows that
$\lambda_{\max}\left(\mathbf{R}^{\scriptscriptstyle{\top}}\mathbf{D}\mathbf{R}\right)\leq
C\kappa\leq\alpha-2$, which is our desired condition as long as $\alpha$ is
sufficiently large. We will show this bound with $\kappa=1$.

Note that
$\mathbf{Q}=\mathbf{R}^{\scriptscriptstyle{\top}}\mathbf{L}\mathbf{R}$ is a
symmetric matrix and has non-negative entries, because the rows of
$\mathbf{R}$, $\mathbf{r}_{i}=\mathbf{u}_{i}\circ\mathbf{u}_{i}$ are the
entry-wise squares of certain vectors, and the entries of $\mathbf{L}$ are all
positive. Moreover, it is readily verified that the all ones vector
$\mathbf{1}_{d}$ (which clearly has non-negative entries), is an eigenvalue of
$\mathbf{Q}$ with unit eigenvalue,

$\displaystyle\mathbf{Q}\mathbf{1}_{d}=\mathbf{1}_{d}.$

In other words, $\mathbf{Q}$ is a symmetric doubly stochastic matrix. In more
detail, we have

$\displaystyle\mathbf{Q}\mathbf{1}_{d}$
$\displaystyle=\mathbf{R}^{\scriptscriptstyle{\top}}\mathbf{L}\mathbf{R}\mathbf{1}_{d}=\sum_{i=1}^{n}\frac{\mathbf{r}_{i}\mathbf{r}_{i}^{\scriptscriptstyle{\top}}}{\|\mathbf{r}_{i}\|_{1}}\mathbf{1}_{d}=\sum_{i=1}^{n}\mathbf{r}_{i}\cdot\frac{\mathbf{r}_{i}^{\scriptscriptstyle{\top}}\mathbf{1}_{d}}{\|\mathbf{r}_{i}\|_{1}}.$

Now, clearly, since $\mathbf{r}_{i}$ have non-negative entries, we have
$\mathbf{r}_{i}^{\scriptscriptstyle{\top}}\mathbf{1}_{d}=\|\mathbf{r}_{i}\|_{1}$.
Therefore, we find

$\displaystyle\mathbf{Q}\mathbf{1}_{d}$
$\displaystyle=\sum_{i=1}^{n}\mathbf{r}_{i}\cdot\frac{\|\mathbf{r}_{i}\|_{1}}{\|\mathbf{r}_{i}\|_{1}}=\sum_{i=1}^{n}\mathbf{r}_{i}=\mathbf{1}_{d}.$

In the last equality, we have used that, since the columns of $\mathbf{U}$ are
orthogonal vectors, we have that $\sum_{i=1}^{n}r_{ij}=1$ for all
$j=1,\ldots,d$.

Hence, the largest eigenvalue of $\mathbf{Q}$ is at least $1$. By the Perron-
Frobenius theorem for non-negative matrices, it follows that the largest
eigenvalue of $\mathbf{Q}$ is paired with an eigenvector $\mathbf{v}$ of non-
negative entries, see e.g., [Mey00, claims 8.3.1 and 8.3.2]. We can write, for
any such vector $\mathbf{v}\geq 0$, that

$\displaystyle\mathbf{Q}\mathbf{v}$
$\displaystyle=\mathbf{R}^{\scriptscriptstyle{\top}}\mathbf{L}\mathbf{R}\mathbf{v}=\sum_{i=1}^{n}\frac{\mathbf{r}_{i}\mathbf{r}_{i}^{\scriptscriptstyle{\top}}}{\|\mathbf{r}_{i}\|_{1}}\mathbf{v}=\sum_{i=1}^{n}\mathbf{r}_{i}\cdot\frac{\mathbf{r}_{i}^{\scriptscriptstyle{\top}}\mathbf{v}}{\|\mathbf{r}_{i}\|_{1}}.$

Now, clearly
$\mathbf{r}_{i}^{\scriptscriptstyle{\top}}\mathbf{v}/\|\mathbf{r}_{i}\|_{1}\leq\|\mathbf{v}\|_{\infty}$.
Since each entry of each $\mathbf{r}_{i}$ is non-negative, we have that
$0\leq(\mathbf{Q}\mathbf{v})_{j}\leq(\sum_{i=1}^{n}r_{ij})\|\mathbf{v}\|_{\infty}$.
As mentioned, we also have that $\sum_{i=1}^{n}r_{ij}=1$. Hence,

$\displaystyle 0\leq(\mathbf{Q}\mathbf{v})_{j}$
$\displaystyle\leq\|\mathbf{v}\|_{\infty},\,\,j=1,\ldots,n.$

Suppose $\mathbf{v}$ is an eigenvector of $\mathbf{Q}$ with eigenvalue
$\lambda\geq 0$, i.e., $\mathbf{Q}\mathbf{v}=\lambda\mathbf{v}$. Based on the
above inequality, we find
$\|\lambda\mathbf{v}\|_{\infty}\leq\|\mathbf{v}\|_{\infty}$, hence
$\lambda\leq 1$. This shows that the largest eigenvalue of $\mathbf{Q}$ is at
most unity. Thus, by the above reasoning
$\lambda_{\max}\left(\mathbf{R}^{\scriptscriptstyle{\top}}\mathbf{D}\mathbf{R}\right)\leq
C$, and thus Condition 2 holds as long as $C+2\leq\alpha$. This finishes the
proof.

###  C.2 Restricted Bai-Silverstein for LESS embeddings

In this section we show that a sub-gaussian random vector sparsified using our
leverage score sparsifier (LESS) satisfies Condition 2 (Restricted Bai-
Silverstein) with $\alpha=O(1)$. We use this fact in Appendix D to prove
Theorem 8.

######  Lemma 28 (Restricted Bai-Silverstein for LESS).

Fix a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ with rank $d$ and let
${\boldsymbol{\xi}}$ be a leverage score sparsifier for $\mathbf{A}$. For any
p.s.d. matrix $\mathbf{B}$ restricted to the span of $\mathbf{A}$ and any
$\mathbf{x}^{\scriptscriptstyle{\top}}=(x_{1},...,x_{n})$ having independent
entries with mean zero, unit variance and
$\mathbb{E}[\mathbf{x}_{i}^{4}]=O(1)$,

$\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}).$

###### Proof.

Let
$\mathbf{U}=\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1/2}$
be the orthonormal basis matrix for the span of $\mathbf{A}$, and let
$\mathbf{U}_{{\boldsymbol{\xi}}}=\operatorname*{\mathop{\mathrm{diag}}}({\boldsymbol{\xi}})\mathbf{U}$.
Note that
$\mathbf{B}=\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}=\mathbf{U}\mathbf{C}\mathbf{U}^{\scriptscriptstyle{\top}}$
for $\mathbf{C}=\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{U}$. It
follows that:

$\displaystyle\mathrm{Var}\big{[}(\mathbf{x}\circ{\boldsymbol{\xi}})^{\scriptscriptstyle{\top}}\mathbf{B}(\mathbf{x}\circ{\boldsymbol{\xi}})\big{]}$
$\displaystyle=\mathrm{Var}[\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}}\mathbf{x}]=\mathrm{Var}\big{[}\mathrm{tr}(\mathbf{U}_{{\boldsymbol{\xi}}}\mathbf{C}\mathbf{U}_{{\boldsymbol{\xi}}}^{\scriptscriptstyle{\top}})\big{]}+\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{]},$