# Algorithmic Gaussianization through Sketching:
Converting Data into Sub-gaussian Random Designs

Michał Dereziński
Department of Electrical Engineering & Computer Science
University of Michigan
<EMAIL_ADDRESS>

###### Abstract

Algorithmic Gaussianization is a phenomenon that can arise when using
randomized sketching or sampling methods to produce smaller representations of
large datasets: For certain tasks, these sketched representations have been
observed to exhibit many robust performance characteristics that are known to
occur when a data sample comes from a sub-gaussian random design, which is a
powerful statistical model of data distributions. However, this phenomenon has
only been studied for specific tasks and metrics, or by relying on
computationally expensive methods. We address this by providing an algorithmic
framework for gaussianizing data distributions via averaging, proving that it
is possible to efficiently construct data sketches that are nearly
indistinguishable (in terms of total variation distance) from sub-gaussian
random designs. In particular, relying on a recently introduced sketching
technique called Leverage Score Sparsified (LESS) embeddings, we show that one
can construct an $n\times d$ sketch of an $N\times d$ matrix $\mathbf{A}$,
where $n\ll N$, that is nearly indistinguishable from a sub-gaussian design,
in time $O(\mathrm{nnz}(\mathbf{A})\log N+nd^{2})$, where
$\mathrm{nnz}(\mathbf{A})$ is the number of non-zero entries in $\mathbf{A}$.
As a consequence, strong statistical guarantees and precise asymptotics
available for the estimators produced from sub-gaussian designs (e.g., for
least squares and Lasso regression, covariance estimation, low-rank
approximation, etc.) can be straightforwardly adapted to our sketching
framework. We illustrate this with a new approximation guarantee for sketched
least squares, among other examples.

##  1 Introduction

In a standard statistical learning setup, we are given a sample of i.i.d.
points $\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{n}$, and our goal is to
perform an estimation or prediction task. For example, in linear regression,
we aim to learn a linear model $\mathbf{w}^{*}$ from labels/responses
$y_{i}=\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{w}^{*}+\xi_{i}$, where
$\xi_{i}$ represents the noise. Naturally, the performance of prediction
models for linear regression, such as ordinary least squares (OLS) and Lasso,
depends greatly on the properties of the sample distribution, as well as on
the distribution of the noise $\xi$, e.g., whether they are heavy-tailed or
not. To that end, there is extensive literature (see Section 4.2) which
provides precise analysis of statistical and machine learning models under
strong distributional assumptions on the data, one of the most common
assumptions being Gaussianity and sub-gaussianity. When, the data and noise
exhibit small sub-gaussian tails (e.g., by having a constant sub-gaussian
Orlicz norm), then estimators such as OLS and Lasso, as well as model
selection methods such as cross-validation, exhibit provably good, even
optimal, performance. While in modern learning tasks these distributional
assumptions are rarely met, we can often benefit from a great abundance of
cheap data available in many domains. To that end, we ask: Can we leverage
this data abundance to algorithmically introduce sub-gaussianity into a data
distribution, and when is it practical?

A key motivation in this context is randomized sketching, which is useful in
situations when running a learning/estimation algorithm directly on the entire
dataset is computationally prohibitive. Instead, we produce a smaller sketch
of the data, e.g., via importance sampling or random projections, and use it
as a surrogate data sample. To illustrate this, suppose that our goal is to
solve a least squares regression task $(\mathbf{A},\mathbf{b})$:

$\displaystyle\text{find}\qquad\mathbf{w}^{*}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}}L(\mathbf{w})\qquad\text{for}\quad
L(\mathbf{w})=\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2},$

where $\mathbf{A}$ is a large tall $N\times d$ matrix and $\mathbf{b}$ is an
$N$-dimensional vector. In the context of sketching, a classical way of
computationally reducing this task to a much smaller and more tractable
instance is via the so-called Gaussian embedding: We apply a randomized linear
transformation to both $\mathbf{A}$ and $\mathbf{b}$, where the transformation
is defined by an $n\times N$ sketching matrix $\mathbf{S}$ with i.i.d.
Gaussian entries $\mathcal{N}(0,1/n)$ to produce a much smaller $n\times d$
sketched regression task
$(\widetilde{\mathbf{X}},\tilde{\mathbf{y}})=(\mathbf{S}\mathbf{A},\mathbf{S}\mathbf{b})$,
where $n\ll N$. Remarkably, for an arbitrary fixed $(\mathbf{A},\mathbf{b})$,
the resulting algorithmically generated data sample
$(\widetilde{\mathbf{X}},\tilde{\mathbf{y}})$ exactly matches the standard
statistical Gaussian linear model: each row vector
$\tilde{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}$ of
$\widetilde{\mathbf{X}}$ is a Gaussian random vector with covariance
$\mathbf{\Sigma}=\frac{1}{n}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$,
whereas each entry of $\tilde{\mathbf{y}}$ is distributed according to a
linear model
$\tilde{y}_{i}=\tilde{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}\mathbf{w}^{*}+\xi_{i}$
with independent mean zero Gaussian noise.111This follows because
$\tilde{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}\mathbf{w}^{*}$ and
$\tilde{y}_{i}-\tilde{\mathbf{x}}_{i}^{\scriptscriptstyle{\top}}\mathbf{w}^{*}$
are uncorrelated and jointly Gaussian, and therefore, independent. Note that
we did not impose any linear noise model on the original problem
$(\mathbf{A},\mathbf{b})$, but rather it arises naturally through the
sketching transformation. In particular, this allows us to derive the exact
expected approximation error of the sketched ordinary least squares (OLS)
estimator:

$\displaystyle\mathbb{E}\big{[}L(\widehat{\mathbf{w}})-L(\mathbf{w}^{*})\big{]}=\frac{d}{n-d-1}\cdot
L(\mathbf{w}^{*})\qquad\text{for}\qquad\widehat{\mathbf{w}}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}}\|\widetilde{\mathbf{X}}\mathbf{w}-\tilde{\mathbf{y}}\|^{2}.$
(1)

This statistical model implies many other strong performance guarantees for a
variety of tasks such as covariance estimation, OLS, Lasso, PCA etc., which
directly apply to data transformed by Gaussian embeddings, regardless of the
distribution of the original data matrix $\mathbf{A}$. We will refer to this
phenomenon as Algorithmic Gaussianization. Unfortunately, Gaussian embeddings
carry a substantial preprocessing cost, compared to, say, uniformly down-
sampling the data, and so, many other sketching techniques have been proposed,
e.g., Sparse Johnson-Lindenstrauss Transforms (SJLT), Subsampled Randomized
Hadamard Transforms (SRHT) and Leverage Score Sampling [42, 24, 10] (see
Section 4.2), which can be viewed as computationally efficient algorithms for
transforming a large dataset into a small data sample. Existing work has shown
various approximation guarantees for these algorithms, via arguments based on
subspace embeddings and the Johnson-Lindenstrauss property [48, 23, 18].
However, even though some level of gaussianization is implied by these
guarantees, exact parallels with Gaussian random designs are not available.

In this work, we establish a framework for studying Algorithmic
Gaussianization. Specifically, in our main result (Theorem 1), we characterize
the central limit behavior of a sum of i.i.d. random vectors by bounding its
total variation distance from the closest sub-gaussian random design. This
approach can be used to show sub-gaussian properties for extremely sparse
sketching transformations, including the so-called Leverage Score Sparsified
(LESS) embeddings (Corollary 1), whose time complexity matches that of other
state-of-the-art techniques [17]. For these transformations, our result yields
a sub-gaussian property known as the Hanson-Wright inequality [41], which is
widely used to establish both asymptotic and non-asymptotic guarantees for
statistical models. We illustrate the strength of our framework on several
examples, including least squares (Theorem 2), Lasso regression (Corollary 2),
and Randomized SVD (Corollary 3). In particular, as a secondary result of
independent interest, we show that the Gaussian approximation error formula
(1) for the sketched OLS estimator can be non-asymptotically extended to
efficient sketches such as LESS embeddings. Namely, we show that for any
$\mathbf{A}$ and $\mathbf{b}$, if $\mathbf{S}$ is a LESS embedding matrix,
then the sketched OLS estimator satisfies:

$\displaystyle\text{(Theorem \ref{t:main-
app})}\qquad\mathbb{E}\big{[}L(\widehat{\mathbf{w}})-L(\mathbf{w}^{*})\big{]}\approx_{1+\epsilon}\frac{d}{n-d-1}\cdot
L(\mathbf{w}^{*})\qquad\text{for}\qquad\epsilon=\tilde{O}(1/\sqrt{d}).$

Here, $a\approx_{1+\epsilon}b$ denotes a $(1+\epsilon)$-approximation
$b/(1+\epsilon)\leq a\leq(1+\epsilon)b$. This type of guarantee for sketched
OLS is a first of its kind (outside of Gaussian embeddings), as prior work
only showed upper bounds (with additional constant and logarithmic factors),
e.g., of the form $\tilde{O}(d/n)\cdot L(\mathbf{w}^{*})$ [19], or asymptotic
results [20]. Remarkably, we show empirically that for LESS embeddings the
above approximation error estimate is extremely accurate, whereas for other
popular sketching methods that are not covered by our theory (such as row
subsampling and SRHT), the approximation error can be much more problem
dependent.

##  2 Main result

In this section, we present our main result, which bounds the number of
samples from a multivariate distribution needed to produce a gaussianized
sample, i.e., one that is close to a sub-gaussian random design in terms of
the total variation distance. Our result relies on the same gaussianization
procedure that arises in the classical central limit theorem, but we
deliberately relax the target distribution to sub-gaussian so that we can
obtain a stronger convergence guarantee.

As there are a number of closely related notions of sub-gaussianity, we
clarify the notation here, with a more detailed discussion in Section 4. We
say that a variable $X$ is $K$-sub-gaussian if its corresponding sub-gaussian
Orlicz norm is bounded by $K$, defined as
$\|X\|_{\psi_{2}}=\inf\\{t>0:\,\mathbb{E}\,\exp(X^{2}/t^{2})\leq 2\\}$. For a
random vector, perhaps the simplest and most popular model of sub-gaussianity
is to assume that it has independent sub-gaussian entries. However, entry-wise
independence is a fairly strong assumption. The multivariate sub-gaussian
norm, i.e.,
$\|\mathbf{x}\|_{\psi_{2}}=\sup_{\mathbf{v}:\|\mathbf{v}\|=1}\|\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{x}\|_{\psi_{2}}$,
is a popular relaxation that allows for sub-gaussian vectors with dependent
entries (we call them $K$-sub-gaussian vectors). However, this notion is too
weak for some applications, e.g., in high-dimensional statistics, or to obtain
our least squares result (Theorem 2). These settings often rely on the
following strictly stronger property, satisfied by random vectors with
independent sub-gaussian entries, called the Hanson-Wright inequality. Below,
$c$ denotes an absolute constant.

######  Lemma 1 (Hanson-Wright inequality, [41])

Let $\mathbf{x}$ have independent $K$-sub-gaussian entries with mean zero and
unit variance. Then, it satisfies the Hanson-Wright inequality with constant
$K$:

$\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}}{K^{4}\|\mathbf{B}\|_{F}^{2}},\frac{t}{K^{2}\|\mathbf{B}\|}\Big{\\}}\bigg{)}\qquad\text{for
any \ $\mathbf{B}$}.$ (2)

###  2.1 Sub-gaussian limit theorem via new Hanson-Wright inequality

There is extensive literature dedicated to extending and adapting the Hanson-
Wright inequality (e.g., [29, 1, 46, 5]). Yet among popular sketching methods,
only Gaussian and sub-gaussian embeddings are known to exhibit this property
with a small constant $K=O(1)$. To address this, we introduce the following
technical result, which shows that, by relying on the central limit
phenomenon, we can establish that a sum of bounded random vectors gets very
close to a random vector that satisfies Hanson-Wright with a small constant,
even though the sum may not satisfy this property directly (proof in Appendix
C). Below, we use $d_{{\mathrm{tv}}}$ to denote the total variation distance
(see Appendix B for notation).

###### Theorem 1

Let $\mathbf{x}_{1},\mathbf{x}_{2},...$ be i.i.d. vectors, where
$\mathbb{E}[\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}]=\mathbf{\Sigma}$
is full-rank and $\|\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}\|_{2}$ is $M$-sub-
gaussian. Define a $k$-gaussianized sample
$\tilde{\mathbf{x}}=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}r_{i}\mathbf{x}_{i}$,
where $r_{i}$ are independent mean zero unit variance $R$-sub-gaussian random
variables. For any $\delta>0$, there is
$L=O\big{(}1+M\log(M/\delta)/\sqrt{k}\big{)}$ and a mean zero isotropic vector
$\mathbf{z}$ satisfying the Hanson-Wright inequality (2) with constant $K=LR$,
such that:

$\displaystyle
d_{{\mathrm{tv}}}\big{(}\,\tilde{\mathbf{x}},\,\widetilde{\mathbf{\Sigma}}^{1/2}\mathbf{z}\,\big{)}\leq\delta\qquad\text{for}\qquad\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{\Sigma}.$

###### Remark 1

If we make a stronger assumption that
$\|\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}\|_{2}\leq M$ almost surely, then we
can obtain a slightly sharper result: $L=O(1+M\sqrt{\log(M/\delta)/k})$. Note
that
$M^{2}\geq\mathbb{E}[\|\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}\|_{2}^{2}]=d$,
where $d$ is the dimension of $\mathbf{x}_{i}$. Either of the boundedness
assumptions includes all distributions over finite populations in general
position.

###### Remark 2

Crucially, random vector $\mathbf{z}$ may have dependent entries (in which
case, it does not satisfy the assumptions of Lemma 1), yet we show that it
behaves much like a vector with independent sub-gaussian entries.

At a high level, Theorem 1 states that as $k$ goes to infinity, the total
variation distance of $\tilde{\mathbf{x}}$ from a certain family of sub-
gaussian distributions goes to zero. Using the total variation distance to
measure central limit behavior (instead of, say, the Wasserstein distance) is
unusual because of how strong of a guarantee it requires. This is possible
here because we are comparing to a sub-gaussian distribution, instead of
directly to a Gaussian. To highlight the benefits of our approach, recall that
small total variation distance means that we can couple the $k$-gaussianized
sample $\tilde{\mathbf{x}}$ with a scaled sub-gaussian sample
$\widetilde{\mathbf{\Sigma}}^{1/2}\mathbf{z}$ so that they are identical with
probability $1-\delta$. Using a union bound, we can do the same for a
collection of $k$-gaussianized samples
$\tilde{\mathbf{x}}_{1},...,\tilde{\mathbf{x}}_{n}$. This effectively allows
us to transfer any property that holds with high probability for the sub-
gaussian random design $\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}$ to the
$k$-gaussianized design $\widetilde{\mathbf{X}}$. The only caveat in this
reduction is that the covariance matrix gets slightly distorted in the
process, but both this distortion and the failure probability are controlled
by $\delta$ which can be made negligibly small. We illustrate this argument in
Section 2.2 when we apply the result to LESS embeddings. In Section 5, we
establish that our result is nearly tight: There are distributions for which
we need to gaussianize $k=\Omega(M^{2})$ samples to obtain
$\tilde{\mathbf{x}}$ that is $\delta$-close to an $O(1)$-sub-gaussian sample
in the sense of Theorem 1, even for $\delta=\Omega(1)$. This, up to
logarithmic factors, matches our upper bound of $\tilde{O}(M^{2})$ samples.

###  2.2 Implementation with LESS embeddings

Let us connect our notion of a $k$-gaussianized sample to the randomized
sketching framework which motivated our approach. Theorem 1 can be applied to
sketching matrices known as Leverage Score Sparsified (LESS) embeddings. These
embeddings are essentially just sparse random matrices (e.g., with random sign
entries) where the sparsity pattern is chosen via an importance sampling
method called leverage score sampling [24]. LESS embeddings have been shown,
for particular tasks (namely, inverse covariance estimation [17] and
stochastic optimization [14]), to produce estimators that are nearly
indistinguishable from their Gaussian counterparts, at a fraction of the cost.
In the below definition, which is a slight variation222We generate a separate
random variable $r_{i}$ for each sampled index, whereas [17] propose to
generate $r_{i}$’s only for each _unique_ sampled index (this is different
since we are sampling with replacement). Our definition of a LESS embedding
fits more naturally with Theorem 1, but all our results could be easily
converted to the version from [17]. Also, unlike [17], we restrict random
variables $r_{i}$ to $\pm 1$ Rademacher variables (instead of arbitrary sub-
gaussian), and require an $\Omega(1/N)$ lower bound on the probabilities
$p_{i}$ to avoid worst-case numerical precision issues in Theorem 2. on the
one proposed by [17], $\mathbf{e}_{i}$ denotes the $i$th standard basis vector
in $\mathbb{R}^{N}$.

######  Definition 1 (LESS: LEverage Score Sparsified embedding [17])

Consider a full rank $N\times d$ matrix $\mathbf{A}$ and let
$l_{i}=\mathbf{a}_{i}^{\scriptscriptstyle{\top}}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{a}_{i}$
for $i=1,...,N$ denote its leverage scores. Let $I_{1},...,I_{k}$ be sampled
i.i.d. from a probability distribution $(p_{1},...,p_{N})$ such that
$p_{i}\geq(l_{i}/d+1/N)/C$ for all $i$ and an absolute constant $C$. We call
sketching matrix $\mathbf{S}$ a LESS embedding of size $n$ for $\mathbf{A}$
with $k$ non-zeros per row, if it consists of $n$ i.i.d. row vectors
distributed as
$\frac{1}{\sqrt{n}}\sum_{i=1}^{k}\frac{r_{i}}{\sqrt{kp_{I_{i}}}}\mathbf{e}_{I_{i}}^{\scriptscriptstyle{\top}}$,
where $r_{i}$ are independent random $\pm 1$ variables.

Unless specified otherwise, we will let the number of non-zeros per row be
$k=d$, as recommended by [17]. Using a standard technique from randomized
linear algebra, called approximate leverage score sampling [24, 22], we can
implement a LESS embedding in time $O({\mathrm{nnz}}(\mathbf{A})\log
N+nd^{2})$, where ${\mathrm{nnz}}(\mathbf{A})$ is the number of non-zero
entries in $\mathbf{A}$. Thus, this sketching method is nearly linear in the
input size, scales well with the sparsity of the data, and its complexity
matches the state-of-the-art sketching complexity up to logarithmic factors
for most problems, such as least squares and Lasso regression. Our Theorem 1
sheds new light on the Gaussian-like performance of LESS embeddings via a
reduction to a sub-gaussian design.

###### Corollary 1

Given any full rank $N\times d$ matrix $\mathbf{A}$ and any corresponding
$n\times N$ LESS embedding random matrix $\mathbf{S}$, for any $\delta>0$,
there is an $n\times d$ matrix $\mathbf{Z}$ consisting of i.i.d. mean zero
isotropic row vectors $\mathbf{z}_{i}^{\scriptscriptstyle{\top}}$, each
satisfying the Hanson-Wright inequality (2) with constant
$K=O(\sqrt{\log(nd/\delta)})$, such that:

$\displaystyle
d_{{\mathrm{tv}}}\big{(}\,\sqrt{n}\,\mathbf{S}\mathbf{A},\,\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}\,\big{)}\leq\delta\qquad\text{for}\qquad\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}.$

###### Remark 3

The claim is stated for standard LESS embeddings with $d$ non-zeros per row.
If we increase the non-zeros per row to $d\log(nd/\delta)$, then we can
further reduce the Hanson-Wright constant to $K=O(1)$.

We now briefly explain how the corollary arises from Theorem 1. Suppose that
random vectors $\mathbf{x}_{i}$ are defined as uniform row samples from the
data matrix $\mathbf{A}$, scaled by $N$ so that
$\mathbf{\Sigma}=\mathbb{E}[\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}]=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$.
Now, let us construct an $n\times d$ random matrix $\widetilde{\mathbf{X}}$
consisting of rows $\tilde{\mathbf{x}}_{1},...,\tilde{\mathbf{x}}_{n}$ which
are independent $k$-gaussianized samples defined as in Theorem 1. This sketch
can be obtained by applying a sparse $n\times N$ sketching matrix
$\sqrt{n}\mathbf{S}$ to $\mathbf{A}$, where each row of $\sqrt{n}\mathbf{S}$
is given by
$\frac{1}{\sqrt{k}}\sum_{i=1}^{k}r_{i}\mathbf{e}_{I_{i}}^{\scriptscriptstyle{\top}}$
and $I_{i}$ are uniformly random indices (in Section 6, we empirically study a
variant of this simple sketch, called LessUniform). Notice that the constant
$M$ in Theorem 1 is bounded by $N\cdot\max_{i}\ell_{i}$, where $\ell_{i}$ is
the statistical leverage score of the $i$th row $\mathbf{a}_{i}$ of matrix
$\mathbf{A}$, i.e.,
$\ell_{i}=\mathbf{a}_{i}^{\scriptscriptstyle{\top}}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{a}_{i}=\|\mathbf{\Sigma}^{-1/2}\mathbf{a}_{i}\|^{2}$,
with $\mathbf{\Sigma}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$. To
avoid this dependence on the maximum leverage score (i.e., the _coherence_ of
matrix $\mathbf{A}$), we can transform the data distribution so that, instead
of sampling each row uniformly at random, it samples each row with the
probability $\propto\ell_{i}/d$, and then scales the row by
$\sqrt{d/\ell_{i}}$ to maintain the same covariance $\mathbf{\Sigma}$ (here,
$d$ is the normalization constant for the distribution, because
$\sum_{i=1}^{N}\ell_{i}=d$). As a result, we obtain LESS embeddings and ensure
that $M=O(\sqrt{d})$. Note that we can rely here on the stronger boundedness
condition as in Remark 1, so that by letting $k=d$ in Theorem 1, we can get
$L=O(\sqrt{\log(d/\delta)})$. It remains to replace $\delta$ with $\delta/n$
and apply a union bound over the $n$ rows, to show the total variation
distance of $\delta$ between the entire sketch $\sqrt{n}\mathbf{S}\mathbf{A}$
and the corresponding sub-gaussian design.

##  3 Applications

We next present some of the implications of Theorem 1. In Section 3.1, we
present our main application, giving a precise characterization of the
expected approximation error of the sketched least squares estimator with LESS
embeddings. Here, this is not merely a corollary of our main result, but
rather a result of independent interest, and to our knowledge, it was not
previously known even for sub-gaussian embeddings. Then, in Section 3.2, we
demonstrate how our results imply new efficient algorithms for approximately
solving Lasso regression and other constrained optimization tasks. Further
applications to low-rank approximation via Randomized SVD, as well as
corollaries for the subspace embedding and Johnson-Lindenstrauss properties of
LESS embeddings, are given in Appendix A.

###  3.1 Sketched least squares with LESS embeddings

As a motivating application of our theory, we study the approximation
properties of the sketched least squares estimator, showing that, when the
sketch is a LESS embedding, then it achieves nearly the same expected
approximation error as for a Gaussian embedding, even down to constant
factors.

Consider an $N\times d$ data matrix $\mathbf{A}$ with full column rank, and an
$N$-dimensional vector $\mathbf{b}$. Recall that our goal is to approximately
solve the least squares task $(\mathbf{A},\mathbf{b})$ with the loss function
$L(\mathbf{w})=\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2}$, i.e., find:

$\displaystyle\mathbf{w}^{*}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}}L(\mathbf{w})=(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{b}.$

Given an $n\times N$ sketching matrix $\mathbf{S}$, the sketched least squares
estimator $\widehat{\mathbf{w}}$ for this model is the solution of the
sketched version of the problem, i.e.,
$(\mathbf{S}\mathbf{A},\mathbf{S}\mathbf{b})$. If $\mathbf{S}$ is a Gaussian
embedding, then via a reduction to a Gaussian random design, the expected
error of this estimator is given by the exact formula (1), but the cost of
constructing such a sketch is prohibitively expensive. If we let $\mathbf{S}$
be a LESS embedding matrix with $d$ non-zeros per row, then the cost of
computing this estimator is $O({\mathrm{nnz}}(\mathbf{A})\log N+nd^{2})$,
which includes the cost of approximating the leverage scores, sketching the
data, and computing $\widehat{\mathbf{w}}$ from the sketch. We show that the
resulting estimator enjoys nearly the same expected error guarantee as the
Gaussian embedding. To avoid numerical precision issues which may occur when
computing the exact expectation (there is a small but non-zero probability
that $\mathbf{S}\mathbf{A}$ is very ill-conditioned), we will restrict the
entries of the sketched solution $\widehat{\mathbf{w}}$ to a finite range
$\mathcal{R}=[-D,D]$, for a sufficiently large
$D={\mathrm{poly}}(N,\kappa(\mathbf{A}),\|\mathbf{A}\|,\|\mathbf{b}\|)$ which
is absorbed by the logarithmic factors in the result (see Appendix D for
details).

###### Theorem 2

Fix $\mathbf{A}$ with rank $d$ and let $\mathbf{S}$ be a LESS embedding with
$n\geq\tilde{O}(d)$ rows. Then, for any $\mathbf{b}\in\mathbb{R}^{N}$, the
estimator
$\widehat{\mathbf{w}}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}\in\mathcal{R}^{d}}\|\mathbf{S}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2}$
satisfies:

$\displaystyle\mathbb{E}\big{[}L(\widehat{\mathbf{w}})-L(\mathbf{w}^{*})\big{]}\,\approx_{1+\epsilon}\frac{d}{n-d-1}\cdot
L(\mathbf{w}^{*})\qquad\text{for}\qquad\epsilon=\tilde{O}(1/\sqrt{d}\,),$

where $L(\mathbf{w})=\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2}$ and
$\mathbf{w}^{*}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}}L(\mathbf{w})$.

Observe that the obtained expression for the expected multiplicative
approximation error, i.e., $\frac{d}{n-d-1}$, is completely independent of the
data matrix $\mathbf{A}$ or the regression vector $\mathbf{b}$, and it matches
the formula (1) for the Gaussian embedding. Comparable non-asymptotic
guarantees for other fast sketching methods, such as SJLT, SRHT, or leverage
score sampling, take the form of $\tilde{O}(d/n)\cdot L(\mathbf{w}^{*})$,
where the precise constants are problem-dependent.

The result essentially holds for any sketching matrix whose rows satisfy the
Hanson-Wright claim of Theorem 1 with a sufficiently small constant, and to
our knowledge, it was not known even for sub-gaussian embeddings. Moreover,
our proof of Theorem 2 is itself of independent interest, since it proceeds
differently from the existing bounds for sketched least squares (see Appendix
D). Our approach is related to the asymptotic analysis of the OLS estimator in
high-dimensional statistics, see e.g., [27], but we must take additional care
to obtain a non-asymptotic bound. We start by rewriting the expected loss
$\mathbb{E}[L(\widehat{\mathbf{w}})]$ as the expectation of the leave-one-out
cross-validation estimator, which is easier to work with than the loss itself.
We then show that all of the leverage scores of a LESS embedding sketch are
almost exactly uniform (in a much stronger sense than for a leverage score
sampling sketch). Finally, we extend a simple formula for the expectation of a
sketch of the so-called _hat matrix_ (which is the orthogonal projection onto
the column-span of $\mathbf{A}$), from Gaussian embeddings to LESS and sub-
gaussian embedding matrices (see Lemma 8).

In Section 6 we empirically show that our estimate for the expected
approximation error of sketched OLS with a LESS embedding is even more
accurate than suggested by the theory, on a range of benchmark regression
tasks. Our experiments also show that this problem-independent expected
approximation error is not shared by other fast sketching methods that are
considered to have strong guarantees for least squares, such as SRHTs.

###  3.2 Lasso and constrained least squares

In this section, we show how our results can be used to extend existing
sketching guarantees in constrained least squares optimization [40] from sub-
gaussian embeddings to LESS embeddings. This general problem setting includes
such standard tasks as Lasso regression, sparse recovery and training a
support vector machine.

We consider the following constrained optimization task for
$\mathbf{A}\in\mathbb{R}^{N\times d}$ and $\mathbf{b}\in\mathbb{R}^{N}$:

$\displaystyle\text{find}\qquad\mathbf{w}^{*}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}\in\mathcal{C}}L(\mathbf{w}),\qquad\text{for}\qquad
L(\mathbf{w})=\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2},$

where $\mathcal{C}$ is some convex subset of $\mathbb{R}^{d}$. Our goal is to
approximately solve this task by solving a sketched version of the problem:

$\displaystyle\widehat{\mathbf{w}}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}\in\mathcal{C}}\|\mathbf{S}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2},$

for an $n\times N$ sketching matrix $\mathbf{S}$. When the sketching matrix is
a sub-gaussian embedding, then [40] showed that the sketch size $n$ needed for
$\widehat{\mathbf{w}}$ to achieve a $(1+\epsilon)$-approximation, i.e.,
$L(\widehat{\mathbf{w}})\leq(1+\epsilon)\cdot L(\mathbf{w}^{*})$, is
controlled by the so-called Gaussian width, which can be viewed as a measure
of the degrees of freedom for the constrained optimization task (in
particular, it can be much smaller than the actual dimension $d$), and is
defined as follows:

$\displaystyle\mathbb{W}:=\mathbb{W}(\mathbf{A}\mathcal{K})=\mathbb{E}_{{\mathbf{g}}}\Big{[}\sup_{\mathbf{u}\in\mathbf{A}\mathcal{K}\cap\mathcal{S}^{N-1}}|{\mathbf{g}}^{\scriptscriptstyle{\top}}\mathbf{u}|\Big{]},$

where $\mathcal{K}$ is the tangent cone of the constraint set $\mathcal{C}$ at
the optimum $\mathbf{w}^{*}$, i.e., the closed convex hull of
$\\{\Delta\in\mathbb{R}^{d}\mid\Delta(\mathbf{w}-\mathbf{w}^{*})\text{ for
}t\geq 0,\mathbf{w}\in\mathcal{C}\\}$. Specifically, [40] showed that if
$\mathbf{S}$ consists of i.i.d. mean zero isotropic sub-gaussian rows and
$n\geq O(\mathbb{W}^{2}/\epsilon^{2})$, then with high probability the
sketched estimate $\widehat{\mathbf{w}}$ is a $(1+\epsilon)$-approximation of
$\mathbf{w}^{*}$. We can use our framework to extend this result to LESS
embeddings (see Appendix E).

###### Theorem 3

If $\mathbf{S}$ is a LESS embedding matrix for $\mathbf{A}$, with size
$n\geq\tilde{O}(\mathbb{W}^{2}/\epsilon^{2})$, then with high probability
$\widehat{\mathbf{w}}$ is a $(1+\epsilon)$-approximation of $\mathbf{w}^{*}$.

Following [40], we illustrate this claim with a corollary for Lasso
regression, where the constraint set is an $\ell_{1}$ ball, i.e.,
$\mathcal{C}=\\{\mathbf{w}:\|\mathbf{w}\|_{1}\leq R\\}$. Crucially, here the
sketch size needed to obtain a good approximation scales with the sparsity of
the solution $\mathbf{w}^{*}$, rather than with the dimension of the data.
First, we formulate an $\ell_{1}$-restricted condition number of $\mathbf{A}$,
which is always smaller than the actual condition number, and which naturally
arises in bounding the Gaussian width:

$\displaystyle\kappa=\max_{i}\|\mathbf{A}_{:,i}\|/\gamma_{s}^{-}(\mathbf{A}),\qquad\text{where}\qquad\gamma_{s}^{-}(\mathbf{A})=\min_{\|\mathbf{v}\|=1,\,\|\mathbf{v}\|_{1}\leq
2\sqrt{s}}\|\mathbf{A}\mathbf{v}\|.$

###### Corollary 2

Suppose that the Lasso solution $\mathbf{w}^{*}$ is $s$-sparse. Then, using a
LESS embedding of size $n\geq\tilde{O}(\kappa^{2}s/\epsilon^{2})$, we can
obtain a $(1+\epsilon)$-approximation of $\mathbf{w}^{*}$ in time
$O({\mathrm{nnz}}(\mathbf{A})\log N+nd^{2})$.

Finally, we note that similar results for approximately solving Lasso and
other constrained least squares problems can be obtained for SRHT embeddings,
also given by [40]. However, there are two key differences: (1) the bounds
obtained for SRHT are weaker than the ones achieved for LESS or sub-gaussian
embeddings (e.g., in addition to the Gaussian width, they also scale with
other problem-dependent quantities such as the Rademacher width, see [40,
Theorem 2]); and (2) the cost of implementing the SRHT is $O(Nd\log N)$, i.e.,
independent of the sparsity of the data matrix $\mathbf{A}$. This is higher
than the cost of LESS when ${\mathrm{nnz}}(\mathbf{A})\ll Nd$.

##  4 Background and related work

In this section, we put our results in context by discussing different types
of sub-gaussianity that have been considered in the literature. Here, we rely
on the hierarchy of sub-gaussian concentration, shown in Figure 1, which is
based on a range of classical results in high-dimensional probability (see
[45] for an overview). We then discuss how these concepts arise in related
work, particularly in the context of randomized sketching.

###  4.1 Sub-gaussian concentration hierarchy for random vectors

While the notion of sub-gaussian concentration for scalar random variables can
be naturally represented by the Orlicz norm $\|\cdot\|_{\psi_{2}}$, the
landscape of sub-gaussian concentration becomes more complex when we consider
multivariate distributions. As discussed in Section 2, the natural
multivariate extension of the Orlicz norm, i.e., the sub-gaussian norm
$\|\mathbf{x}\|_{\psi_{2}}=\sup_{\mathbf{v}:\|\mathbf{v}\|=1}\|\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{x}\|_{\psi_{2}}$,
while useful for some applications (such as Theorem 3), does not suffice for
other Gaussian-like guarantees (such as Theorem 2). A more general approach to
quantifying sub-gaussianity of random vectors is to analyze the concentration
of real-valued functions of these vectors. In this context, let $\mathcal{F}$
be some subset of functions $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ with
bounded Lipschitz constant:
$\|f\|_{{\mathrm{Lip}}}=\sup_{\mathbf{u},\mathbf{v}}|f(\mathbf{u})-f(\mathbf{v})|/\|\mathbf{u}-\mathbf{v}\|$.

###### Definition 2

We say that a $d$-dimensional random vector $\mathbf{x}$ has the concentration
property over $\mathcal{F}$ with constant $K$ if
$\|f(\mathbf{x})-\mathbb{E}\,f(\mathbf{x})\|_{\psi_{2}}\leq
K\cdot\|f\|_{{\mathrm{Lip}}}$ for all $f\in\mathcal{F}$.

For the sake of simplicity, we focus on random vectors that are mean zero,
i.e., $\mathbb{E}[\mathbf{x}]=\mathbf{0}$, and isotropic, i.e.,
$\mathbb{E}[\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}]=\mathbf{I}$.
Gaussian vectors (as well as, for example, random vectors on the sphere)
satisfy the concentration property over all Lipschitz functions, with an
absolute constant $K$. However, as this condition is very restrictive, smaller
function families have been considered. If we restrict $\mathcal{F}$ to convex
functions, then vectors with independent bounded entries satisfy the
concentration property (but vectors with i.i.d. sub-gaussian entries may not).
On the other end of the concentration hierarchy depicted in Figure 1 is the
family of linear functions, i.e.,
$\mathcal{F}=\\{f:f(\mathbf{x})=\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{x},\text{
for some }\mathbf{v}\in\mathbb{R}^{d}\\}$, for which the concentration
constant of $\mathbf{x}$ becomes simply its sub-gaussian norm. Of primary
interest to this work is the family of Euclidean functions [46], defined as
$\mathcal{F}=\\{f:f(\mathbf{x})=\sqrt{\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{x}}\
\text{ for psd }\,\mathbf{B}\\}$, which falls between convex and linear in the
hierarchy, and can be used to characterize vectors with i.i.d. sub-gaussian
entries. Any mean zero isotropic vector that satisfies the Hanson-Wright
inequality (2) also satisfies Euclidean concentration.

###### Proposition 1

[45, Section 6.3] If an isotropic mean zero vector $\mathbf{x}$ satisfies the
Hanson-Wright inequality (2) with constant $K$, then it is $O(K)$-sub-gaussian
and satisfies Euclidean concentration with constant $O(K^{2})$.

Examples
$\mathbf{x}\in\mathbb{R}^{d}$ Gaussian vectors Uniform on the sphere i.i.d.
bounded entries i.i.d. sub-gaussian entries Hanson-Wright vectors sub-gaussian
vectors Concentration
$\mathcal{F}\subseteq\\{\mathbb{R}^{d}\\!\rightarrow\\!\mathbb{R}\\}$
Lipschitz functions Convex functions Euclidean functions
$f(\mathbf{x})=\sqrt{\mathbf{x}^{\scriptscriptstyle{\top}}\\!\mathbf{B}\mathbf{x}}$
Linear functions
$f(\mathbf{x})=\mathbf{v}^{\scriptscriptstyle{\top}}\\!\mathbf{x}$ Figure 1:
Hierarchy of sub-gaussian concentration for mean zero isotropic random
vectors, with arrows indicating implications. Right column represents
different function classes $\mathcal{F}$ used in Definition 2, and left column
provides examples. “Hanson-Wright vectors” refers to vectors satisfying the
Hanson-Wright inequality (2), as in Theorem 1, whereas “sub-gaussian vectors”
refers to vectors with bounded multivariate sub-gaussian norm
$\|\cdot\|_{\psi_{2}}$.

###  4.2 Related work

Randomized sketching has emerged out of theoretical computer science, as part
of the broader area of Randomized Numerical Linear Algebra (RandNLA; [48, 23,
18]), including computationally efficient techniques such as the Subsampled
Randomized Hadamard Transform (SRHT; [2]), Leverage Score Sampling [24], the
CountSketch [10], Sparse Johnson-Lindenstrauss Transforms (SJLT; [39]),
Leverage Score Sparsified embeddings (LESS; [17]) and Determinantal Point
Processes (DPPs; [11]). More recently, there has been an increased interest in
the statistical analysis of sketching and importance sampling techniques, for
example in the context of linear regression [35], kernel ridge regression [3]
experimental design [12, 16], model averaging [47], and the bootstrap [33].
Some works have shown deep connections between the performance of certain
sketching methods and Gaussian/sub-gaussian embeddings: for DPPs, e.g., in the
context of low-rank approximation and stochastic optimization [15, 13, 38];
for LESS embeddings in the context of inverse covariance estimation [17] and
stochastic optimization [14] (see Section 2.2); and for SRHT, a similar
phenomenon has been studied asymptotically in linear regression [20] and
stochastic optimization [31]. Despite this extensive literature, no general
equivalence result was known between any of these fast sketching methods and
sub-gaussian embeddings. Our results show that such equivalence can be
established for LESS embeddings.

Statistical analysis of random design models has a long history, with
important connections to both asymptotic and non-asymptotic random matrix
theory (e.g., see [4, 45]), which we could hardly summarize here. Naturally, a
variety of random designs have been considered, of which the most relevant to
our work are Gaussian and sub-gaussian designs. These models have proven
extremely useful in understanding the performance of a variety of linear
regression estimators [21, 27, 6, 37], sample covariance estimators [30, 32]
and others. Here, a sub-gaussian random design typically refers to a matrix
$\mathbf{Z}\mathbf{\Sigma}^{1/2}$, where $\mathbf{Z}$ consists of independent
rows which have either bounded sub-gaussian norm or i.i.d. sub-gaussian
entries, or more generally, satisfy a concentration property for some family
of Lipschitz functions (e.g., see [34]). These distinctions are important
here, because our main result relates a $k$-gaussianized sample with a sub-
gaussian vector that may have dependent entries, and yet, similarly to vectors
with independent entries, it satisfies the Hanson-Wright inequality.

Our results are related to the study of the rates of convergence in the
multivariate central limit theorem (CLT), with the key difference that we
relax our notion of gaussianity, which allows us to use a stronger notion of
distance (i.e., the total variation distance, together with the Hanson-Wright
constant). Nevertheless, it is helpful to compare our approach with the
convergence rates for the multivariate CLT in terms of the Wasserstein
distance [9, 7]. Using the setup from Theorem 1, under the assumption that
$\|\mathbf{\Sigma}^{-1/2}\mathbf{x}\|\leq M$ almost surely, the currently best
known Wasserstein CLT of order 2 yields
$W_{2}(\tilde{\mathbf{x}},\mathcal{N}(\mathbf{0},\mathbf{\Sigma}))=O(M\sqrt{d\log(k)/k})$
[25], whereas our bound on the Hanson-Wright constant is
$O(1+M\sqrt{\log(M)/k})$. The guarantees are not directly comparable, but one
could argue that both bounds are useful primarily when $k$ is large enough to
absorb the dependence on $M$ and $d$, i.e., $k=O(M^{2}d\log(M))$ for the
Wasserstein bound and $k=O(M^{2}\log(M))$ for Theorem 1. While the Wasserstein
distance can be used to bound the difference between expectations of Lipschitz
functions, it is non-trivial to effectively bound the Lipschitz constant for
most of the quantities considered in this work, such as the least squares
approximation error, and using this approach would require a more complex
case-by-case analysis.

##  5 Lower bound

In this section we demonstrate the sharpness of our main result, Theorem 1.
Specifically, we lower bound the sub-gaussian norm of a gaussianized vector
(which is even stronger than lower bounding the Hanson-Wright constant), and
all vectors close to it in total variation distance. Consider the following
isotropic sparse random vector $\mathbf{x}$ with a parameter $M$ which
controls the sparsity (and sub-gaussianity) of $\mathbf{x}$:

$\displaystyle\mathbf{x}=M\cdot[t_{1}b_{1},...,t_{d}b_{d}]^{\scriptscriptstyle{\top}}\quad$
$\displaystyle\text{where}\quad t_{i}\sim\mathrm{Rademacher},\quad
b_{i}\sim\mathrm{Bernoulli}(1/M^{2}).$ (3)

Note that $\mathbf{x}$ is isotropic, i.e.,
$\mathbf{\Sigma}=\mathbb{E}[\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}]=\mathbf{I}$.
Next, let $\mathbf{x}_{1},\mathbf{x}_{2},...$ be i.i.d. samples distributed
identically to $\mathbf{x}$. We will study the $k$-gaussianized sample from
this distribution, i.e.,
$\tilde{\mathbf{x}}=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}r_{i}\mathbf{x}_{i}$,
where $r_{i}\sim\mathrm{Rademacher}$ are independent random sign variables.
First, observe that there is an absolute constant $c>1$ such that, for any
$k$:

$\displaystyle\|\tilde{\mathbf{x}}\|_{\psi_{2}}\approx_{O(1)}\|\mathbf{x}\|_{\psi_{2}}\approx_{O(1)}M/\sqrt{\ln
M}.$ (4)

(Note that replacing $\mathbf{x}$ by $\mathbf{\Sigma}^{-1/2}\mathbf{x}$ does
not affect the claim since $\mathbf{\Sigma}=\mathbf{I}$.) In other words,
gaussianization does not directly reduce the sub-gaussian norm (or the Hanson-
Wright constant) of the sample, which is why in Theorem 1 we must construct a
separate random variable $\mathbf{z}$ which is close to $\mathbf{x}$ in total
variation distance. The second approximation in (4) follows by using the fact
that for $b\sim\mathrm{Bernoulli}(\delta)$ we have:
$\|b\|_{\psi_{2}}\approx_{2}1/\sqrt{\ln(1/\delta)}$. The first approximation
in (4) can be shown easily by relying on an alternate definition of sub-
gaussianity for mean zero random variables [45, Section 2.5.1]:
$\|X\|_{\psi_{2}}\approx_{O(1)}\inf\\{t>0:\mathbb{E}\exp(\lambda
X)\leq\exp(t^{2}\lambda^{2}),\ \forall_{\lambda>0}\\}$. Then, we can invoke
the independence between $\mathbf{x}_{i}$’s to obtain that a bound on the sub-
gaussian norm of $\tilde{\mathbf{x}}$ implies the same bound for $\mathbf{x}$
up to a constant factor:

$\displaystyle\mathbb{E}\,\exp\big{(}\lambda\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{x}\big{)}=\Big{(}\prod_{i=1}^{k}\mathbb{E}\,\exp\big{(}\lambda
r_{i}\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{x}_{i}\big{)}\Big{)}^{1/k}\\!=\Big{(}\mathbb{E}\exp\big{(}\sqrt{k}\lambda\mathbf{v}^{\scriptscriptstyle{\top}}\tilde{\mathbf{x}}\big{)}\Big{)}^{1/k}\\!\leq\exp(t^{2}\lambda^{2}k)^{1/k}=\exp(t^{2}\lambda^{2}).$

Having established that it is not practical to bound the sub-gaussian norm of
the $k$-gaussianized sample $\tilde{\mathbf{x}}$ directly, we now show that,
given our relaxed notion of sub-gaussianity defined through total variation
distance, the bound in our main result (Theorem 1) is optimal up to
logarithmic factors. To do so, we first point out the relationship between the
sub-gaussian norm of the vector $\mathbf{x}$ itself and the sub-gaussian norm
of its Euclidean length:

$\displaystyle\big{\|}\|\mathbf{x}\|_{2}\big{\|}_{\psi_{2}}\approx_{O(1)}\sqrt{d}+\|\mathbf{x}\|_{\psi_{2}}.$

This follows by observing that
$\mathbb{E}[\exp(\|\mathbf{x}\|^{2}/t^{2})]=\big{(}\mathbb{E}[\exp(M^{2}b_{1}/t^{2})]\big{)}^{d}$.
As a consequence, when $\mathbf{x}$ is sufficiently sparse (e.g.,
$M\geq\sqrt{d\log d}$), then
$\|\|\mathbf{x}\|_{2}\|_{\psi_{2}}\approx_{O(1)}\|\mathbf{x}\|_{\psi_{2}}$.

The below lower bound shows that, for the example distribution in (3), every
isotropic random vector $\mathbf{z}$ that is sufficiently close to
$\tilde{\mathbf{x}}$ has to have sub-gaussian norm at least
$(1+\|\mathbf{x}\|_{\psi_{2}}/\sqrt{k})\cdot c\|\mathbf{r}\|_{\psi_{2}}$. When
$M\geq\sqrt{d\log d}$, this nearly matches our upper bound of
$\tilde{O}(1+\|\|\mathbf{x}\|_{2}\|_{\psi_{2}}/\sqrt{k})\cdot\|\mathbf{r}\|_{\psi_{2}}$.

###### Theorem 4

Assume that $\mathbf{x}_{1},\mathbf{x}_{2},...$ are i.i.d. samples distributed
according to (3), and let $\tilde{\mathbf{x}}$ be the $k$-gaussianized sample.
Also, let $\delta\leq\frac{1}{2}\cdot\frac{dk}{3M^{2}+dk}$ and consider any
isotropic random vector $\mathbf{z}$ and covariance matrix
$\widetilde{\mathbf{\Sigma}}$ such that
$d_{{\mathrm{tv}}}\big{(}\,\widetilde{\mathbf{\Sigma}}^{-1/2}\,\tilde{\mathbf{x}},\,\mathbf{z}\,\big{)}\leq\delta$
and $\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{\Sigma}$. There is
an absolute constant $c>0$ such that:

$\displaystyle\|\mathbf{z}\|_{\psi_{2}}\geq\bigg{(}1+\frac{\|\mathbf{x}\|_{\psi_{2}}}{\sqrt{k}}\bigg{)}\cdot
c\|\mathbf{r}\|_{\psi_{2}}.$

##  6 Experiments

In this section, we aim to evaluate the degree of Algorithmic Gaussianization
for various sketching methods, using the sketched least squares task from
Section 3.1 as an example. In this task, we are given an $N\times d$ matrix
$\mathbf{A}$ and an $N$-dimensional vector $\mathbf{b}$ which define the
regression loss $L(\mathbf{w})=\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2}$. We
then use an $n\times N$ sketching matrix $\mathbf{S}$ to generate a sketched
least squares estimate
$\widehat{\mathbf{w}}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}}\|\mathbf{S}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2}$
of the exact solution $\mathbf{w}^{*}$. Our goal is to compare the normalized
expected approximation error
$\mathbb{E}[L(\widehat{\mathbf{w}})-L(\mathbf{w}^{*})]/L(\mathbf{w}^{*})$ of
the sketched estimate to the problem-independent expression $\frac{d}{n-d-1}$
achieved by a Gaussian embedding.

Figure 2: Comparison of expected approximation error for sketched least
squares on Libsvm regression datasets [8], with shading indicating standard
error of the mean.

In Figure 2, we plot the empirically estimated expected approximation error
with varying sketch size for several sketching techniques on four benchmark
datasets. We compare them against the expected error for Gaussian embeddings,
which is given by (1). Two of the sketching techniques, Subsampled Randomized
Hadamard Transform (SRHT, [2]) and Leverage Score Sparsified (LESS, [17])
embeddings, have strong theoretical guarantees for the expected error, and
have comparable nearly-linear time complexity.333Recall that the complexity of
sketched least squares with LESS is $O({\mathrm{nnz}}(\mathbf{A})\log
N+nd^{2})$, compared to $O(Nd\log N+nd^{2})$ for SRHT. So, the complexity of
LESS is better than SRHT when $\mathbf{A}$ is a sparse matrix, i.e., when
${\mathrm{nnz}}(\mathbf{A})\ll Nd$, and otherwise, the complexities are
comparable. However, there are some key differences: LESS has a precise non-
asymptotic guarantee given in this work (Theorem 2), whereas SRHT only has a
non-asymptotic upper bound [42] and a precise asymptotic guarantee under some
additional assumptions [20]. We compare these two sketches with two much
cheaper baselines: (1) uniform sampling, and (2) a simplified variant of LESS,
called LessUniform [14], which eliminates the leverage score approximation
preprocessing step (but still uses $d$ non-zeros per row). The overall cost
for both LessUniform and for uniform sampling is the same as the cost of
solving the sketched sub-problem, i.e., $O(nd^{2})$, where $n$ is the sketch
size. Our Theorems 1 and 2 can also be applied to LessUniform, but the Hanson-
Wright constant becomes dependent on the maximum leverage score of matrix
$\mathbf{A}$ (i.e., the coherence of $\mathbf{A}$, see Section 2.2).

From Figure 2, we first confirm that LESS enjoys a Gaussian-like problem-
independent expected approximation error that matches our theory. Remarkably,
we can verify this for all datasets and all sketch sizes, downto the precision
of our empirical mean estimates. Thus, these results suggest that the Gaussian
error estimate for LESS embeddings is even more accurate and broadly
applicable than promised by our theory. Next, we observe that for small sketch
sizes, SRHT does not always exhibit Gaussian-like error (see left two plots),
but as the sketch size increases, the gap decreases. This appears roughly in
line with the existing asymptotic theory for SRHTs. Finally, LessUniform
performs much better than uniform sampling, despite having the same time
complexity. The error for LessUniform appears Gaussian-like in three out of
four cases. In the case of dataset cpusmall, the problem exhibits very high
coherence, and as a result both uniform sampling and LessUniform perform
worse. This aligns with Theorem 1 applied to LessUniform, since the Hanson-
Wright constant scales with the largest leverage score of $\mathbf{A}$.

##  7 Conclusions

We provided the first general characterization of Algorithmic Gaussianization,
which refers to the phenomenon that many algorithmic techniques used to
construct small random representations (sketches) of large datasets produce
data samples that are more Gaussian-like than the original data. In our main
result, we showed that a simple averaging scheme (akin to a sparse sketch) can
be used to produce a sample that is nearly indistinguishable from a sub-
gaussian random design. We used this to show a reduction between a fast
sketching technique called LESS embeddings and a sub-gaussian random matrix
whose rows satisfy the classical Hanson-Wright inequality. Our lower bound
showed that the result is nearly tight up to logarithmic factors, and we
demonstrated how our techniques can be used to provide improved guarantees for
sketching based estimators in least squares, Lasso regression, and low-rank
approximation.

#### Acknowledgments.

Thanks to Edgar Dobriban and Michael Mahoney for valuable discussions
regarding this paper.

##  Appendix A More applications: Randomized SVD and low-distortion embeddings

In this section, we provide further examples of how Theorem 1 and Corollary 1
can be applied in combination with existing results for sub-gaussian random
designs. In each of these examples, we can extend existing results from sub-
gaussian embeddings (which are dense, and therefore not very efficient) to
analogous results for LESS embeddings (which are sparse, and thus can be
implemented much more efficiently).

###  A.1 Randomized SVD

An important application of randomized sketching is low-rank approximation,
where our goal is to estimate a small number of top principal directions of an
$N\times d$ data matrix $\mathbf{A}$. One of the most popular techniques in
this area is the Randomized SVD algorithm [26], where we construct the
approximation by projecting the dataset onto a subspace spanned by the rows of
the sketch $\mathbf{S}\mathbf{A}$, with $\mathbf{S}$ denoting the $n\times N$
sketching matrix, and $n\ll d$. The error of such an approximation is often
measured by the sum of squared lengths of the residuals from the projection:
$\|\mathbf{A}-\mathbf{A}\cdot\mathrm{Proj}_{\mathbf{S}\mathbf{A}}\|_{F}^{2}$,
where $\mathrm{Proj}_{\mathbf{S}\mathbf{A}}$ denotes the projection onto the
row-span of $\mathbf{S}\mathbf{A}$ and $\|\cdot\|_{F}$ is the Frobenius norm.
Prior work has shown strong approximation guarantees for Gaussian embeddings
in this context, such as the following, given by [26], comparing the expected
error to the best rank $k$ approximation, for some $k<n-2$:

$\displaystyle\mathbb{E}\Big{[}\big{\|}\mathbf{A}-\mathbf{A}\cdot\mathrm{Proj}_{\mathbf{S}\mathbf{A}}\big{\|}_{F}^{2}\Big{]}\leq\Big{(}1+\frac{k}{n-k-1}\Big{)}\cdot\min_{\mathbf{B}:\,\mathrm{rank}(\mathbf{B})=k}\|\mathbf{A}-\mathbf{B}\|_{F}^{2}.$

Those guarantees were later extended (and, in some regimes, improved) for sub-
gaussian embeddings [15] and Determinantal Point Processes [13]. The results
for Gaussian and sub-gaussian embeddings rely heavily on the Hanson-Wright
property of the rows of the sketching matrices. Thus, our techniques can be
used to extend them to LESS embeddings by applying the coupling argument based
on Corollary 1. Consequently, we can obtain bounds for the expected
approximation error of Randomized SVD with LESS embeddings that, in a certain
regime of small sketch sizes, nearly match those achieved by Gaussian
embeddings.

###### Corollary 3

Consider an $N\times d$ matrix $\mathbf{A}$ with stable rank
$r=\|\mathbf{A}\|_{F}^{2}/\|\mathbf{A}\|^{2}$. Given a LESS embedding
$\mathbf{S}$ of size $n\leq r/2$, and a Gaussian embedding
$\widetilde{\mathbf{S}}$ of the same size, we have for
$\epsilon=\tilde{O}(1/\sqrt{r})$:

$\displaystyle\mathbb{E}\Big{[}\big{\|}\mathbf{A}-\mathbf{A}\cdot\mathrm{Proj}_{\mathbf{S}\mathbf{A}}\big{\|}_{F}^{2}\Big{]}\approx_{1+\epsilon}\mathbb{E}\Big{[}\big{\|}\mathbf{A}-\mathbf{A}\cdot\mathrm{Proj}_{\widetilde{\mathbf{S}}\mathbf{A}}\big{\|}_{F}^{2}\Big{]}.$

Proof  We rely on Theorem 2 of [15], characterizing the expected approximation
error for a sub-gaussian embedding of size $n$ via an implicit analytic
formula:

$\displaystyle\mathbb{E}\Big{[}\big{\|}\mathbf{A}-\mathbf{A}\cdot\mathrm{Proj}_{\widetilde{\mathbf{S}}\mathbf{A}}\big{\|}_{F}^{2}\Big{]}\approx_{1+\epsilon}n\lambda_{n},\quad\text{where}\quad\lambda_{n}=f_{\mathbf{A}}^{-1}(n)\quad\text{for}\quad
f_{\mathbf{A}}(\lambda)=\mathrm{tr}\,\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}+\lambda\mathbf{I})^{-1}.$
(5)

Here, $f_{\mathbf{A}}^{-1}(n)$ denotes the function inverse of
$f_{\mathbf{A}}$ at $n$. The function $f_{\mathbf{A}}(\lambda)$ is known in
the literature as $\lambda$-statistical dimension, a notion of degrees of
freedom that arises in the analysis of ridge regression, e.g., see [16].
Naturally, their result applies to the Gaussian embedding
$\widetilde{\mathbf{S}}\mathbf{A}$, but also to any random design
$\mathbf{Z}\mathbf{\Sigma}^{1/2}$ (where
$\mathbf{\Sigma}=\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$) such that
the rows satisfy the Hanson-Wright inequality with a constant $K$ (with
$\epsilon$ having a polynomial dependence on $K$). So, we can apply our
reduction to LESS embeddings by coupling $\sqrt{n}\mathbf{S}\mathbf{A}$ with
$\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}$ for some
$\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{\Sigma}$, as in
Corollary 1. Letting $\mathrm{Er}(\mathbf{S}\mathbf{A})$ denote
$\|\mathbf{A}-\mathbf{A}\cdot\mathrm{Proj}_{\mathbf{S}\mathbf{A}}\|_{F}^{2}$
and $\mathcal{E}$ be the $(1-\delta)$-probability event where
$\sqrt{n}\mathbf{S}\mathbf{A}=\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}$, we
have:

$\displaystyle\big{|}\mathbb{E}\big{[}\mathrm{Er}(\mathbf{S}\mathbf{A})\big{]}-\mathbb{E}\big{[}\mathrm{Er}(\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2})\big{]}\big{|}$
$\displaystyle=\delta\cdot\big{|}\mathbb{E}\big{[}\mathrm{Er}(\mathbf{S}\mathbf{A})-\mathrm{Er}(\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2})\mid\neg\mathcal{E}\big{]}\big{|}\leq\delta\cdot\|\mathbf{A}\|_{F}^{2}.$

Since $\|\mathbf{A}\|_{F}^{2}\leq
r\cdot\|\mathbf{A}\|^{2}\leq\frac{r}{r-n}\cdot\min_{\mathbf{Z}}\mathrm{Er}(\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2})$,
setting $\delta\leq\epsilon/2$, we get
$\mathbb{E}\big{[}\mathrm{Er}(\mathbf{S}\mathbf{A})\big{]}\approx_{1+\epsilon}\mathbb{E}\big{[}\mathrm{Er}(\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2})\big{]}$.
Next, since $\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{\Sigma}$, we
can use the implicit error formula for sketching
$\widetilde{\mathbf{\Sigma}}^{1/2}$ with $\mathbf{Z}$, and relate that back to
the original problem as follows:

$\displaystyle\mathbb{E}\big{[}\mathrm{Er}(\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2})\big{]}\approx_{1+\delta}\mathbb{E}\big{[}\|\widetilde{\mathbf{\Sigma}}^{1/2}(\mathbf{I}-\mathrm{Proj}_{\mathbf{Z}\widetilde{\mathbf{\Sigma}}})\|_{F}^{2}\big{]}\approx_{1+\epsilon}n\tilde{\lambda}_{n},$

where $\tilde{\lambda}_{n}$ is defined as in (5), but with
$\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$ replaced by
$\widetilde{\mathbf{\Sigma}}$. To close the loop, we must relate the implicit
analytic expression based on $\widetilde{\mathbf{\Sigma}}$ back to
$\lambda_{n}$. This can be done easily by bounding the derivative of
$f_{\mathbf{A}}^{-1}$, showing that
$\tilde{\lambda}_{n}\approx_{1+O(\delta)}\lambda_{n}$. Now, it remains to
observe that since both the LESS embedding error
$\mathbb{E}\big{[}\mathrm{Er}(\mathbf{S}\mathbf{A})\big{]}$ and the Gaussian
embedding error
$\mathbb{E}\big{[}\mathrm{Er}(\widetilde{\mathbf{S}}\mathbf{A})\big{]}$ are
approximated by the same quantity $n\lambda_{n}$ up to a $1+O(\epsilon)$
factor, they are also approximated by each other, concluding the proof.

###  A.2 Low-distortion embeddings

Another important property of sub-gaussian random matrixes is that they can be
used to construct low-dimensional embeddings that preserve the geometry of a
high-dimensional space, such as the subspace embedding property and the
Johnson-Lindenstrauss property, which are central to many applications of
sketching (e.g., see [48]). We briefly mention some classical examples and
discuss how they can be extended to gaussianized matrices, and in particular,
to LESS embeddings.

We say that an $n\times N$ matrix $\mathbf{S}$ is an $\epsilon$-low-distortion
embedding for some set of points $\mathcal{X}\subseteq\mathbb{R}^{N}$ if there
is a fixed scalar $\alpha>0$ such that
$\|\alpha\mathbf{S}\mathbf{v}\|\approx_{1+\epsilon}\|\mathbf{v}\|$ for all
$\mathbf{v}\in\mathcal{X}$, i.e., the embedding approximately preserves the
Euclidean norm of $\mathbf{v}$. When $\mathbf{S}$ is a sketching matrix
applied to an $N\times d$ data matrix $\mathbf{A}$, we are typically
interested in vectors $\mathbf{v}$ from the column span of $\mathbf{A}$ (i.e.,
such that $\mathbf{v}=\mathbf{A}\mathbf{x}$ for some
$\mathbf{x}\in\mathbb{R}^{d}$), denoted by $\mathrm{span}(\mathbf{A})$.
Relying on Corollary 1 and standard properties of sub-gaussian random matrices
[45], we can establish the following low-distortion embedding properties for
LESS embeddings.

###### Corollary 4

Let $\mathbf{A}$ be $N\times d$, and $\mathbf{S}$ be a LESS embedding for
$\mathbf{A}$ of size $n$ with $d\log(nd/(\epsilon\delta))$ non-zeros per row.
Then, there is an absolute constant $C$ such that the following claims are
true for $\mathbf{S}$:

  1. 1.

_(Johnson-Lindenstrauss property)_ For any finite set
$\mathcal{X}\subseteq\mathrm{span}(\mathbf{A})$, if $n\geq
C\log(|\mathcal{X}|/\delta)/\epsilon^{2}$, then with probability $1-\delta$,
matrix $\mathbf{S}$ is an $\epsilon$-low-distortion embedding for
$\mathcal{X}$.

  2. 2.

_(Subspace embedding property)_ If $n\geq C(d+\log(1/\delta))/\epsilon^{2}$,
then with probability $1-\delta$, matrix $\mathbf{S}$ is an $\epsilon$-low-
distortion embedding for $\mathcal{X}=\mathrm{span}(\mathbf{A})$.

Proof  Whenever, $\mathcal{X}\subseteq\mathrm{span}(\mathbf{A})$, the low-
distortion property can be formulated as a property of the sketch
$\mathbf{S}\mathbf{A}$, treated as a linear transformation: for any
$\mathbf{v}\in\mathrm{span}(\mathbf{A})$, we have
$\mathbf{v}=\mathbf{A}\mathbf{x}$ for some $\mathbf{x}\in\mathbb{R}^{d}$, and
we need to ensure that
$\|\alpha\mathbf{S}\mathbf{A}\mathbf{x}\|\approx_{1+\epsilon}\|\mathbf{A}\mathbf{x}\|$.
Here, we can once again rely on Corollary 1, coupling
$\sqrt{n}\mathbf{S}\mathbf{A}$ with a sub-gaussian design
$\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}$ for some
$\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$.
If $\mathbf{Z}$ is an $\epsilon$-low-distortion embedding for
$\widetilde{\mathcal{X}}=\\{\widetilde{\mathbf{\Sigma}}^{1/2}\mathbf{x}:\
\mathbf{A}\mathbf{x}\in\mathcal{X}\\}$ with constant $\alpha$, then with
probability $1-\delta$ (specifically, when
$\sqrt{n}\mathbf{S}\mathbf{A}=\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}$),
for any $\mathbf{A}\mathbf{x}\in\mathcal{X}$, we have:

$\displaystyle\|\alpha\sqrt{n}\mathbf{S}\mathbf{A}\mathbf{x}\|=\|\alpha\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}\mathbf{x}\|\approx_{1+\epsilon}\|\widetilde{\mathbf{\Sigma}}^{1/2}\mathbf{x}\|\approx_{1+\delta}\|\mathbf{A}\mathbf{x}\|,$

where in the last step we used that
$\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A}$.
Thus, as long as $\delta\leq\epsilon$, we reduced the problem to showing the
corresponding property for $\mathbf{Z}$. Both the Johnson-Lindenstrauss and
subspace embedding properties can be shown for sub-gaussian embeddings (with
sizes as given in the statement) using standard arguments (e.g., see the proof
of Theorem 4.6.1 in [45]), thus concluding the proof.
We note that low-distortion embedding properties have been shown for other
sketching operators, including SRHTs [43] and other sparse sketching matrices
like CountSketch [10] and SJLTs [39]. These results often circumvent sub-
gaussian analysis by relying on matrix concentration inequalities such as the
matrix Bernstein inequality (Lemma 3) or other more involved arguments,
however this inevitably incurs additional overhead factors in the required
sketch size $n$. For example, an SRHT requires sketch size of at least $n\geq
O(d\log(d)/\epsilon^{2})$ to establish the subspace embedding property,
instead of $O(d/\epsilon^{2})$ shown above (the additional log factor is a
well-known limitation of the matrix concentration-style analysis). The
$O(d\log(d)/\epsilon^{2})$ guarantee was also previously shown for LESS
embeddings with $d$ non-zeros per row, also using matrix concentration
inequalities instead of sub-gaussian concentration [17, Lemma 12]. Thus, the
above corollary suggests that, when requiring a subspace embedding, we can
trade an additional log-factor in the sketch size for a log-factor in the
density of a LESS embedding matrix. This could be a useful trade-off when we
wish to minimize the storage space of the sketched data.

##  Appendix B Notation and preliminaries for the proofs

#### Notation.

We say that $a\approx_{\alpha}b$ for $\alpha\geq 1$, if $b/\alpha\leq
a\leq\alpha b$. We analogously define $\mathbf{A}\approx_{\alpha}\mathbf{B}$
for positive semidefinite matrices using the Loewner ordering. We use
$\|\mathbf{A}\|_{F}=\sqrt{\mathrm{tr}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})}$
to denote the Frobenius norm and $\|\mathbf{A}\|$ to denote the spectral norm.
Also, we let $a\lesssim b$ denote that there exists an absolute constant $C>0$
such that $a\leq Cb$, and $a={\mathrm{poly}}(b,c)$ means that $a$ is bounded
by a polynomial function of $b$ and $c$. Moreover, we use $\tilde{O}(\cdot)$
to refer to big-O notation where polylogarithmic terms are ignored. For random
variables/vectors $X$ and $Y$ defined over the same domain with measures $\mu$
and $\nu$, respectively, we define the total variation distance between them
as $d_{{\mathrm{tv}}}(X,Y)=\sup_{E\in\mathcal{B}}|\mu(E)-\nu(E)|$, where
$\mathcal{B}$ denotes all measurable events. Note that the total variation
distance can also be defined as the infimum over $\delta$ such that there
exists a coupling between $X$ and $Y$ for which $\mathrm{Pr}(X\neq Y)=\delta$.
Finally, we define the sub-gaussian Orlicz norm as:
$\|X\|_{\psi_{2}}=\inf\\{t>0:\,\mathbb{E}\,\exp(X^{2}/t^{2})\leq 2\\}$.

#### Matrix concentration inequalities.

In the proof of our main results, we use the following versions of the matrix
Bernstein concentration inequality for the sums of independent symmetric
random matrices.

######  Lemma 2 (Sub-exponential matrix Bernstein)

[44, Theorem 6.2] For $i=1,2,...$, consider a finite sequence $\mathbf{M}_{i}$
of $d\times d$ independent and symmetric random matrices such that

$\mathbb{E}[\mathbf{M}_{i}]=\mathbf{0},\quad\mathbb{E}[\mathbf{M}_{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{M}_{i}\Big{)}\geq
t\bigg{\\}}$ $\displaystyle\leq
d\cdot\exp\left(\frac{-t^{2}/2}{\sigma^{2}+Rt}\right).$

######  Lemma 3 (Bounded matrix Bernstein)

[44, Theorem 6.1] For $i=1,2,...$, consider a finite sequence $\mathbf{M}_{i}$
of $d\times d$ independent and symmetric random matrices such that

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

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

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

#### Linear algebraic identities.

Our proof of the least squares approximation guarantee for LESS embeddings
relies on the following standard rank-one update formula for the matrix
inverse.

######  Lemma 4 (Sherman-Morrison formula)

For a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ and
$\mathbf{u},\mathbf{v}\in\mathbb{R}^{n}$ such that both $\mathbf{A}$ and
$\mathbf{A}+\mathbf{u}\mathbf{v}^{\scriptscriptstyle{\top}}$ are invertible,
we have:

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

##  Appendix C Sub-gaussian limit theorem: Proof of Theorem 1

In this section, we prove our main result, Theorem 1, showing that a
$k$-gaussianized sample
$\tilde{\mathbf{x}}=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}r_{i}\mathbf{x}_{i}$ is
close in total variation distance to a random vector that satisfies the
Hanson-Wright inequality with a small constant.

We start by using the matrix Bernstein inequality to verify how well the
sample covariance
$\frac{1}{k}\sum_{i=1}^{k}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}$
can be bounded in terms of the true covariance $\mathbf{\Sigma}$. For this we
establish the following lemma, proven in Appendix G. The lemma is used to show
the matrix moment condition required by the sub-exponential matrix Bernstein
(Lemma 2).

###### Lemma 5

There is an absolute constant $C$ such that any $d$-dimensional random vector
$\mathbf{x}$ with covariance
$\mathbb{E}[\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}]=\mathbf{\Sigma}$,
such that $\|\mathbf{\Sigma}^{-1/2}\mathbf{x}\|$ is $M$-sub-gaussian,
satisfies the following sub-exponential matrix moment bound:

$\displaystyle\bigg{\|}\mathbb{E}\Big{[}\Big{(}\mathbf{\Sigma}^{-1/2}\mathbf{x}\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{\Sigma}^{-1/2}-\mathbf{I}\Big{)}^{p}\Big{]}\bigg{\|}\leq\big{(}CM^{2}(p+\log
M)\big{)}^{p-1}.$

We can now use Lemma 5 in conjunction with the sub-exponential matrix
Bernstein inequality (Lemma 2) by setting
$\mathbf{M}_{i}=\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{\Sigma}^{-1/2}-\mathbf{I}$.
From Lemma 5, we can set $R=CM^{2}\log M$ and
$\mathbf{A}_{i}^{2}=R\mathbf{I}$, with $\sigma^{2}=kR$ (appropriately
adjusting constant $C$), concluding that:

$\displaystyle\mathrm{Pr}\bigg{\\{}\lambda_{\max}\Bigg{(}\sum_{i=1}^{k}\mathbf{M}_{i}\Big{)}\geq
kt\Big{\\}}\leq d\cdot\exp\Big{(}\frac{-t^{2}k}{2(1+t)CM^{2}\log M}\Big{)}.$

Now, define the $k\times d$ matrix $\mathbf{U}$ with $i$th row
$\frac{1}{\sqrt{k}}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{\Sigma}^{-1/2}$.
The above concentration inequality implies that if $k\geq
4CM^{2}\log^{2}(M/\delta)/t$ for $t\geq 1$, then
$\|\mathbf{U}\|^{2}=\|\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}\|\leq
1+t$ with probability $1-\delta$. We mention that, if we assumed
$\|\mathbf{\Sigma}^{-1/2}\mathbf{x}\|\leq M$ almost surely, as in Remark 1,
then we can use the bounded matrix Bernstein (Lemma 3) to obtain a slightly
sharper guarantee of $k\geq O(M^{2}\log(d/\delta)/t)$, without relying on
Lemma 5.

We next define the random variable $\mathbf{z}$ and covariance matrix
$\widetilde{\mathbf{\Sigma}}$ discussed in the theorem:

$\displaystyle\mathbf{z}$
$\displaystyle=\mathbf{1}_{\mathcal{E}}\cdot\widetilde{\mathbf{\Sigma}}^{-1/2}\tilde{\mathbf{x}}+\mathbf{1}_{\neg\mathcal{E}}\cdot\alpha\mathbf{V}^{\scriptscriptstyle{\top}}\mathbf{r},\qquad\widetilde{\mathbf{\Sigma}}=\frac{1}{\mathrm{Pr}\\{\mathcal{E}\\}}\,\mathbb{E}[\mathbf{1}_{\mathcal{E}}\tilde{\mathbf{x}}\tilde{\mathbf{x}}^{\scriptscriptstyle{\top}}],$
$\displaystyle\text{for}\qquad\mathcal{E}$
$\displaystyle=\big{[}\|\mathbf{U}\|\leq
L\big{]},\qquad\alpha=\max\Big{\\{}1,\sqrt{\frac{d}{k}}\,\Big{\\}},\qquad\mathbf{V}=\mathbf{G}(\mathbf{G}^{\scriptscriptstyle{\top}}\mathbf{G})^{\dagger/2},$

where $\mathbf{r}$ is the vector of $R$-sub-gaussian random variables $r_{i}$
from the definition of $\tilde{\mathbf{x}}$; $\mathbf{G}$ denotes a $k\times
d$ matrix with i.i.d. Gaussian entries and $(\cdot)^{\dagger}$ is the Moore-
Penrose pseudoinverse (so that either
$\mathbf{V}^{\scriptscriptstyle{\top}}\mathbf{V}=\mathbf{I}$ or
$\mathbf{V}\mathbf{V}^{\scriptscriptstyle{\top}}=\mathbf{I}$); and
$\mathcal{E}$ denotes an event in the probability space of
$\tilde{\mathbf{x}}$, with $\mathbf{1}_{\mathcal{E}}$ being the characteristic
function of $\mathcal{E}$. Here, $L=O(1+M\log(M/\delta)/\sqrt{k})$ is chosen
so that $\mathrm{Pr}\\{\neg\mathcal{E}\\}\leq\delta/L^{2}$.444If
$\|\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}\|\leq M$ a.s., as in Remark 1, then we
can let $L=O(1+M\sqrt{\log(M/\delta)/k})$ by using bounded matrix Bernstein
(Lemma 3). We can do this by simply adjusting the constants, because
$M\log(ML/\delta)=O(M\log(M/\delta))$. Also, we can easily make sure that
$L\geq\max\\{2,\alpha\\}$ (recall that $\alpha\leq 1+\sqrt{d/k}\leq
1+M/\sqrt{k}$). Note that $\mathbf{z}$ is an isotropic random vector since:

$\displaystyle\mathbb{E}[\mathbf{z}\mathbf{z}^{\scriptscriptstyle{\top}}]=\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbb{E}[\mathbf{1}_{\mathcal{E}}\tilde{\mathbf{x}}\tilde{\mathbf{x}}^{\scriptscriptstyle{\top}}]\widetilde{\mathbf{\Sigma}}^{-1/2}+\mathrm{Pr}\\{\neg\mathcal{E}\\}\cdot\mathbb{E}[\alpha^{2}\mathbf{V}^{\scriptscriptstyle{\top}}\mathbf{V}]=\mathrm{Pr}\\{\mathcal{E}\\}\mathbf{I}+\mathrm{Pr}\\{\neg\mathcal{E}\\}\,\mathbf{I}=\mathbf{I}.$

Also, note that we can write
$\mathbf{z}=\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\mathbf{r}$,
where

$\displaystyle\widetilde{\mathbf{U}}=\mathbf{1}_{\mathcal{E}}\cdot\mathbf{U}\mathbf{\Sigma}^{1/2}\widetilde{\mathbf{\Sigma}}^{-1/2}+\mathbf{1}_{\neg\mathcal{E}}\cdot\alpha\mathbf{V},$

and since $\mathbf{r}$ is mean zero and independent of
$\widetilde{\mathbf{U}}$, then $\mathbf{z}$ is also mean zero.

Next we show that $\widetilde{\mathbf{\Sigma}}$ approximates
$\mathbf{\Sigma}$. First, by definition we immediately have
$\widetilde{\mathbf{\Sigma}}\preceq\frac{1}{1-\delta/2}\mathbf{\Sigma}\preceq(1+\delta)\mathbf{\Sigma}$.
Next, observe that we have:

$\displaystyle\|\mathbf{\Sigma}^{-\frac{1}{2}}(\mathbf{\Sigma}-\widetilde{\mathbf{\Sigma}})\mathbf{\Sigma}^{-\frac{1}{2}}\|$
$\displaystyle=\|\mathbf{\Sigma}^{-\frac{1}{2}}\mathbb{E}\big{[}\tilde{\mathbf{x}}\tilde{\mathbf{x}}^{\scriptscriptstyle{\top}}\cdot\mathbf{1}_{\neg\mathcal{E}}\big{]}\mathbf{\Sigma}^{-\frac{1}{2}}\|=\big{\|}\mathbb{E}\big{[}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{r}\mathbf{r}^{\scriptscriptstyle{\top}}\mathbf{U}\cdot\mathbf{1}_{\neg\mathcal{E}}\big{]}\big{\|}$
$\displaystyle=\|\mathbb{E}[\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}\cdot\mathbf{1}_{\neg\mathcal{E}}]\|\leq\mathbb{E}\big{[}\|\mathbf{U}\|^{2}\cdot\mathbf{1}_{\neg\mathcal{E}}\big{]}$
$\displaystyle\overset{(*)}{=}\int_{0}^{\infty}\mathrm{Pr}\big{\\{}\|\mathbf{U}\|^{2}\cdot\mathbf{1}_{\neg\mathcal{E}}>x\big{\\}}dx$
$\displaystyle\leq\frac{\delta}{L^{2}}\cdot
L^{2}+\int_{L^{2}}^{\infty}\mathrm{Pr}\big{\\{}\|\mathbf{U}\|^{2}>x\big{\\}}dx,$

where $(*)$ is the integral formula for the expectation of a non-negative
random variable via its cdf, and the last step uses the observation that we
have
$\mathrm{Pr}\\{\|\mathbf{U}\|^{2}\cdot\mathbf{1}_{\neg\mathcal{E}}>x\\}\leq\mathrm{Pr}\\{\neg\mathcal{E}\\}\leq\delta/L^{2}$.
To bound the integral, note that for $x\geq L^{2}$:

$\displaystyle\mathrm{Pr}\big{\\{}\|\mathbf{U}\|^{2}>x\big{\\}}\leq
d\cdot\exp\Big{(}\frac{-(x-1)k}{4CM^{2}\log
d}\Big{)}\leq(\delta/L^{2})^{(x-1)/L^{2}}.$

Thus, using the formula
$\int(\delta/L^{2})^{x/L^{2}}dx=-L^{2}(\delta/L^{2})^{x/L^{2}}/\ln(L^{2}/\delta)$,
we obtain:

$\displaystyle\int_{L^{2}}^{\infty}\mathrm{Pr}\big{\\{}\|\mathbf{U}\|^{2}>x\big{\\}}dx$
$\displaystyle\leq(\delta/L^{2})^{-1/L^{2}}\int_{L^{2}}^{\infty}(\delta/L^{2})^{x/L^{2}}=(\delta/L^{2})^{-1/L^{2}}\cdot
L^{2}(\delta/L^{2})/\ln(L^{2}/\delta)\leq\delta.$

We conclude that
$\|\mathbf{\Sigma}^{-\frac{1}{2}}(\mathbf{\Sigma}-\widetilde{\mathbf{\Sigma}})\mathbf{\Sigma}^{-\frac{1}{2}}\|\leq\delta+\delta\leq
2\delta$, obtaining:

$\displaystyle(1-2\delta)\cdot\mathbf{\Sigma}\preceq\widetilde{\mathbf{\Sigma}}\preceq(1+\delta)\mathbf{\Sigma},$

so $\widetilde{\mathbf{\Sigma}}\approx_{1+O(\delta)}\mathbf{\Sigma}$. Note
that the above also implies that
$\|\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}^{1/2}\|^{2}=\|\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}\widetilde{\mathbf{\Sigma}}^{-1/2}\|\leq\frac{1}{1-2\delta}$.
As a consequence we have:

$\displaystyle\|\widetilde{\mathbf{U}}\|^{2}$
$\displaystyle=\mathbf{1}_{\mathcal{E}}\cdot\|\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}^{1/2}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{\Sigma}^{1/2}\widetilde{\mathbf{\Sigma}}^{-1/2}\|+\mathbf{1}_{\neg\mathcal{E}}\cdot\|\alpha^{2}\mathbf{V}^{\scriptscriptstyle{\top}}\mathbf{V}\|$
$\displaystyle\leq\max\big{\\{}\|\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}^{1/2}\|^{2}\|\mathbf{1}_{\mathcal{E}}\mathbf{U}\|^{2},\alpha^{2}\big{\\}}\leq\frac{L^{2}}{1-2\delta}\leq
2L^{2}.$

We are now ready to establish the Hanson-Wright inequality for $\mathbf{z}$.
Specifically, consider a $d\times d$ psd matrix $\mathbf{B}$. We study the
concentration of the quadratic form
$\mathbf{z}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{z}$ around its mean
$\mathrm{tr}\mathbf{B}$. We start with the following decomposition:

$\displaystyle|\mathbf{z}^{\scriptscriptstyle{\top}}\mathbf{B}\mathbf{z}-\mathrm{tr}\mathbf{B}|$
$\displaystyle=|\mathbf{r}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\mathbf{r}-\mathrm{tr}\mathbf{B}|\leq|\mathbf{r}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\mathbf{r}-\mathrm{tr}(\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}})|+|\mathrm{tr}(\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}})-\mathrm{tr}\mathbf{B}|.$
(6)

Since $\widetilde{\mathbf{U}}$ and $\mathbf{r}$ are independent, and we have
$\mathbb{E}[\mathbf{r}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\mathbf{r}]=\mathrm{tr}(\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}})$,
we can show concentration for the first term in (6) by using the classical
Hanson-Wright inequality (Lemma 1) applied to the random vector $\mathbf{r}$
and the matrix
$\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}$:

$\displaystyle\mathrm{Pr}\Big{\\{}|\mathbf{r}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\mathbf{r}-\mathrm{tr}(\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}})|\geq
t\mid\widetilde{\mathbf{U}}\Big{\\}}$ $\displaystyle\leq
2\exp\Big{(}-c\min\Big{\\{}\frac{t^{2}}{R^{4}\|\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}\|_{F}^{2}},\frac{t}{R^{2}\|\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}}\|}\Big{\\}}\Big{)}$
$\displaystyle\leq
2\exp\Big{(}-c\min\Big{\\{}\frac{t^{2}}{4L^{4}R^{4}\|\mathbf{B}\|_{F}^{2}},\frac{t}{2L^{2}R^{2}\|\mathbf{B}\|}\Big{\\}}\Big{)},$

where we used that $\|\widetilde{\mathbf{U}}\|^{2}\leq 2L^{2}$. Next, to show
concentration for the second term in (6), we decompose it further as follows,
letting
$\widetilde{\mathbf{B}}=\mathbf{\Sigma}^{1/2}\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbf{B}\widetilde{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}^{1/2}$:

$\displaystyle|\mathrm{tr}(\widetilde{\mathbf{U}}\mathbf{B}\widetilde{\mathbf{U}}^{\scriptscriptstyle{\top}})-\mathrm{tr}\mathbf{B}|$
$\displaystyle\leq\mathbf{1}_{\mathcal{E}}\cdot|\mathrm{tr}(\mathbf{U}\widetilde{\mathbf{B}}\mathbf{U}^{\scriptscriptstyle{\top}})-\mathrm{tr}\mathbf{B}|+\mathbf{1}_{\neg\mathcal{E}}\cdot|\alpha^{2}\mathrm{tr}\mathbf{V}\mathbf{B}\mathbf{V}^{\scriptscriptstyle{\top}}-\mathrm{tr}\mathbf{B}|$
$\displaystyle\leq|\mathrm{tr}(\mathbf{U}\widetilde{\mathbf{B}}\mathbf{U}^{\scriptscriptstyle{\top}})-\mathrm{tr}\widetilde{\mathbf{B}}|+|\mathrm{tr}\widetilde{\mathbf{B}}-\mathrm{tr}\mathbf{B}|+|\alpha^{2}\mathrm{tr}\mathbf{V}\mathbf{B}\mathbf{V}^{\scriptscriptstyle{\top}}-\mathrm{tr}\mathbf{B}|.$
(7)

Note that
$\mathrm{tr}(\mathbf{U}\widetilde{\mathbf{B}}\mathbf{U})-\mathrm{tr}\widetilde{\mathbf{B}}=\frac{1}{k}\sum_{i=1}^{k}X_{i}$
for
$X_{i}=\|\widetilde{\mathbf{B}}^{1/2}\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}\|^{2}-\mathrm{tr}\widetilde{\mathbf{B}}$.
We are going to again use Bernstein’s inequality, this time the scalar
version, i.e., Lemma 2 with $d=1$.555Here, again, if we have the boundedness
assumption from Remark 1, then we can use bounded Bernstein (scalar version of
Lemma 3). Note that we have $\mathbb{E}[X_{i}]=0$, and also in Lemma 10 (along
the same lines as Lemma 5, see Appendix G) we show that:

$\displaystyle\mathbb{E}[|X_{i}|^{p}]\leq(CM^{2}(p+\log
M))^{p-1}\|\widetilde{\mathbf{B}}\|^{p-2}\|\widetilde{\mathbf{B}}\|_{F}^{2}.$

Using Lemma 2 with $R=O(\|\widetilde{\mathbf{B}}\|M^{2}\log M)$ and
$\sigma^{2}=k\cdot O(\|\widetilde{\mathbf{B}}\|_{F}^{2}M^{2}\log M)$, we
obtain (for some absolute constants $c,c^{\prime}>0$):

$\displaystyle\mathrm{Pr}\Big{\\{}|\mathrm{tr}(\mathbf{U}\widetilde{\mathbf{B}}\mathbf{U}^{\scriptscriptstyle{\top}})-\mathrm{tr}\widetilde{\mathbf{B}}|\geq
t\Big{\\}}$ $\displaystyle\leq
2\exp\Big{(}-c\min\Big{\\{}\frac{t^{2}k}{\|\widetilde{\mathbf{B}}\|_{F}^{2}M^{2}\log
M},\frac{tk}{\|\widetilde{\mathbf{B}}\|M^{2}\log M}\Big{\\}}\Big{)}$
$\displaystyle\leq
2\exp\Big{(}-c^{\prime}\min\Big{\\{}\frac{t^{2}}{L^{2}\|\mathbf{B}\|_{F}^{2}},\frac{t}{L^{2}\|\mathbf{B}\|}\Big{\\}}\Big{)},$

where we used the facts that $\|\widetilde{\mathbf{B}}\|\leq 2\|\mathbf{B}\|$
and $L^{2}\geq M^{2}\log(M)/k$. We now turn to the second term in (7). Note
that this is in fact a deterministic quantity which can be bounded as follows:
$|\mathrm{tr}\widetilde{\mathbf{B}}-\mathrm{tr}\mathbf{B}|\leq
4\delta\cdot\mathrm{tr}\mathbf{B}$. Now, assuming that $\delta\leq 1/(4d)$
(which is w.l.o.g. by adjusting the constants), we have that for any $0\leq
t\leq 4\delta\mathrm{tr}\mathbf{B}$:

$\displaystyle\min\Big{\\{}\frac{t^{2}}{\|\mathbf{B}\|_{F}^{2}},\frac{t}{\|\mathbf{B}\|}\Big{\\}}$
$\displaystyle\leq\min\Big{\\{}\frac{(\mathrm{tr}\mathbf{B})^{2}}{d^{2}\|\mathbf{B}\|_{F}^{2}},\frac{\mathrm{tr}\mathbf{B}}{d\|\mathbf{B}\|}\Big{\\}}\leq\min\Big{\\{}\frac{\|\mathbf{B}\|^{2}}{\|\mathbf{B}\|_{F}^{2}},\frac{\|\mathbf{B}\|}{\|\mathbf{B}\|}\Big{\\}}\leq
1,$

where we used that $\mathrm{tr}\mathbf{B}/d\leq\|\mathbf{B}\|$. This means
that
$\mathrm{Pr}\\{\|\mathrm{tr}\widetilde{\mathbf{B}}-\mathrm{tr}\mathbf{B}|\geq
t\\}\leq\mathbf{1}_{[t\leq 4\delta\cdot\mathrm{tr}\mathbf{B}]}\leq
2\exp(-c\min\\{\frac{t^{2}}{\|\mathbf{B}\|_{F}^{2}},\frac{t}{\|\mathbf{B}\|}\\})$,
so the desired concentration inequality trivially holds. Finally, it remains
to establish the concentration for
$|\alpha^{2}\mathrm{tr}(\mathbf{V}\mathbf{B}\mathbf{V}^{\scriptscriptstyle{\top}})-\mathrm{tr}\mathbf{B}|$,
the final term in (7). If $k\geq d$, then
$\mathbf{V}^{\scriptscriptstyle{\top}}\mathbf{V}=\mathbf{I}$ and
$\alpha^{2}=1$, so the term is $0$. Now, suppose that $k<d$. In this case,
$\mathbf{V}^{\scriptscriptstyle{\top}}\mathbf{V}$ is a projection onto a
uniformly random $k$-dimensional subspace of $\mathbb{R}^{d}$. Here, we will
use a simple form of the Johnson-Lindenstrauss lemma for uniformly random
projections [45, Lemma 5.3.2], which states that for any fixed unit vector
$\mathbf{b}\in\mathbb{R}^{d}$:

$\displaystyle\mathrm{Pr}\Big{\\{}\big{|}\|\alpha\mathbf{V}\mathbf{b}\|-1\big{|}\geq\epsilon\Big{\\}}\leq
2\exp(-c\epsilon^{2}k).$

Let
$\mathbf{B}=\sum_{i=1}^{d}\lambda_{i}\mathbf{b}_{i}\mathbf{b}_{i}^{\scriptscriptstyle{\top}}$
be the eigendecomposition of $\mathbf{B}$. The Johnson-Lindenstrauss
concentration lemma implies that with probability $1-2d\exp(-c\epsilon^{2}k)$
we have $\max_{i}|\|\alpha\mathbf{V}\mathbf{b}_{i}\|-1|\leq\epsilon$, which in
turn implies that:

$\displaystyle|\alpha^{2}\mathrm{tr}(\mathbf{V}\mathbf{B}\mathbf{V}^{\scriptscriptstyle{\top}})-\mathrm{tr}\mathbf{B}|\leq\sum_{i=1}^{d}\lambda_{i}\big{|}\|\alpha\mathbf{V}\mathbf{b}_{i}\|^{2}-1\big{|}\leq\mathrm{tr}\mathbf{B}\cdot\max_{i}\big{|}\|\alpha\mathbf{V}\mathbf{b}_{i}\|^{2}-1\big{|}\leq\epsilon(\epsilon+2)\cdot\mathrm{tr}\mathbf{B}.$

Setting $t=\epsilon(\epsilon+2)\mathrm{tr}\mathbf{B}$ and solving for
$\epsilon$, we convert this to a concentration inequality:

$\displaystyle\mathrm{Pr}\Big{\\{}|\alpha^{2}\mathrm{tr}(\mathbf{V}\mathbf{B}\mathbf{V}^{\scriptscriptstyle{\top}})-\mathrm{tr}\mathbf{B}|\geq
t\Big{\\}}$ $\displaystyle\leq
2d\exp\Big{(}-c\min\Big{\\{}\frac{t^{2}k}{\mathrm{tr}\mathbf{B}},\frac{tk}{\mathrm{tr}\mathbf{B}}\Big{\\}}\Big{)}$
$\displaystyle\leq
2\exp\Big{(}-c\min\Big{\\{}\frac{t^{2}k}{\|\mathbf{B}\|d\log
d},\frac{tk}{\|\mathbf{B}\|d\log d}\Big{\\}}\Big{)}$ $\displaystyle\leq
2\exp\Big{(}-c\min\Big{\\{}\frac{t^{2}}{L^{2}\|\mathbf{B}\|_{F}^{2}},\frac{t}{L^{2}\|\mathbf{B}\|}\Big{\\}}\Big{)}.$

Putting everything together, we combine the four inequalities to obtain that
for any psd matrix $\mathbf{B}$:

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

We can easily extend this to arbitrary matrices $\mathbf{B}$. First, observe
that if $\mathbf{B}$ is symmetric, then we can write it as
$\mathbf{B}=\mathbf{B}_{+}-\mathbf{B}_{-}$ where $\mathbf{B}_{+}$ and
$\mathbf{B}_{-}$ are psd, and it suffices to apply the concentration
inequality to both $\mathbf{B}_{+}$ and $\mathbf{B}_{-}$. Next, consider an
arbitrary $\mathbf{B}$. Then, it suffices to apply the result to
$(\mathbf{B}+\mathbf{B}^{\scriptscriptstyle{\top}})/2$, which is symmetric.

##  Appendix D Sketched least squares: Proof of Theorem 2

In this section, we prove Theorem 2, giving a precise estimate for the
expected approximation error for sketched least squares with LESS embeddings.
Recall that we are given an $N\times d$ data matrix $\mathbf{A}$ and an
$N$-dimensional vector $\mathbf{b}$. We let
$L(\mathbf{w})=\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2}$ and
$\mathbf{w}^{*}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}}L(\mathbf{w})$.
Given an $n\times N$ LESS embedding matrix $\mathbf{S}$, we define the
sketched least squares estimator as:

$\displaystyle\widehat{\mathbf{w}}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}\in\mathcal{R}^{d}}\|\mathbf{S}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2},$

where $\mathcal{R}=[-D,D]$ is only used to protect against the corner cases
where $\mathbf{S}\mathbf{A}$ is very ill-conditioned, which arise with a small
but non-zero probability when computing the exact expectation of the loss. In
the proof, we will use our reduction between a LESS embedding and a sub-
gaussian design to define an event that holds with a high $1-\delta$
probability (Lemma 7). This event ensures certain regularity conditions on the
sketch, such as near-uniformity of the leverage scores, and avoids the
aforementioned corner cases. We then proceed to analyze the conditional
expectation, conditioned on the high-probability event. To obtain the final
unconditional expectation, we must choose $\delta$ sufficiently small to
absorb the worst-case loss that might occur if the high-probability event
fails. This is where we incur a dependence on the range diameter $D$, which is
in general unavoidable when computing the expected loss, but it only affects
the logarithmic factors in the result.

Let us briefly comment on those logarithmic dependencies. If we let the LESS
embedding have $d$ non-zeros per row, as recommended by [17], then our
reduction to a sub-gaussian sketch (which is used in Lemma 7) incurs a
logarithmic dependence in the Hanson-Wright constant, i.e.,
$K=O(\sqrt{\log(nd/\delta)})$, and also, the lower bound on the sketch size
required for the claim becomes $n\geq O(d\log(d/\delta))$. This also results
in a polylogarithmic factor in the error $\epsilon$. We note that these
logarithmic terms can be avoided if we slightly increase the density of the
sketch, from $d$ to $d\log(nd/\delta)$ non-zeros per row, following the
recommendation in Remark 3. In that case, the result can be established for
$n\geq O(d)$ and $\epsilon=O(1/\sqrt{d})$.

Instead of computing the expectation of $L(\widehat{\mathbf{w}})$ directly, we
will define an unbiased estimate of the expected loss via leave-one-out cross-
validation (CV):

$\displaystyle
L_{{\mathrm{CV}}}=\sum_{i=1}^{n}\big{(}\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\widehat{\mathbf{w}}^{-i}-\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{b}\big{)}^{2},\qquad\text{for}\quad\widehat{\mathbf{w}}_{-i}=\operatorname*{\mathop{\mathrm{argmin}}}_{\mathbf{w}\in\mathcal{R}^{d}}\|\mathbf{S}_{-i}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2},$

where $\mathbf{s}_{i}^{\scriptscriptstyle{\top}}$ is the $i$th row of
$\mathbf{S}$ and $\mathbf{S}_{-i}$ denotes $\mathbf{S}$ without the $i$th row.
We use the following standard “shortcut formula” for the CV estimate, which
holds as long as the scalar range $\mathcal{R}=[-D,D]$ is large enough so that
each $\widehat{\mathbf{w}}^{-i}$ is the same as the unconstrained least
squares solution (we include the proof in Appendix H).

###### Lemma 6

Suppose that
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{S}_{-i}\mathbf{A})^{-1}\preceq
O(1)\cdot(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ for all $i$.
Then, for sufficiently large $D$, the sketched estimator is given by
$\widehat{\mathbf{w}}=(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{b}$,
and the leave-one-out cross-validation loss can be computed as:

$\displaystyle
L_{{\mathrm{CV}}}=\sum_{i=1}^{n}\bigg{(}\frac{\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\widehat{\mathbf{w}}-\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{b}}{1-\ell_{i}(\mathbf{S}\mathbf{A})}\bigg{)}^{2},$

where
$\ell_{i}(\mathbf{S}\mathbf{A})=\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}$
is the $i$th leverage score of the sketch.

Crucially, the CV loss is an unbiased estimate of the expected loss based on
the sketch of size $n-1$:

$\displaystyle\mathbb{E}[L_{{\mathrm{CV}}}]=\mathbb{E}[L(\widehat{\mathbf{w}}^{-n})],$

Thus, for the rest of the proof, we will actually estimate
$\mathbb{E}[L(\widehat{\mathbf{w}}^{-n})]$ rather than
$\mathbb{E}[L(\widehat{\mathbf{w}})]$, but the difference between the
estimates is sufficiently small to be absorbed by $\epsilon$ since
$\frac{d}{n-d}\approx_{1+\epsilon}\frac{d}{n-d-1}$.

A key part of the analysis is showing that the leverage scores of the sketched
matrix are nearly uniform with high probability. This is established in the
following lemma.

###### Lemma 7

With probability $1-\delta$, for every $i\in[n]$ we have:

$\displaystyle(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{A})^{-1}\preceq
O(1)\cdot(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\qquad\text{and}\qquad\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}\,\approx_{1+\epsilon}\frac{d}{n-d},$

for $\epsilon=\tilde{O}(1/\sqrt{d})$. In particular, this implies that
$\ell_{i}(\mathbf{S}\mathbf{A})\approx_{1+2\epsilon}\frac{d}{n}$ for all $i$.

We now define a high probability event $\mathcal{E}$ that ensures uniformity
of the leverage scores:

$\displaystyle\mathcal{E}\ =\
\Big{[}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{A})^{-1}\preceq
O(1)\cdot(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1},\quad\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}\,\approx_{1+\epsilon}\frac{d}{n-d}\quad\forall_{i}\Big{]}.$

In particular, the event implies that
$1-\ell_{i}(\mathbf{S}\mathbf{A})\approx_{1+\epsilon^{\prime}}1-\frac{d}{m}$
for $\epsilon^{\prime}=\tilde{O}(\sqrt{d}/n)$, which allows us to simplify the
CV loss by approximating all the leverage scores with a uniform estimate. In
what follows, we will use $\mathbb{E}_{\mathcal{E}}$ to denote expectation
conditioned on the event $\mathcal{E}$:

$\displaystyle\mathbb{E}_{\mathcal{E}}[L_{{\mathrm{CV}}}]$
$\displaystyle\approx_{1+\epsilon^{\prime}}\mathbb{E}_{\mathcal{E}}\bigg{[}\sum_{i=1}^{n}\bigg{(}\frac{\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{A}\widehat{\mathbf{w}}-\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{b}}{1-d/n}\bigg{)}^{2}\bigg{]}=\Big{(}\frac{n}{n-d}\Big{)}^{2}\cdot\mathbb{E}_{\mathcal{E}}\big{[}\|\mathbf{S}(\mathbf{A}\widehat{\mathbf{w}}-\mathbf{b})\|^{2}\big{]}.$

Thus, it remains to approximate
$\mathbb{E}_{\mathcal{E}}\big{[}\min_{\mathbf{w}}\|\mathbf{S}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2}\big{]}$
in terms of
$L(\mathbf{w}^{*})=\min_{\mathbf{w}}\|\mathbf{A}\mathbf{w}-\mathbf{b}\|^{2}$.
Here, we are going to take advantage of the similar form of the two
quantities, in that they can be both written as squared distance between a
vector and a subspace. In the case of $L(\mathbf{w}^{*})$, it is the squared
distance between $\mathbf{b}$ and the column-span of $\mathbf{A}$, and it can
be written as:

$\displaystyle
L(\mathbf{w}^{*})=\mathbf{b}^{\scriptscriptstyle{\top}}(\mathbf{I}-\mathbf{H})\mathbf{b},\qquad\text{for}\quad\mathbf{H}=\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}\mathbf{A},$

where $\mathbf{H}$ is the so-called hat matrix, an orthogonal projection onto
the column-span of $\mathbf{A}$. In the sketched version, we get a squared
distance between $\mathbf{S}\mathbf{b}$ and the column-span of
$\mathbf{S}\mathbf{A}$, which can be written similarly, using a sketched
version of the hat matrix:

$\displaystyle\min_{\mathbf{w}}\|\mathbf{S}(\mathbf{A}\mathbf{w}-\mathbf{b})\|^{2}=\mathbf{b}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}\mathbf{b},\qquad\text{for}\quad\widehat{\mathbf{H}}=\mathbf{S}\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{A})^{-1}\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}.$

Thus, both quantities are quadratic forms applied to the vector $\mathbf{b}$.
We relate these two quadratic forms in expectation through the following
result.

###### Lemma 8

The sketched hat matrix $\widehat{\mathbf{H}}$ satisfies the following
approximate expectation formula:

$\displaystyle\mathbb{E}_{\mathcal{E}}\big{[}\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}\big{]}\approx_{1+\epsilon^{\prime}}\Big{(}1-\frac{d}{n}\Big{)}\cdot(\mathbf{I}-\mathbf{H}),\qquad\text{for}\quad\epsilon^{\prime}=\tilde{O}(\sqrt{d}/n).$

###### Remark 4

Here, $\approx_{1+\epsilon}$ refers to upper/lower bounds in terms of the
positive semidefinite ordering. This result extends an exact expectation
formula that holds when $\mathbf{S}$ is a Gaussian embedding.

Putting everything together, we obtain the following approximation of the
expected CV loss conditioned on $\mathcal{E}$:

$\displaystyle\mathbb{E}_{\mathcal{E}}[L_{\mathrm{CV}}]$
$\displaystyle\approx_{1+\epsilon^{\prime}}\Big{(}\frac{n}{n-d}\Big{)}^{2}\mathbf{b}^{\scriptscriptstyle{\top}}\mathbb{E}_{\mathcal{E}}\big{[}\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}\big{]}\mathbf{b}$
$\displaystyle\approx_{1+\epsilon^{\prime}}\Big{(}\frac{n}{n-d}\Big{)}^{2}\Big{(}1-\frac{d}{n}\Big{)}\cdot
L(\mathbf{w}^{*})=\frac{n}{n-d}\cdot L(\mathbf{w}^{*}).$

From this it follows that

$\mathbb{E}_{\mathcal{E}}[L_{{\mathrm{CV}}}]-L(\mathbf{w}^{*})\approx_{1+\epsilon}\frac{d}{n-d}\cdot
L(\mathbf{w}^{*})\approx_{1+\epsilon}\frac{d}{n-d-1}\cdot L(\mathbf{w}^{*}),$

for $\epsilon=\tilde{O}(1/\sqrt{d})$. Note that for this we have to switch
from $\epsilon^{\prime}=\tilde{O}(\sqrt{d}/n)$ to
$\epsilon=\tilde{O}(1/\sqrt{d})$, because the right-hand side gets scaled by
$d/n$. Finally, we convert from conditional expectation as follows:

$\displaystyle(1-\delta)\mathbb{E}_{\mathcal{E}}[L_{{\mathrm{CV}}}]\leq\mathbb{E}[L_{{\mathrm{CV}}}]$
$\displaystyle\leq\mathbb{E}_{\mathcal{E}}[L_{{\mathrm{CV}}}]+\delta\sup_{\mathbf{S}}L_{{\mathrm{CV}}}$
$\displaystyle\leq\mathbb{E}_{\mathcal{E}}[L_{{\mathrm{CV}}}]+\delta\cdot
n\cdot\sup_{\mathbf{S}}\|\mathbf{s}_{n}\|^{2}\|\mathbf{A}\widehat{\mathbf{w}}^{-n}-\mathbf{b}\|^{2}$
$\displaystyle\leq\mathbb{E}_{\mathcal{E}}[L_{{\mathrm{CV}}}]+\delta\cdot{\mathrm{poly}}(N,D),$

where in the last step, to avoid the corner cases where sketch
$\mathbf{S}\mathbf{A}$ is extremely ill-conditioned, we used that
$\widehat{\mathbf{w}}$, as well as $\mathbf{A}$ and $\mathbf{b}$, have entries
in $[-D,D]$, and that, by our definition, a LESS embedding matrix $\mathbf{S}$
has entries uniformly bounded by $O(\sqrt{Nk})$. Choosing sufficiently small
$\delta\leq\epsilon L(\mathbf{w}^{*})/{\mathrm{poly}}(N,D)$, we can obsorb
that bound into a logarithmic factor that appears in the upper bound on
$\epsilon$ (if we use $d$ non-zeros per row in LESS). Alternatively, we can
increase the number of non-zeros per row in LESS to $d\log(nd/\delta)$, in
which case, we avoid the dependence of $\epsilon$ on $D$. We note that another
simple way to impose a strict bound on the worst-case loss is to use a small
$\ell_{2}$ regularizer when computing the sketched least squares, ensuring
that the sketched problem is never ill-conditioned. Extending our analysis to
the regularized least squares setting is an interesting direction for future
work.

###  D.1 Proof of Lemma 8

Let
$\mathbf{U}=\mathbf{A}(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1/2}$
and note that $\mathbf{U}$ is an $N\times d$ matrix that satisfies
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{I}$. We can rewrite
both $\mathbf{H}$ and $\widehat{\mathbf{H}}$ by replacing $\mathbf{A}$ with
$\mathbf{U}$, namely:
$\mathbf{H}=\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}$ and
$\widehat{\mathbf{H}}=\mathbf{S}\mathbf{U}(\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}$.
To show the claim, it suffices to show that for every vector
$\mathbf{v}\in\mathbb{R}^{N}$, the following approximation holds:

$\displaystyle\mathbb{E}_{\mathcal{E}}[\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}\mathbf{v}]\,\approx_{1+\epsilon^{\prime}}\Big{(}1-\frac{d}{n}\Big{)}\,\mathbf{v}^{\scriptscriptstyle{\top}}(\mathbf{I}-\mathbf{H})\mathbf{v}.$
(8)

First, if $\mathbf{v}$ lies in the column-span of $\mathbf{U}$, i.e.,
$\mathbf{v}=\mathbf{U}\mathbf{x}$ for some $\mathbf{x}$, then clearly the
right-hand side is zero, and also:

$\displaystyle\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}\mathbf{v}=\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}\mathbf{x}-\mathbf{x}^{\scriptscriptstyle{\top}}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\widehat{\mathbf{H}}\mathbf{S}\mathbf{U}\mathbf{x}=0=\Big{(}1-\frac{d}{n}\Big{)}\,\mathbf{v}^{\scriptscriptstyle{\top}}(\mathbf{I}-\mathbf{H})\mathbf{v}.$

Thus, it suffices to show (8) for $\mathbf{v}$ orthogonal to the column-span
of $\mathbf{U}$, i.e., one that satisfies
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{v}=\mathbf{0}$. In that case,

$\displaystyle\Big{|}\mathbf{v}^{\scriptscriptstyle{\top}}\mathbb{E}_{\mathcal{E}}[\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}]\mathbf{v}-\Big{(}1-\frac{d}{n}\Big{)}\mathbf{v}^{\scriptscriptstyle{\top}}(\mathbf{I}-\mathbf{H})\mathbf{v}\Big{|}=\Big{|}\mathbf{v}^{\scriptscriptstyle{\top}}\mathbb{E}_{\mathcal{E}}[\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}]\mathbf{v}-\mathbf{v}^{\scriptscriptstyle{\top}}\mathbb{E}_{\mathcal{E}}[\mathbf{S}^{\scriptscriptstyle{\top}}\widehat{\mathbf{H}}\mathbf{S}]\mathbf{v}-\Big{(}1-\frac{d}{n}\Big{)}\Big{|}.$

We now focus on analyzing
$\mathbb{E}_{\mathcal{E}}[\mathbf{S}^{\scriptscriptstyle{\top}}\widehat{\mathbf{H}}\mathbf{S}]$.
To that end, let us use the shorthands
$\mathbf{Q}=(\gamma\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U})^{-1}$
and
$\mathbf{Q}_{-i}=(\gamma\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{U})^{-1}$,
where $\gamma=\frac{n}{n-d}$, so that
$\widehat{\mathbf{H}}=\mathbf{S}\mathbf{U}\gamma\mathbf{Q}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}$.
We also use the formula
$\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}=\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}/\gamma_{i}$,
where
$\gamma_{i}=1+\gamma\mathbf{s}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{s}_{i}$,
which is a consequence of the Sherman-Morrison formula (Lemma 4). Recall that,
from Lemma 7, the event $\mathcal{E}$ implies that $\gamma_{i}\approx\gamma$.
Also, let us use a simplifying shorthand of
${\tilde{\mathbf{s}}}_{i}=\sqrt{n}\,\mathbf{s}_{i}$, so that
$\mathbb{E}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]=\mathbf{I}$.
We can now rewrite the expectation as follows (where $i$ is any fixed index):

$\displaystyle\mathbb{E}_{\mathcal{E}}[\mathbf{S}^{\scriptscriptstyle{\top}}\widehat{\mathbf{H}}\mathbf{S}]$
$\displaystyle=\mathbb{E}_{\mathcal{E}}[\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\mathbf{U}\gamma\mathbf{Q}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}]=\mathbb{E}_{\mathcal{E}}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\gamma\mathbf{Q}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}]=\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma}{\gamma_{i}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}^{\scriptscriptstyle{\top}}\mathbf{S}\Big{]}$
$\displaystyle=\frac{1}{n}\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma}{\gamma_{i}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]+\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma}{\gamma_{i}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\Big{]}$
$\displaystyle=\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma_{i}-1}{\gamma_{i}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]+\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma}{\gamma_{i}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\Big{]}$
$\displaystyle=\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma-1}{\gamma}\cdot{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\Big{]}+\mathbb{E}_{\mathcal{E}}\big{[}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\big{]}$
$\displaystyle\quad+\mathbb{E}_{\mathcal{E}}\Big{[}\Big{(}\frac{1}{\gamma}-\frac{1}{\gamma_{i}}\Big{)}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\Big{]}+\mathbb{E}_{\mathcal{E}}\Big{[}\Big{(}\frac{\gamma}{\gamma_{i}}-1\Big{)}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\Big{]}=:\mathbf{T}_{1}+\mathbf{T}_{2}+\mathbf{T}_{3}+\mathbf{T}_{4}.$

Thus, it follows that:

$\displaystyle\mathbb{E}_{\mathcal{E}}\big{[}\mathbf{S}^{\scriptscriptstyle{\top}}(\mathbf{I}-\widehat{\mathbf{H}})\mathbf{S}\big{]}$
$\displaystyle=\mathbb{E}_{\mathcal{E}}\big{[}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]-\mathbb{E}_{\mathcal{E}}\Big{[}\frac{\gamma-1}{\gamma}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\Big{]}-(\mathbf{T}_{2}+\mathbf{T}_{3}+\mathbf{T}_{4})$
$\displaystyle=\Big{(}1-\frac{d}{n}\Big{)}\mathbb{E}_{\mathcal{E}}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]-(\mathbf{T}_{2}+\mathbf{T}_{3}+\mathbf{T}_{4}).$

In the following steps, we will show that
$\mathbb{E}_{\mathcal{E}}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]\approx\mathbf{I}$
and that $\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{T}_{k}\mathbf{v}\approx
0$ for $k\in\\{2,3,4\\}$. First, denoting
$\delta=\mathrm{Pr}(\neg\mathcal{E})$, we have:

$\displaystyle\|\mathbb{E}_{\mathcal{E}}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]-\mathbf{I}\|=\frac{1}{1-\delta}\|\mathbb{E}[\mathbf{1}_{\neg\mathcal{E}}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}]\|\leq\frac{\delta}{1-\delta}\cdot
O(Nk),$

where we used that $\|{\tilde{\mathbf{s}}}_{i}\|^{2}\leq O(Nk)$. Next, observe
that for $\mathbf{v}$ orthogonal to the column-span of $\mathbf{U}$, we have
$\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{T}_{2}\mathbf{v}=\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v}=0$
because $\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{U}=\mathbf{0}$.
Furthermore, using Lemma 7, conditioned on $\mathcal{E}$ we know that
$|\gamma_{i}-\gamma|=\tilde{O}(\sqrt{d}/n)$, which allows us to control the
remaining terms. In the following, we also use the fact that
$\mathbb{E}_{\mathcal{E}}[X]\leq\frac{1}{1-\delta}\mathbb{E}[X]$ for any non-
negative $X$:

$\displaystyle\big{|}\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{T}_{3}\mathbf{v}\big{|}$
$\displaystyle\leq\mathbb{E}_{\mathcal{E}}\Big{[}\frac{|\gamma_{i}-\gamma|}{\gamma\gamma_{i}}(\mathbf{v}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i})^{2}\Big{]}$
$\displaystyle\leq\tilde{O}(\sqrt{d}/n)\cdot\mathbb{E}_{\mathcal{E}}[(\mathbf{v}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i})^{2}]$
$\displaystyle\leq\tilde{O}(\sqrt{d}/n)\cdot\mathbf{v}^{\scriptscriptstyle{\top}}\mathbb{E}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}]\mathbf{v}=\tilde{O}(\sqrt{d}/n)\cdot\|\mathbf{v}\|^{2}.$

Similarly, we bound the last term, using the Cauchy-Schwartz inequality:

$\displaystyle\big{|}\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{T}_{4}\mathbf{v}\big{|}$
$\displaystyle\leq\tilde{O}(\sqrt{d}/n)\cdot\mathbb{E}_{\mathcal{E}}\big{[}|\mathbf{v}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i}|\cdot|{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v}|\big{]}$
$\displaystyle\leq\tilde{O}(\sqrt{d}/n)\cdot\sqrt{\mathbb{E}_{\mathcal{E}}\big{[}(\mathbf{v}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i})^{2}\big{]}}\cdot\sqrt{\mathbb{E}_{\mathcal{E}}\big{[}({\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v})^{2}\big{]}}.$

As before, we have
$\mathbb{E}_{\mathcal{E}}[(\mathbf{v}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i})^{2}]\leq
O(1)\cdot\|\mathbf{v}\|^{2}$. For the last term in the product, we use a LESS
embedding property of the random vector vector
$\mathbf{U}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{i}$ that was shown
in [17, Lemma 28], which they called the Bai-Silverstein property. In
particular, this property implies that for any vector $\mathbf{v}$, we have
$\mathbb{E}[({\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{v})^{2}]\leq
O(1)\cdot\|\mathbf{v}\|^{2}$ (this property also holds for any sub-gaussian
vector, so it could be easily inferred from our Theorem 1, after conditioning
on a high-probability event). We also use that the event $\mathcal{E}$ implies
that $\|\mathbf{Q}_{-i}\|=O(1)$. It follows that:

$\displaystyle\mathbb{E}_{\mathcal{E}}\big{[}({\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v})^{2}\big{]}$
$\displaystyle\leq\frac{1}{1-\delta}\mathbb{E}\big{[}({\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}\mathbf{U}\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v})^{2}\mid\|\mathbf{Q}_{-i}\|=O(1)\big{]}$
$\displaystyle\leq
O(1)\cdot\mathbb{E}\big{[}\|\mathbf{Q}_{-i}\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v}\|^{2}\mid\|\mathbf{Q}_{-i}\|=O(1)\big{]}\leq
O(1)\cdot\mathbb{E}\big{[}\|\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{v}\|^{2}\big{]}$
$\displaystyle=O(1)\cdot\mathbb{E}\bigg{[}\Big{\|}\frac{1}{n}\sum_{j\neq
i}\mathbf{U}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{j}{\tilde{\mathbf{s}}}_{j}^{\scriptscriptstyle{\top}}\mathbf{v}\Big{\|}^{2}\bigg{]}\leq
O(1)\cdot\frac{1}{n}\mathbb{E}\big{[}\|\mathbf{U}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{j}{\tilde{\mathbf{s}}}_{j}^{\scriptscriptstyle{\top}}\mathbf{v}\|^{2}\big{]},$

where in the last step we used that ${\tilde{\mathbf{s}}}_{j}$ are
independent, and $\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{v}=\mathbf{0}$.
To bound
$\mathbb{E}[\|\mathbf{U}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{j}{\tilde{\mathbf{s}}}_{j}^{\scriptscriptstyle{\top}}\mathbf{v}\|^{2}]$,
we will use the definition of a LESS embedding (Definition 1), where we let
$I_{i}$ be the $k$ random indices that indicate the non-zero entries of
${\tilde{\mathbf{s}}}_{j}=\sum_{i=1}^{k}\frac{r_{i}\mathbf{e}_{I_{i}}}{\sqrt{kp_{I_{i}}}}$
and $r_{i}$’s are the random sign variables.

$\displaystyle\mathbb{E}\big{[}\|\mathbf{U}^{\scriptscriptstyle{\top}}{\tilde{\mathbf{s}}}_{j}{\tilde{\mathbf{s}}}_{j}^{\scriptscriptstyle{\top}}\mathbf{v}\|^{2}\big{]}\leq
k\,\mathbb{E}\bigg{[}\frac{\|\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{e}_{I_{1}}\|^{2}(\mathbf{e}_{I_{1}}^{\scriptscriptstyle{\top}}\mathbf{v})^{2}}{k^{2}p_{I_{1}}^{2}}\bigg{]}+k^{2}\,\mathbb{E}\bigg{[}\frac{\|\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{e}_{I_{1}}\|^{2}}{kp_{I_{1}}}\bigg{]}\mathbb{E}\bigg{[}\frac{(\mathbf{e}_{I_{2}}^{\scriptscriptstyle{\top}}\mathbf{v})^{2}}{kp_{I_{2}}}\bigg{]}\leq\big{(}Cd/k+Cd\big{)}\cdot\|\mathbf{v}\|^{2},$

where, to eliminate the cross-terms when expanding the square, we used
independence of $I_{i}$’s and that
$\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{v}=\mathbf{0}$. In the second
step, we used that
$p_{I_{i}}\geq\|\mathbf{U}^{\scriptscriptstyle{\top}}\mathbf{e}_{I_{i}}\|^{2}/(Cd)$.
Thus,
$|\mathbf{v}^{\scriptscriptstyle{\top}}\mathbf{T}_{4}\mathbf{v}|\leq\tilde{O}(d/n^{1.5})\|\mathbf{v}\|^{2}=\tilde{O}(\sqrt{d}/n)\|\mathbf{v}\|^{2}$.
We conclude that for sufficiently small $\delta$, claim (8) is satisfied with
$\epsilon=\tilde{O}(\sqrt{d}/n)$, which completes the proof.

###  D.2 Proof of Lemma 7

We note that the first claim, i.e, that
$(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{S}_{-i}\mathbf{S}_{-i}\mathbf{A})^{-1}\preceq
O(1)\cdot(\mathbf{A}^{\scriptscriptstyle{\top}}\mathbf{A})^{-1}$ follows
easily for LESS embeddings with $d$ non-zeros per row of size $n\geq
O(d\log(d/\delta))$, by relying on the subspace embedding property, shown in
[17, Lemma 12], which essentially follows from the matrix Bernstein inequality
(Lemma 3). An alternative way of showing this claim is to rely on the subspace
embedding guarantee for sub-gaussian embeddings, and then convert that to a
result for LESS embeddings via our Corollary 1. This approach is given in
Corollary 4. Interestingly, this leads to a slightly different trade-off in
terms of the logarithmic factors. Namely, if we consider a LESS embedding with
$d\log(nd/\delta)$ non-zeros per row, slightly more than recommended by [17],
then the subspace embedding guarantee holds for $n\geq O(d)$ samples, rather
than $n\geq O(d\log d)$. Such a guarantee is new for LESS embeddings, since it
cannot be obtained via the matrix Bernstein inequality.

We now focus on establishing the second part of the claim, i.e., the
uniformity of the leverage scores, which is were the Hanson-Wright inequality
is essential. Since we want to show a statement that holds with high
probability, we can use Corollary 1 to show this claim for an $n\times d$ sub-
gaussian design $\mathbf{X}=\mathbf{Z}\mathbf{\Sigma}$ with $\mathbf{Z}$
consisting of mean zero isotropic rows that satisfy the Hanson-Wright
inequality (2) with constant $\tilde{O}(1)$, and then convert the claim to
LESS embeddings. In this case, our goal is to show that for
$n\geq\tilde{O}(d)$, with probability $1-\delta$:

$\displaystyle\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\big{(}\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i}\big{)}^{-1}\mathbf{x}_{i}\approx_{1+\epsilon}\frac{d}{n-d}.$

Note that, by Hanson-Wright, with high probability we immediately have that
$\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\big{(}\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i}\big{)}^{-1}\mathbf{x}_{i}\approx\mathrm{tr}\,\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}$
(this will be formalized later on), so the main challenge is in showing that
the trace functional is also concentrated. Below, we provide a slightly more
general result, regarding the concentration properties of a certain class of
trace functionals for a sub-gaussian design, which may be of independent
interest. (In the proof we use $a\lesssim b$ to denote $a\leq cb$ for an
absolute constant $c$).

###### Lemma 9

Suppose that $\mathbf{X}$ consists of rows
$\mathbf{x}_{1}^{\scriptscriptstyle{\top}},\mathbf{x}_{2}^{\scriptscriptstyle{\top}},...,\mathbf{x}_{n}^{\scriptscriptstyle{\top}}$
that are i.i.d. sampled from a $d$-variate distribution with covariance
$\mathbb{E}[\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}]=\mathbf{\Sigma}$
such that $\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}$ satisfies the Hanson-Wright
inequality (2) with constant $K=\tilde{O}(1)$ and $n\geq\tilde{O}(d)$.
Consider
$F(\mathbf{X})=\mathrm{tr}\,\mathbf{B}(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}+\mathbf{C})^{-1}$
for some psd matrices $\mathbf{B}$ and $\mathbf{C}$. Then, there is a scalar
$\tilde{F}$ such that with probability $1-\delta$:

$\displaystyle|F(\mathbf{X})-\tilde{F}|\leq
C\|\mathbf{M}\|\cdot\frac{\mathrm{r}(\mathbf{M})+K^{4}}{\sqrt{n}}\log
1/\delta,$

where
$\mathbf{M}=\mathbf{B}^{1/2}(\mathbf{\Sigma}+\mathbf{C})^{-1}\mathbf{B}^{1/2}$,
$\mathrm{r}(\mathbf{M})=\mathrm{tr}(\mathbf{M})/\|\mathbf{M}\|$, and $C$ is an
absolute constant.

Proof  Note that $\mathbf{C}$ may not be positive definite (in fact, of our
primary interest is when $\mathbf{C}=\mathbf{0}$), so function $F(\mathbf{X})$
may not be well defined everywhere because of the inverse. In particular, it
may not have a bounded expectation. So, to define our scalar $\tilde{F}$, we
first construct a high-probability event $\mathcal{E}$ which will ensure that
$F(\mathbf{X})$ is not only well-defined, but also sufficiently well behaved.
Without loss of generality, assume that $n$ is even. Now, define:

$\displaystyle\mathcal{E}_{1}=\Big{[}\frac{2}{n}\sum_{i=1}^{n/2}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\succeq\frac{1}{C}\mathbf{\Sigma}\Big{]},\qquad\mathcal{E}_{2}=\Big{[}\frac{2}{n}\sum_{i=n/2+1}^{n}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\succeq\frac{1}{C}\mathbf{\Sigma}\Big{]},$

with $\mathcal{E}=\mathcal{E}_{1}\wedge\mathcal{E}_{2}$. Since
$n\geq\tilde{O}(d)$, by relying on the subspace embedding property (as
discussed at the beginning of the section), each of the events
$\mathcal{E}_{1}$ and $\mathcal{E}_{2}$ holds with probability at least
$1-\delta$. We will use $\tilde{F}=\mathbb{E}_{\mathcal{E}}[F(\mathbf{X})]$ as
our scalar. Note that, conditioned on $\mathcal{E}$, the function $F$
satisfies $0\leq F(\mathbf{X})\leq C\,\mathrm{tr}\,\mathbf{M}$, so in
particular the conditional expectation is well-defined.

Next, we proceed with the following decomposition of the trace (following
standard literature in random matrix theory). Let
$\mathbf{Q}=(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X}+\mathbf{C})^{-1}$
and
$\mathbf{Q}_{-i}=(\frac{1}{n}\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i}+\mathbf{C})^{-1}$
where $\mathbf{X}_{-i}$ denotes $\mathbf{X}$ without the $i$ row. Also, let
$\mathbb{E}_{i}$ denote expectation conditioned on
$\mathbf{x}_{1},...,\mathbf{x}_{i}$. We have, conditioned on $\mathcal{E}$
that:

$\displaystyle F(\mathbf{X})-\tilde{F}$
$\displaystyle=\mathbb{E}_{n}[\mathrm{tr}\mathbf{B}\mathbf{Q}\mid\mathcal{E}]-\mathbb{E}_{0}[\mathrm{tr}\mathbf{B}\mathbf{Q}\mid\mathcal{E}]$
$\displaystyle=\sum_{i=1}^{n}(\mathbb{E}_{i}-\mathbb{E}_{i-1})[\mathrm{tr}\mathbf{B}\mathbf{Q}\mid\mathcal{E}]$
$\displaystyle=\sum_{i=1}^{n}(\mathbb{E}_{i}-\mathbb{E}_{i-1})[\mathrm{tr}\mathbf{B}(\mathbf{Q}-\mathbf{Q}_{-i})\mid\mathcal{E}]+(\mathbb{E}_{i}-\mathbb{E}_{i-1})[\mathrm{tr}\mathbf{B}\mathbf{Q}_{-i}\mid\mathcal{E}].$

We proceed to bound each of the two terms. The latter can be bounded
straightforwardly. Note that one of $\mathcal{E}_{1}$ or $\mathcal{E}_{2}$ is
independent of $\mathbf{x}_{i}$. Without loss of generality suppose that it is
$\mathcal{E}_{1}$. Then, we have:

$\displaystyle\big{|}(\mathbb{E}_{i}-\mathbb{E}_{i-1})[\mathrm{tr}\mathbf{B}\mathbf{Q}_{-i}\mid\mathcal{E}]\big{|}\leq\frac{\delta}{1-\delta}\cdot\big{|}(\mathbb{E}_{i}-\mathbb{E}_{i-1})[\mathrm{tr}\mathbf{B}\mathbf{Q}_{-i}\mid\mathcal{E}_{1}\wedge\neg\mathcal{E}_{2}]\big{|}\lesssim\delta\cdot\mathrm{tr}\,\mathbf{M}.$

To bound the first term, we use the Sherman-Morrison rank-one update formula
(Lemma 4), observing that:

$\displaystyle|\mathrm{tr}\mathbf{B}(\mathbf{Q}-\mathbf{Q}_{-i})|=\frac{\frac{1}{n}\mathrm{tr}\,\mathbf{B}\mathbf{Q}_{-i}\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{Q}_{-i}}{1+\frac{1}{n}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{Q}_{-i}\mathbf{x}_{i}}\leq\frac{1}{n}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}\mathbf{Q}_{-i}\mathbf{B}\mathbf{Q}_{-i}\mathbf{x}_{i}=\frac{1}{n}\mathbf{z}_{i}^{\scriptscriptstyle{\top}}\tilde{\mathbf{M}}_{i}\mathbf{z}_{i},$

where we let $\mathbf{z}_{i}=\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}$ and
$\tilde{\mathbf{M}}_{i}=\mathbf{\Sigma}^{1/2}\mathbf{Q}_{-i}\mathbf{B}\mathbf{Q}_{-i}\mathbf{\Sigma}^{1/2}$.
Now, integrating out the Hanson-Wright inequality for $\mathbf{z}_{i}$ into a
moment bound, letting
$Z_{i}=\sqrt{\mathbf{z}_{i}^{\scriptscriptstyle{\top}}\tilde{\mathbf{M}}_{i}\mathbf{z}_{i}}$
and $M_{i}=\mathbb{E}[Z_{i}\mid\mathcal{E},\tilde{\mathbf{M}}_{i}]$, we have:

$\displaystyle\mathbb{E}\big{[}|\mathbf{z}_{i}^{\scriptscriptstyle{\top}}\tilde{\mathbf{M}}_{i}\mathbf{z}_{i}|^{p}\mid\mathcal{E}\big{]}^{1/p}$
$\displaystyle=\mathbb{E}\big{[}(M_{i}-M_{i}+Z_{i})^{2p}\mid\mathcal{E}]^{1/p}$
$\displaystyle\lesssim\mathbb{E}\big{[}\mathbb{E}[M_{i}^{2p}\mid\tilde{\mathbf{M}}_{i}]\mid\mathcal{E}\big{]}^{1/p}+\mathbb{E}\big{[}\mathbb{E}[|Z_{i}-M_{i}|^{2p}\mid\tilde{\mathbf{M}}_{i}]\mid\mathcal{E}\big{]}^{1/p}$
$\displaystyle\lesssim\mathbb{E}\big{[}(\mathrm{tr}\tilde{\mathbf{M}}_{i})^{p}\mid\mathcal{E}\big{]}^{1/p}+\mathbb{E}\big{[}(C\sqrt{p}K^{2}\sqrt{\|\mathbf{M}\|})^{2p}\mid\mathcal{E}\big{]}^{1/p}$
$\displaystyle\lesssim\mathrm{tr}\mathbf{M}+CpK^{4}\|\mathbf{M}\|.$

where the last step follows because conditioned on $\mathcal{E}$ we have
$\mathrm{tr}\,\tilde{\mathbf{M}}_{i}\lesssim\mathrm{tr}\,\mathbf{M}$ and
$\|\tilde{\mathbf{M}}_{i}\|\lesssim\|\mathbf{M}\|$. Using Burkholder’s
inequality (e.g., see [28]) for the martingale difference sequence
$X_{i}=(\mathbb{E}_{i}-\mathbb{E}_{i-1})[\mathrm{tr}\mathbf{B}\mathbf{Q}\mid\mathcal{E}]$,
we get:

$\displaystyle\mathbb{E}\big{[}|F(X)-\tilde{F}|^{p}\big{]}^{1/p}$
$\displaystyle=\mathbb{E}\bigg{[}\Big{(}\sum_{i=1}^{n}X_{i}\Big{)}^{p}\bigg{]}^{1/p}\leq
Cp\cdot\Bigg{(}\mathbb{E}\bigg{[}\Big{(}\sum_{i=1}^{n}\mathbb{E}_{i-1}[X_{i}^{2}]\Big{)}^{p/2}\bigg{]}+\sum_{i=1}^{n}\mathbb{E}\big{[}X_{i}^{p}\big{]}\Bigg{)}^{1/p}$
$\displaystyle\lesssim\frac{p}{n}\cdot\Bigg{(}\Big{(}n((\mathrm{tr}\,\mathbf{M})^{2}+K^{4}\|\mathbf{M}\|)\Big{)}^{p/2}+n(\mathrm{tr}\,\mathbf{M}+pK^{4}\|\mathbf{M}\|)^{p}\Bigg{)}^{1/p}$
$\displaystyle\lesssim
p\big{(}\mathrm{tr}\,\mathbf{M}+K^{4}\|\mathbf{M}\|\big{)}\cdot\Big{(}\frac{1}{\sqrt{n}}+\frac{p}{n^{1-1/p}}\Big{)}.$

Applying Markov’s inequality with $p=\log(1/\delta)$, we obtain the desired
result.

Note that in Lemma 9 we did not explicitly state what is the scalar
$\tilde{F}$ around which the trace functional $F(\mathbf{X})$ concentrates.
Normally, one would expect this quantity to be the expectation of
$F(\mathbf{X})$, however since this expectation may not exist, we instead
realy on the notion of a nearly-unbiased estimator. For the trace functional
we are interested in, i.e.,
$F(\mathbf{X})=\mathrm{tr}\,\mathbf{\Sigma}(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}$,
a nearly-unbiased estimator has been derived by [17], and we use that result
in the following corollary.

###### Corollary 5

Suppose that $\mathbf{X}$ consists of rows
$\mathbf{x}_{1}^{\scriptscriptstyle{\top}},\mathbf{x}_{2}^{\scriptscriptstyle{\top}},...,\mathbf{x}_{n}^{\scriptscriptstyle{\top}}$
that are i.i.d. sampled from a $d$-variate distribution with covariance
$\mathbb{E}[\mathbf{x}_{i}\mathbf{x}_{i}^{\scriptscriptstyle{\top}}]=\mathbf{\Sigma}$
such that $\mathbf{\Sigma}^{-1/2}\mathbf{x}_{i}$ satisfies the Hanson-Wright
inequality (2) with $K=\tilde{O}(1)$ and $n\geq\tilde{O}(d)$. Then, with
probability $1-\delta$:

$\displaystyle\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}\approx_{1+\epsilon}\frac{d}{n-d}\quad\text{where}\quad\epsilon\leq\tilde{O}\bigg{(}\frac{\log
1/\delta}{\sqrt{n}}\bigg{)}.$

Proof  Here, we use the result from [17] which shows that
$\mathbf{\Sigma}^{1/2}(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}\mathbf{\Sigma}^{1/2}$
is an $(\epsilon^{\prime},\delta)$-unbiased estimator of
$\frac{n}{n-d}\mathbf{I}$, where $\epsilon^{\prime}=O(K^{4}\sqrt{d}/n)$.
Namely, this means that there is an event $\mathcal{E}$ with probability
$1-\delta$ such that
$\mathbb{E}[\mathbf{\Sigma}^{1/2}(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}\mathbf{\Sigma}^{1/2}\mid\mathcal{E}]\approx_{1+\epsilon}\frac{n}{n-d}\mathbf{I}$
and conditioned on $\mathcal{E}$ we have
$\mathbf{\Sigma}^{1/2}(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}\mathbf{\Sigma}^{1/2}\preceq
O(1)\cdot\frac{n}{n-d}\mathbf{I}$. As a consequence, following the argument
analogous to Lemma 34 in [17], we can easily show that if
$F(\mathbf{X})=\mathrm{tr}\mathbf{\Sigma}(\frac{1}{n}\mathbf{X}^{\scriptscriptstyle{\top}}\mathbf{X})^{-1}$,
then our scalar $\tilde{F}$ used in Lemma 9 also satisfies
$\tilde{F}\approx_{1+\epsilon}\frac{n}{n-d}\cdot d$. So, combining this with
Lemma 9, and observing that in this case we have $\mathbf{M}=\mathbf{I}$
completes the proof.

To conclude the proof of Lemma 7, we combine the corollary with the Hanson-
Wright inequality, which states that, conditioned on
$\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i}$:

$\displaystyle\mathrm{Pr}$
$\displaystyle\Big{\\{}|\mathbf{x}_{i}^{\scriptscriptstyle{\top}}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\mathbf{x}_{i}-\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}|\geq\epsilon\cdot\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\Big{\\}}$
$\displaystyle\leq
2\exp\Big{(}-c\min\Big{\\{}\frac{\epsilon^{2}\cdot(\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1})^{2}}{K^{4}\|\mathbf{\Sigma}^{1/2}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\mathbf{\Sigma}^{1/2}\|_{F}^{2}},\frac{\epsilon\cdot\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}}{K^{2}\|\mathbf{\Sigma}^{1/2}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\mathbf{\Sigma}^{1/2}\|}\Big{\\}}\Big{)}$
$\displaystyle\leq\exp\Big{(}-c\epsilon^{2}\frac{\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}}{K^{4}\|\mathbf{\Sigma}^{1/2}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\mathbf{\Sigma}^{1/2}\|}\Big{)}.$

Using that $K=\tilde{O}(1)$, as well as that with high probability,
$\mathrm{tr}\mathbf{\Sigma}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\approx\frac{d}{n-d}$
and
$\|\mathbf{\Sigma}^{1/2}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\mathbf{\Sigma}^{1/2}\|=O(1/n)$,
we conclude that with high probability, we also have:

$\displaystyle\mathbf{x}_{i}^{\scriptscriptstyle{\top}}(\mathbf{X}_{-i}^{\scriptscriptstyle{\top}}\mathbf{X}_{-i})^{-1}\mathbf{x}_{i}\approx_{1+\epsilon}\frac{d}{n-d}\qquad\text{for}\qquad\epsilon=\tilde{O}(1/\sqrt{d}).$

By coupling a LESS embedding with the sub-gaussian design, we obtain the
desired claim.

##  Appendix E Constrained least squares: Proof of Theorem 3

The claim can be shown by using our reduction from a LESS embedding to a sub-
gaussian embedding. Namely, let $\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}$
be the sub-gaussian design that is coupled with the LESS embedding
$\sqrt{n}\mathbf{S}\mathbf{A}$ in Corollary 1. Recall that we have
$\widetilde{\mathbf{\Sigma}}\approx_{1+\delta}\mathbf{\Sigma}$. To simplify
the reduction, let us also define the following sketching matrix:

$\displaystyle\widetilde{\mathbf{S}}=\frac{1}{\sqrt{n}}\big{(}\mathbf{Z}\widetilde{\mathbf{\Sigma}}^{1/2}\mathbf{\Sigma}^{-1/2}\mathbf{U}^{\scriptscriptstyle{\top}}+\mathbf{G}(\mathbf{I}-\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}})\big{)},$

where $\mathbf{U}=\mathbf{A}\mathbf{\Sigma}^{-1/2}$ and $\mathbf{G}$ is an
$n\times N$ Gaussian matrix. It is easy to verify that
$\sqrt{n}\widetilde{\mathbf{S}}\mathbf{A}=\mathbf{Z}\mathbf{\Sigma}^{1/2}=\sqrt{n}\mathbf{S}\mathbf{A}$
with probability $1-\delta$, since
$\mathbf{I}-\mathbf{U}\mathbf{U}^{\scriptscriptstyle{\top}}$ is the projection
onto the complement of the column-span of $\mathbf{A}$. Moreover, we have
$\mathbb{E}[\widetilde{\mathbf{S}}^{\scriptscriptstyle{\top}}\widetilde{\mathbf{S}}]\approx_{1+\delta}\mathbf{I}$,
with each row of $\sqrt{n}\widetilde{\mathbf{S}}$ being mean zero and
$\tilde{O}(1)$-sub-gaussian. Thus, it suffices to verify that the analysis of
[40] still works when the sketching matrix is only approximately isotropic.
Note that the difference here is very small since we can easily let
$\delta\ll\epsilon$. It suffices to show that their Proposition 1 can be
adapted. This proposition is itself a corollary of Theorem D from [36]. We can
easily obtain a variant that suits our setup.

###### Proposition 2

Let ${\tilde{\mathbf{s}}}_{1},...,{\tilde{\mathbf{s}}}_{n}$ be i.i.d.
$N$-dimensional samples from a zero mean $K$-sub-gaussian distribution with
covariance
$\mathbb{E}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]\approx_{1+\delta}\mathbf{I}$.
For any subset $\mathcal{Y}\subseteq\mathcal{S}^{N-1}$, if $n\geq
CK^{2}\mathbb{W}^{2}(\mathcal{Y})/\epsilon^{2}$ and $\delta\leq c\epsilon$,
then:

$\displaystyle\sup_{y\in\mathcal{Y}}\Big{|}\mathbf{y}^{\scriptscriptstyle{\top}}\Big{(}\frac{1}{n}\sum_{i=1}^{n}{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}-\mathbf{I}\Big{)}\mathbf{y}\Big{|}\leq\epsilon,$

with probability at least $1-2\exp(-c\epsilon^{2}n/K^{4})$.

The proposition follows from Theorem D of [36], with functions
$f_{y}(\mathbf{s})=\frac{\mathbf{s}^{\scriptscriptstyle{\top}}\mathbf{y}}{\|\mathbf{y}\|_{\mathbf{M}}}$,
where
$\mathbf{M}=\mathbb{E}[{\tilde{\mathbf{s}}}_{i}{\tilde{\mathbf{s}}}_{i}^{\scriptscriptstyle{\top}}]$.
Remainder of the analysis in [40] stays the same.

##  Appendix F Lower bound: Proof of Theorem 4

First, observe that since $\mathbf{z}$ is isotropic, it has to satisfy
$\|\mathbf{z}\|_{\psi_{2}}=\Omega(1)=\Omega(\|\mathbf{r}\|_{\psi_{2}})$, so
without loss of generality, assume that
$k\leq\|\mathbf{x}\|_{\psi_{2}}^{2}\leq M^{2}$. Note that for any $\mathbf{z}$
specified in the theorem, there is a random vector $\tilde{\mathbf{z}}$ that
can be coupled with $\tilde{\mathbf{x}}$ so that conditioned on an event
$\mathcal{E}$ we have $\tilde{\mathbf{z}}=\tilde{\mathbf{x}}$ and: