Algorithm

# Robust and Sample Optimal Algorithms for
PSD Low-Rank Approximation

Ainesh Bakshi
CMU
<EMAIL_ADDRESS>Ainesh Bakshi and David Woodruff acknowledge support in
part from NSF No. CCF-1815840. Nadiia Chepurko
MIT
<EMAIL_ADDRESS>David P. Woodruff11footnotemark: 1
CMU
<EMAIL_ADDRESS>
###### Abstract

Recently, Musco and Woodruff (FOCS, 2017) showed that given an $n\times n$
positive semidefinite (PSD) matrix $\mathbf{A}$, it is possible to compute a
$(1+\epsilon)$-approximate relative-error low-rank approximation to
$\mathbf{A}$ by querying $\widetilde{O}(nk/\epsilon^{2.5})$ entries of
$\mathbf{A}$ in time
$\widetilde{O}(nk/\epsilon^{2.5}+nk^{\omega-1}/\epsilon^{2(\omega-1)})$. They
also showed that any relative-error low-rank approximation algorithm must
query ${\Omega}(nk/\epsilon)$ entries of $\mathbf{A}$, this gap has since
remained open. Our main result is to resolve this question by obtaining an
optimal algorithm that queries $\widetilde{O}(nk/\epsilon)$ entries of
$\mathbf{A}$ and outputs a relative-error low-rank approximation in
$\widetilde{O}(n\cdot(k/\epsilon)^{\omega-1})$ time. Note, our running time
improves that of Musco and Woodruff, and matches the information-theoretic
lower bound if the matrix-multiplication exponent $\omega$ is $2$.

We then extend our techniques to negative-type distance matrices. Here, our
input is a pair-wise distance matrix $\mathbf{A}$ corresponding to a point set
$\mathcal{P}=\\{x_{1},x_{2},\ldots,x_{n}\\}$ such that
$\mathbf{A}_{i,j}=\|x_{i}-x_{j}\|^{2}_{2}$. Bakshi and Woodruff (NeurIPS,
2018) showed a bi-criteria, relative-error low-rank approximation for
negative-type metrics. Their algorithm queries
$\widetilde{O}(nk/\epsilon^{2.5})$ entries and outputs a rank-$(k+4)$ matrix.
We show that the bi-criteria guarantee is not necessary and obtain an
$\widetilde{O}(nk/\epsilon)$ query algorithm, which is optimal. Our algorithm
applies to all distance matrices that arise from metrics satisfying negative-
type inequalities, including $\ell_{1},\ell_{2},$ spherical metrics,
hypermetrics and effective resistances on a graph. We also obtain faster
algorithms for ridge regression.

Next, we introduce a new robust low-rank approximation model which captures
PSD matrices that have been corrupted with noise. We assume that the Frobenius
norm of the corruption is bounded. Here, we relax the notion of approximation
to additive-error, since it is information-theoretically impossible to obtain
a relative-error approximation in this setting. While a sample complexity
lower bound precludes sublinear algorithms for arbitrary PSD matrices, we
provide the first sublinear time and query algorithms when the corruption on
the diagonal entries is bounded. As a special case, we show sample-optimal
sublinear time algorithms for low-rank approximation of correlation matrices
corrupted by noise.

###### Contents

  1. 1 Introduction
    1. 1.1 Our Results
  2. 2 Technical Overview
    1. 2.1 Sample-Optimal Low-Rank Approximation
    2. 2.2 Robust Low-Rank Approximation
  3. 3 Preliminaries and Notation
  4. 4 Relative Error PSD Low-Rank Approximation
    1. 4.1 Structured Projections to Low-Rank Approximation
    2. 4.2 Spectral Regression
    3. 4.3 Sample-Optimal Algorithm
    4. 4.4 Negative-Type Distances
    5. 4.5 Ridge Regression
  5. 5 Robust Low-Rank Approximation
    1. 5.1 Lower Bound for Robust PSD Low-Rank Approximation
    2. 5.2 Robust Sublinear Low-Rank Approximation Algorithms

##  1 Introduction

Low-rank approximation is one of the most common dimensionality reduction
techniques, whereby one replaces a large matrix $\mathbf{A}$ with a low-rank
factorization $\mathbf{U}\cdot\mathbf{V}\approx\mathbf{A}$. Such a
factorization provides a compact way of storing $\mathbf{A}$ and allows one to
multiply $\mathbf{A}$ quickly by a vector. It is used as an algorithmic
primitive in clustering [DFK+04, McS01], recommendation systems [DKR02], web
search [AFKM01, Kle99], and learning mixtures of distributions [AM05, KSV05],
and has numerous other applications.

A large body of recent work has looked at relative-error low-rank
approximation, whereby given an $n\times n$ matrix $\mathbf{A}$, an accuracy
parameter $\epsilon>0$, and a rank parameter $k$, one seeks to output a
rank-$k$ matrix $\mathbf{B}$ for which

$\displaystyle\|\mathbf{A}-\mathbf{B}\|_{F}^{2}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2},$
(1.1)

where for a matrix $\mathbf{C}$,
$\|\mathbf{C}\|_{F}^{2}=\sum_{i,j}\mathbf{C}_{i,j}^{2}$, and $\mathbf{A}_{k}$
denotes the best rank-$k$ approximation to $\mathbf{A}$ in Frobenius norm.
$\mathbf{A}_{k}$ can be computed exactly using the singular value
decomposition, but takes time $O(n^{\omega})$, where $\omega$ is the matrix
multiplication constant. We refer the reader to the survey [Woo14] and
references therein.

For worst-case matrices, it is not hard to see that any algorithm achieving
(1.1) must spend at least $\Omega(\textrm{nnz}(\mathbf{A}))$ time, where
$\textrm{nnz}(\mathbf{A})$ denotes the number of non-zero entries (sparsity)
of $\mathbf{A}$. Indeed, without reading most of the non-zero entries of
$\mathbf{A}$, one could fail to read a single large entry, thus making one’s
output matrix $\mathbf{B}$ an arbitrarily bad approximation.

A flurry of recent work [KP16, MW17, BW18, CLW18, Tan19, RSML18, GLT18,
IVWW19, SW19, GSLW19] has looked at the possibility of achieving sublinear
time algorithms (classical and quantum) for low-rank approximation. In
particular, Musco and Woodruff [MW17] consider the important case of positive-
semidefinite (PSD) matrices. PSD matrices include as special cases covariance
matrices, correlation matrices, graph Laplacians, kernel matrices and random
dot product models. Further, the special case where the input itself is low-
rank (PSD Matrix Completion) has applications in quantum state tomography
[GLF+10]. Subsequently, Bakshi and Woodruff [BW18] considered low-rank
approximation of the closely related family of Negative-type (Euclidean
Squared) distance matrices. Negative-type metrics include as special cases
$\ell_{1}$ and $\ell_{2}$ metrics, spherical metrics and hypermetrics, as well
as effective resistances in graphs [DL09, TD87, CRR+96, CKM+11]. Negative-type
metrics have found various applications in algorithm design and optimization
[ALN08, SS11, KMP14, MST15].

Musco and Woodruff show that it is possible to output a low-rank matrix
$\mathbf{B}$ in factored form achieving (1.1) in
$\widetilde{O}(nk/\epsilon^{2.5}+nk^{\omega-1}/\epsilon^{2(\omega-1)})$ time,
while reading only $\widetilde{O}(nk/\epsilon^{2.5})$ entries of $\mathbf{A}$.
They also showed a lower bound that any algorithm achieving (1.1) must read
$\Omega(nk/\epsilon)$ entries, and closing the gap between these bounds has
remained an open question. Similarly, Bakshi and Woodruff exploit the
structure of Negative-type metrics to reduce to the PSD case and obtain a bi-
criteria algorithm that requires $\widetilde{O}(nk/\epsilon^{2.5})$ queries.
The gap in the sample complexity and the requirement of a bi-criteria
guarantee remained open. We resolve these questions here, and describe our
novel technical contributions in Section 2.

Next we consider PSD matrices that have been corrupted by a small amount of
noise. A drawback of algorithms achieving (1.1) is that they cannot tolerate
any amount of unstructured noise. For instance, if one slightly corrupts a few
off-diagonal entries, making the input matrix $\mathbf{A}$ no longer PSD, then
it is impossible to detect such corruptions in sublinear time, making the
relative-error guarantee (1.1) information-theoretically impossible. Motivated
by this, we also introduce a new framework where an adversary corrupts the
input by adding a noise matrix $\mathbf{N}$ to a psd matrix $\mathbf{A}$. We
assume that the Frobenius norm of the corruption is bounded relative to the
Frobenius norm of $\mathbf{A}$, i.e.,
$\|\mathbf{N}\|^{2}_{F}\leq\eta\|\mathbf{A}\|^{2}_{F}$. We also assume the
corruption is well-spread, i.e., each row of $\mathbf{N}$ has
$\ell^{2}_{2}$-norm at most a fixed constant factor larger than
$\ell^{2}_{2}$-norm of the corresponding row of $\mathbf{A}$.

This model captures small perturbations to PSD matrices that we may observe in
real-world datasets, as a consequence of round-off or numerical errors in
tasks such as computing Laplacian pseudoinverses, and systematic measurement
errors when computing a covariance matrix. One important application captured
by our model is low-rank approximation of corrupted correlation matrices.
Finding a low-rank approximation of such matrices occurs when measured
correlations are asynchronous or incomplete, or when models are stress-tested
by adjusting individual correlations. Low-rank approximation of correlation
matrices also has many applications in finance [Hig02].

Given that it is information-theoretically impossible to obtain the relative-
error guarantee (1.1) in the _robust model_ , we relax our notion of
approximation to the following well-studied additive-error guarantee:

$\displaystyle\|\mathbf{A}-\mathbf{B}\|_{F}^{2}\leq\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}+(\epsilon+\eta)\|\mathbf{A}\|_{F}^{2}.$
(1.2)

This additive-error guarantee was introduced by the seminal work of Frieze et.
al. [FKV04], and triggered a long line of work on low-rank approximation from
a computational perspective. Frieze et al. showed that it is possible to
achieve (1.2) in $O(\textrm{nnz}(\mathbf{A}))$ time. Further, given access to
an oracle for computing row norms of $\mathbf{A}$, 1.2 is achievable in
sublinear time. More recently, the same notion of approximation was used to
obtain sublinear sample complexity and running time algorithms for distance
matrices [BW18, IVWW19], and a quantum algorithm for recommendation systems
[KP16], which was subsequently dequantized [Tan19].

This raises the question of how robust are our sublinear low-rank
approximation algorithms for structured matrices, if we relax to additive-
error guarantees and allow for corruption. In particular, can we obtain
additive-error low-rank approximation algorithms for PSD matrices that achieve
sublinear time and sample complexity in the presence of noise? We characterize
when such robust algorithms are achievable in sublinear time.

###  1.1 Our Results

Problem | Prior Work | Our Results | Query Lower
---|---|---|---
Query | Run Time | Query | Run Time | Bound
PSD LRA | $O\left(\frac{nk}{\epsilon^{2.5}}\right)$ | $O\left(\frac{nk^{\omega-1}}{\epsilon^{2\omega-2}}+\frac{nk}{\epsilon^{2.5}}\right)$ | $O^{*}\left(\frac{nk}{\epsilon}\right)$ | $O^{\dagger}\left(\frac{nk^{\omega-1}}{\epsilon^{\omega-1}}\right)$ | $\Omega\left(\frac{nk}{\epsilon}\right)$
[MW17] |  Thm. 4.1 | [MW17]
PSD LRA PSD Output | $O\left(nk\left(\frac{k}{\epsilon^{2}}+\frac{1}{\epsilon^{3}}\right)\right)$ | $O\left(nk^{\omega-1}\left(\frac{k}{\epsilon^{\omega}}+\frac{1}{\epsilon^{3\omega-3}}\right)\right)$ | $O^{*}\left(\frac{nk}{\epsilon}\right)$ | $O^{\dagger}\left(\frac{nk^{\omega-1}}{\epsilon^{\omega-1}}\right)$ | $\Omega\left(\frac{nk}{\epsilon}\right)$
[MW17] |  Thm. 4.1 | [MW17]
Negative-Type LRA | $O\left(\frac{nk}{\epsilon^{2.5}}\right)$ | $O\left(\frac{nk^{\omega-1}}{\epsilon^{2\omega-2}}+\frac{nk}{\epsilon^{2.5}}\right)$ | $O^{*}\left(\frac{nk}{\epsilon}\right)$ | $O^{\dagger}\left(\frac{nk^{\omega-1}}{\epsilon^{\omega-1}}\right)$ | $\Omega\left(\frac{nk}{\epsilon}\right)$
Bi-criteria, [BW18] |  No Bi-criteria, Thm. 4.28 | [BW18]
Coreset Ridge Regression | $O\left(\frac{ns_{\lambda}^{2}}{\epsilon^{4}}\right)$ | $O\left(\frac{ns_{\lambda}^{\omega}}{\epsilon^{\omega}}\right)$ | $O^{*}\left(\frac{ns_{\lambda}}{\epsilon^{2}}\right)$ | $O^{\dagger}\left(\frac{ns^{\omega-1}_{\lambda}}{\epsilon^{2\omega-2}}\right)$ | $\Omega\left(\frac{ns_{\lambda}}{\epsilon^{2}}\right)$
[MW17] |  Thm. 4.31 |  Thm 4.33
Table 1: Comparison with prior work. The notation $O^{*}$ and $O^{\dagger}$
represent existence of matching lower bounds for query complexity and running
time (assuming the fast matrix multiplication exponent $\omega$ is $2$)
respectively. The notation $s_{\lambda}$ is used to denote the statistical
dimension of ridge regression. All bounds are stated ignoring polylogarithmic
factors in $n,k$ and $\epsilon$.

We begin with stating our results for low-rank approximation for structured
matrices. Our main result is an optimal algorithm for low-rank approximation
of PSD matrices:

Theorem 4.1 (Informal Sample-Optimal PSD LRA.) Given a PSD matrix
$\mathbf{A}$, there exists an algorithm that queries
$\widetilde{O}(nk/\epsilon)$ entries in $\mathbf{A}$ and outputs a rank $k$
matrix $\mathbf{B}$ such that with probability $99/100$,
$\|\mathbf{A}-\mathbf{B}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$,
and the algorithm runs in time $\widetilde{O}(n\cdot(k/\epsilon)^{\omega-1})$.

######  Remark 1.1.

Our algorithm matches the sample complexity lower bound of Musco and Woodruff,
up to logarithmic factors, which shows that any randomized algorithm that
outputs a $(1+\epsilon)$-relative-error low-rank approximation for a PSD
matrix $\mathbf{A}$ must read $\Omega(nk/\epsilon)$ entries. Our running time
also improves that of Musco and Woodruff and is optimal if the matrix
multiplication exponent $\omega$ is $2$.

######  Remark 1.2.

We can extend our algorithm such that the low-rank matrix $\mathbf{B}$ we
output is also PSD with the same query complexity and running time. In
comparison, the algorithm of Musco and Woodruff accesses
$\widetilde{O}(nk/\epsilon^{3}+nk^{2}/\epsilon^{2})$ entries in $\mathbf{A}$
and runs in time
$\widetilde{O}(n(k/\epsilon)^{\omega}+nk^{\omega-1}/\epsilon^{3(\omega-1)})$.

At the core of our analysis is a sample optimal algorithm for Spectral
Regression: $\min_{\mathbf{X}}\|\mathbf{D}\mathbf{X}-\mathbf{E}\|^{2}_{2}$. We
show that when $\mathbf{D}$ has orthonormal columns and $\mathbf{E}$ is
arbitrary, we can sketch the problem by sampling rows proportional to the
leverage scores of $\mathbf{D}$ and approximately preserve the minimum cost.
This is particularly surprising since our sketch only computes sampling
probabilities by reading entries in $\mathbf{D}$, while being completely
agnostic to the entries in $\mathbf{E}$. Here, we also prove a spectral
approximate matrix product guarantee for our one-sided leverage score sketch,
which may be of independent interest. We note that such a guarantee for
leverage score sampling does not appear in prior work, and we discuss the
technical challenges we need to overcome in the subsequent section.

The techniques we develop for PSD low-rank approximation also extend to
computing a low-rank approximation for distance matrices that arise from
negative-type (Euclidean-squared) metrics. Here, our input is a pair-wise
distance matrix $\mathbf{A}$ corresponding to a point set
$\mathcal{P}=\\{x_{1},x_{2},\ldots,x_{n}\\}\in\mathbb{R}^{d}$ such that
$\mathbf{A}_{i,j}=\|x_{i}-x_{j}\|^{2}_{2}$. We obtain an optimal algorithm for
computing a low-rank approximation of such matrices:

Theorem 4.28 (Informal Sample-Optimal LRA for Negative-Type Metrics.) Given a
negative-type distance matrix $\mathbf{A}$, there exists an algorithm that
queries $\widetilde{O}(nk/\epsilon)$ entries in $\mathbf{A}$ and outputs a
rank $k$ matrix $\mathbf{B}$ such that with probability $99/100$,
$\|\mathbf{A}-\mathbf{B}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$,
and the algorithm runs in time $\widetilde{O}(n\cdot(k/\epsilon)^{\omega-1})$.

######  Remark 1.3.

Prior work of Bakshi and Woodruff [BW18] obtains a
$\widetilde{O}(nk/\epsilon^{2.5})$ query algorithm that outputs a rank-$(k+4)$
matrix $\mathbf{B}$ such that
$\|\mathbf{A}-\mathbf{B}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$.
We show that the bi-criteria guarantee is not necessary, thereby resolving an
open question in their paper.

#### Structured Regression.

The sample-optimal algorithm for PSD Low-Rank Approximation also leads to a
faster algorithm for Ridge Regression, when the design matrix is PSD. Given a
PSD matrix $\mathbf{A}$, a vector $y$ and a regularization parameter
$\lambda$, we consider the following optimization problem:
$\min_{x\in\mathbb{R}^{n}}\|\mathbf{A}x-y\|^{2}_{2}+\lambda\|x\|^{2}_{2}$.
This problem is often referred to as Ridge Regression and has been the focus
of numerous theoretical and practical works (see [Gru17] and references
therein).

Theorem 4.31 (Informal Ridge Regression.) Given a PSD matrix $\mathbf{A}$, a
regularization parameter $\lambda$ and statistical dimension
$s_{\lambda}=\mathrm{Tr}\left((\mathbf{A}^{2}+\lambda\mathbf{I})^{-1}\mathbf{A}^{2}\right)$,
there exists an algorithm that queries
$\widetilde{O}(ns_{\lambda}/\epsilon^{2})$ entries of $\mathbf{A}$ and with
probability $99/100$ outputs a $(1+\epsilon)$ approximate solution to the
Ridge Regression objective and runs in
$\widetilde{O}(n(s_{\lambda}/\epsilon^{2})^{\omega-1})$ time.

######  Remark 1.4.

Our result improves on prior work by Musco and Woodruff [MW17], who obtain an
algorithm that queries $\widetilde{O}(ns_{\lambda}^{2}/\epsilon^{4})$ entries
in $\mathbf{A}$ and runs in
$\widetilde{O}(n(s_{\lambda}/\epsilon^{2})^{\omega})$ time.

######  Remark 1.5.

Since our algorithm works for all $y$ simultaneously, we obtain a low-rank
coreset of the design matrix (in factored form) that preserves the Ridge
Regression cost up to a $(1+\epsilon)$ factor. Further, in Theorem 4.33, we
prove a matching lower bound on the query complexity for any coreset
construction.

#### Robust Low-Rank Approximation.

Next, we consider a robust form of low-rank approximation problem, where the
input is a PSD matrix corrupted by noise. In this setting, we have query
access to the corrupted matrix $\mathbf{A}+\mathbf{N}$, where $\mathbf{A}$ is
PSD and $\mathbf{N}$ is such that
$\|\mathbf{N}\|^{2}_{F}\leq\eta\|\mathbf{A}\|^{2}_{F}$. Further, for all
$i\in[n]$ $\|\mathbf{N}_{i,*}\|^{2}_{2}\leq c\|\mathbf{A}_{i,*}\|^{2}_{2}$,
for a fixed constant $c$. The diagonal of a PSD matrix carries crucial
information since the largest diagonal entry upper bounds all off-diagonal
entries. Therefore, a reasonable adversarial strategy is to corrupt the
largest diagonal entries and make them close to the small diagonal entries,
which enables the resulting matrix to have large off-diagonal entries that are
hard to find. Capturing this intuition we parameterize our algorithms and
lower bounds by the largest ratio between a diagonal entry of $\mathbf{A}$ and
$\mathbf{A}+\mathbf{N}$, denoted by
$\phi_{\max}=\max_{j\in[n]}\mathbf{A}_{j,j}/|(\mathbf{A}+\mathbf{N})_{j,j}|$.

Theorem 5.4. (Informal lower bound.) Let $\epsilon>\eta>0$. Given
$\mathbf{A}+\mathbf{N}$ such that $\mathbf{A}$ is PSD and $\mathbf{N}$ is a
corruption matrix as defined above, any randomized algorithm that with
probability at least $2/3$ outputs a rank-$k$ approximation up to additive
error $(\epsilon+\eta)\|\mathbf{A}\|^{2}_{F}$ must read
$\Omega\left(\phi_{\max}^{2}nk/\epsilon\right)$ entries of
$\mathbf{A}+\mathbf{N}$.

######  Remark 1.6.

Any algorithm must incur additive error $\eta\|\mathbf{A}\|^{2}_{F}$, since
$\mathbf{A}$ is not even identifiable below additive-error
$\eta\|\mathbf{A}\|^{2}_{F}$.

######  Remark 1.7.

In our hard instance, $\phi^{2}_{\max}$ can be as large as $\epsilon n/k$,
which implies a sample-complexity lower bound of $\Omega(n^{2})$. While this
lower bound precludes sublinear algorithms for arbitrary PSD matrices, we
observe that in many applications $\phi_{\max}$ can be significantly smaller.
For instance, if $\mathbf{A}$ is a correlation matrix, we know that the true
diagonal entries of $\mathbf{A}+\mathbf{N}$ are $1$ and can ignore any
corruption on them to bound $\phi_{\max}$ by $1$.

Motivated by the aforementioned observation, we introduce algorithms for
robust low-rank approximation, parameterized by the corruption on the diagonal
entries. We obtain the following theorem:

Theorem 5.11 (Informal Robust LRA.) Given $\mathbf{A}+\mathbf{N}$, which
satisfies our noise model, there exists an algorithm that queries
$\widetilde{O}(\phi_{\max}^{2}nk/\epsilon)$ entries in $\mathbf{A}+\mathbf{N}$
and computes a rank $k$ matrix $\mathbf{B}$ such that with probability at
least $99/100$,
$\|\mathbf{A}-\mathbf{B}\|^{2}_{F}\leq\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}+(\epsilon+\sqrt{\eta})\|\mathbf{A}\|^{2}_{F}$.

######  Remark 1.8.

While the sample complexity of this algorithm matches the sample complexity in
the lower bound, it incurs additive-error $\sqrt{\eta}\|\mathbf{A}\|^{2}_{F}$
as opposed to $\eta\|\mathbf{A}\|^{2}_{F}$. An interesting open question here
is whether we can achieve additive-error
$o(\sqrt{\eta}\|\mathbf{A}\|^{2}_{F})$, though we note that when
$\eta^{2}\leq\epsilon$, this just changes the additive error guarantee of our
low-rank approximation by a constant factor.

######  Remark 1.9.

Our techniques extend to low-rank approximation of correlation matrices, and
we obtain a sample complexity of $\widetilde{O}(nk/\epsilon)$, which is
optimal. In fact, the hard instance in [MW17] implies an $\Omega(nk/\epsilon)$
lower bound on the sample complexity, even in the presence of no noise.
Surprisingly, corrupting a correlation matrix does not increase the sample
complexity and only incurs an additive error of
$\sqrt{\eta}\|\mathbf{A}\|^{2}_{F}$ (see Corollary 5.18 for a formal
statement).

##  2 Technical Overview

In this section, we provide an overview of our techniques and supply intuition
for our proofs. As a first step, it is easy to see that the
$\Omega(\textsf{nnz}(\mathbf{A}))$ lower bound for general matrices does not
apply to PSD matrices, since it proceeds by hiding arbitrarily large entries.
Observe that, reading the diagonal of a PSD matrix certifies an upper bound on
all entries of the matrix and thus off-diagonals cannot be arbitrarily large.
With this intuition in mind, we focus on sublinear algorithms.

###  2.1 Sample-Optimal Low-Rank Approximation

At a high level, our algorithm consists of two stages: first, we use the
existing machinery developed by Musco and Woodruff [MW17] to obtain weak
projection-cost preserving sketches for $\mathbf{A}$. Our sketches are smaller
than those obtained by Musco and Woodruff, albeit satisfying a weaker
guarantee. Recall that such sketches reduce the dimensionality of the columns
(rows), while ensuring that the norm of all projections of rows of the form
$\mathbf{I}-\mathbf{P}$, for rank $k$ $\mathbf{P}$, space are simultaneously
preserved. Constructing such sketches for both the column and row spaces of
$\mathbf{A}$ results in a much smaller matrix, which we can afford to query.

At this point our approach diverges from that of Musco and Woodruff, since it
is not possible to follow their strategy and recover a
$(1+\epsilon)$-relative-error low-rank approximation from the weaker sketch we
constructed above. However, we show that our sketch has enough information to
extract a structured subspace (represented by an orthonormal basis) for
$\mathbf{A}$ such that the projection onto the orthogonal complement of this
subspace is comparable to the optimal low-rank approximation cost, in spectral
norm. Note, this guarantee is stronger than the span of the structured
subspace containing a low-rank approximation comparable to the optimal in
Frobenius norm, and indeed the latter does not suffice. In the second stage,
we show that we can recover a rank-$k$ matrix in the span of the structured
subspace such that it is a $(1+\epsilon)$-relative-error low-rank
approximation for $\mathbf{A}$. Further, we show that all these steps can be
performed in sublinear time and by reading only $\widetilde{O}(nk/\epsilon)$
entries in $\mathbf{A}$ (see Theorem 4.1 for a precise statement).

We begin by providing a bird’s eye view of the Musco-Woodruff algorithm and
how to adapt parts of it to obtain weak projection-cost preserving sketches
(PCPs). For ease of exposition, we ignore polylogarithmic factors in the
subsequent discussion. Their algorithm begins with computing the so-called
ridge leverage scores (Definition 4.7) for $\mathbf{A}^{1/2}$, which
approximate the ridge leverage scores of $\mathbf{A}$ up to a
$\sqrt{n/k}$-factor. The ridge leverage scores of $\mathbf{A}^{1/2}$ can be
approximated efficiently since we can compute the row norms of
$\mathbf{A}^{1/2}$ by simply reading the diagonal of $\mathbf{A}$. It is well
known [CMM17] that sampling $k/\epsilon^{2}_{0}$ columns of $\mathbf{A}$
proportional to its ridge leverage scores results in a sketch $\mathbf{C}$
that preserves the cost of all rank-$k$ projections $\mathbf{P}$:

$\|\mathbf{C}-\mathbf{P}\mathbf{C}\|^{2}_{F}=(1\pm\epsilon_{0})\|\mathbf{A}-\mathbf{P}\mathbf{A}\|^{2}_{F}$
(2.1)

In prior work, $\mathbf{C}$ is referred to as a projection-cost preserving
sketch (PCP). PCP constructions are useful since a low-rank approximation for
$\mathbf{C}$ translates to a low-rank approximation for $\mathbf{A}$, while
$\mathbf{C}$ has much smaller dimension. Observe that oversampling columns of
$\mathbf{A}$ proportional to the ridge leverage scores of $\mathbf{A}^{1/2}$,
by a $\sqrt{n/k}$ factor, suffices to obtain the guarantee of 2.1 (see Lemma
4.11 for a precise statement). Note, $\mathbf{C}$ may have
$\Omega(n^{1.5}/\epsilon^{2}_{0})$ non-zeros but the algorithm need not query
any entries in $\mathbf{C}$ at this stage. Musco and Woodruff then construct a
row PCP for $\mathbf{C}$ by sampling $\sqrt{nk}/\epsilon^{2.5}_{0}$ rows of
$\mathbf{C}$ proportional to the rank-$(k/\epsilon_{0})$ ridge leverage scores
of $\mathbf{A}$. The resulting matrix $\mathbf{R}$ is a
$\sqrt{nk}/\epsilon^{2.5}_{0}\times\sqrt{nk}/\epsilon^{2.5}_{0}$ matrix such
that for any rank-$k$ projection $\mathbf{P}$,

$\|\mathbf{R}-\mathbf{R}\mathbf{P}\|^{2}_{F}+O(\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F})=(1\pm\epsilon_{0})\|\mathbf{C}-\mathbf{C}\mathbf{P}\|^{2}_{F}$
(2.2)

Since $\mathbf{R}$ is a much smaller matrix, they run an input-sparsity time
algorithm to compute a low-rank approximation for it [CW13]. Using standard
regression techniques (described in [CMM17, MW17, BW18]) along with equations
2.1 and 2.2, setting $\epsilon_{0}=\epsilon$ results in a $(1+\epsilon)$-low-
rank approximation of $\mathbf{A}$ by querying $O(nk/\epsilon^{5})$ entries.
Musco and Woodruff instead use a more complicated algorithm to get a
$1/\epsilon^{2.5}$ dependence.

Diverging from Musco and Woodruff, we begin by showing that a Spectral-
Frobenius projection for $\mathbf{A}$ suffices to obtain a low-rank
approximation with $O(nk/\epsilon)$ queries. A projection matrix $\mathbf{P}$
is a $(\epsilon,k)$ Spectral-Frobenius projection for a matrix $\mathbf{A}$ if
$\|\mathbf{A}-\mathbf{A}\mathbf{P}\|_{2}^{2}\leq\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|_{F}^{2}$
(see Definition 4.13). Assuming we are handed a $(k/\epsilon)$-dimensional
structured subspace such that the projection on to this subspace is given by
$\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$, where $\mathbf{Q}$ is a $n\times
k/\epsilon$ matrix with orthonormal columns satisfying
$\|\mathbf{A}-\mathbf{Q}\mathbf{Q}^{\top}\mathbf{A}\|^{2}_{2}\leq\epsilon/k\cdot\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}$,
we show that we can extract a rank-$k$ relative-error approximation by reading
only $nk/\epsilon$ entries in $\mathbf{A}$ (Theorem 4.16 in Section 4.1). We
provide an overview of the proof here.

The SF projection property implies
$\|\mathbf{A}-\mathbf{Q}\mathbf{Q}^{\top}\mathbf{A}_{k}\mathbf{Q}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$.
Therefore, it suffices to solve the following optimization problem:

$\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}$ (2.3)

since $\mathbf{X}=\mathbf{Q}^{\top}\mathbf{A}_{k}\mathbf{Q}$ is always
feasible. While we are now optimizing over a $k/\epsilon\times k/\epsilon$
matrix $\mathbf{X}$, with rank at most $k$, the problem still seems
intractable to solve optimally in sublinear time and queries to $\mathbf{A}$.
The key idea here is that $\mathbf{Q}$ has orthonormal columns and thus the
leverage scores are precomputed for us. We can then sample columns and rows
proportional to the leverage scores of $\mathbf{Q}$ and consider a
significantly smaller sketched problem. Therefore, we create sampling matrices
$\SS$ and $\mathbf{T}$ that sample $\textrm{poly}(k/\epsilon)$ rows
proportional to the leverage scores of $\mathbf{Q}$ and consider the resulting
optimization problem:

$\min_{\textrm{rank}(X)\leq
k}\|\SS\mathbf{A}\mathbf{T}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}$
(2.4)

Here, we are faced with an intriguing phenomenon: our sketched optimization
problem does not have the property that the minimum cost for Equation 2.4 is a
$(1+\epsilon)$-approximation to the minimum cost for Equation 2.3. The reason
is that our sketch incurs a fixed additive shift term, which we cannot
approximate in sublinear time. We note that this is the bottleneck in
approximating the cost of the optimal low-rank approximation, and as mentioned
in [MW17], it is open to estimate this cost in $o(n^{3/2})$ time.

However, we can apply the structural result in Lemma 4.5 twice, to show that
the optimal solution to Equation 2.4, when plugged in to Equation 2.3 obtains
a $(1+\epsilon)$-approximation to the minimum cost. Formally, $\SS$ and
$\mathbf{T}$ have the property that if
$\widehat{\mathbf{X}}=\arg\min_{\mathbf{X}}\|\SS\mathbf{A}\mathbf{T}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}$,
then

$\|\mathbf{A}-\mathbf{Q}\widehat{\mathbf{X}}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+O(\epsilon))\min_{\textrm{rank}(X)\leq
k}\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}.$

The optimization problem in Equation 2.4 is called Generalized Low-Rank
Approximation and admits a closed form solution [FT07] (Theorem 4.15).
Further, since the problem now has all dimensions independent of $n$, we can
afford to explicitly compute $\SS\mathbf{A}\mathbf{T}$ by querying the
corresponding entries in $\mathbf{A}$. The resulting closed-form solution can
also be computed in $\textrm{poly}(k/\epsilon)$ time (and queries) which only
contributes a lower order term. We obtain one factor for the low-rank
approximation for $\mathbf{A}$ by simply computing an orthonormal basis for
$\mathbf{Q}\widehat{\mathbf{X}}$. In order to compute the second factor, we
set up and approximately solve a regression problem, the details of which can
be found in Algorithm 4.1. Efficiently solving such a regression problem is
now standard in low-rank approximation literature [CMM17, MW17, BW18].
Therefore, we can output a low-rank approximation to $\mathbf{A}$ by querying
only $\widetilde{O}(nk/\epsilon)$ entries.

Now, we focus on obtaining a Spectral-Frobenius projection for $\mathbf{A}$.
Our starting point here is to observe that the PCP construction from equations
(2.1) and (2.2) above allows to preserve the projection of columns of
$\mathbf{A}$ on all $(k/\epsilon)$-dimensional subspaces, albeit up to a
constant factor by setting $\epsilon_{0}=O(1)$. Therefore, a natural approach
is to set the error parameter $\epsilon_{0}$ in the PCP constructions to be a
small fixed constant, say $0.1$, and the rank parameter $k$ to be
$k/\epsilon$, where $\epsilon$ is the desired input accuracy. Further, we
observe that the guarantee obtained in Equation 2.1 can be strengthen to a
mixed Spectral-Frobenius PCP guarantee (also introduced by [MW17]): for all
rank-$(k/\epsilon)$ projection matrices $\mathbf{P}$, the column PCP
$\mathbf{C}$ satisfies :

$(1-0.1)\|\mathbf{A}-\mathbf{P}\mathbf{A}\|^{2}_{2}-\frac{\epsilon}{10k}\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}\leq\|\mathbf{C}-\mathbf{P}\mathbf{C}\|^{2}_{2}\leq(1+0.1)\|\mathbf{A}-\mathbf{P}\mathbf{A}\|^{2}_{2}+\frac{\epsilon}{10k}\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}$
(2.5)

For a formal statement refer to Lemma 4.12. Sampling rows of $\mathbf{C}$
proportional to the same distribution results in a row PCP for $\mathbf{C}$
such that for all rank-$(k/\epsilon)$ projections $\mathbf{P}$,

$(1-0.1)\|\mathbf{C}-\mathbf{C}\mathbf{P}\|^{2}_{2}-\frac{\epsilon}{10k}\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}\leq\|\mathbf{R}-\mathbf{R}\mathbf{P}\|^{2}_{2}\leq(1+0.1)\|\mathbf{C}-\mathbf{C}\mathbf{P}\|^{2}_{2}+\frac{\epsilon}{10k}\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}$
(2.6)

We then use an input-sparsity spectral-low-rank approximation algorithm by
[CEM+15] (Lemma 4.6), to obtain a low-dimensional subspace, represented by a
$\sqrt{nk/\epsilon}\times k/\epsilon$ matrix $\mathbf{Z}$ with orthonormal
columns such that

$\|\mathbf{R}-\mathbf{R}\mathbf{Z}\mathbf{Z}^{\top}\|^{2}_{2}\leq\frac{\epsilon}{k}\|\mathbf{R}-\mathbf{R}_{k/\epsilon}\|^{2}_{F}$
(2.7)

Observe the projection matrix $\mathbf{Z}\mathbf{Z}^{\top}$ as a
$(\epsilon,k)$-Spectral-Frobenius projection for $\mathbf{R}$. Now, using the
fact that $\mathbf{R}$ is a Spectral-Frobenius PCP for $\mathbf{C}$ and
plugging in $\mathbf{P}=\mathbf{Z}\mathbf{Z}^{\top}$ in Equation 2.6, we can
bound $\|\mathbf{C}-\mathbf{C}\mathbf{Z}\mathbf{Z}^{\top}\|^{2}_{2}$ as
follows :

$\begin{split}\|\mathbf{C}-\mathbf{C}\mathbf{Z}\mathbf{Z}^{\top}\|^{2}_{2}\leq\frac{10}{9}\|\mathbf{R}-\mathbf{R}\mathbf{Z}\mathbf{Z}^{\top}\|^{2}_{2}+\frac{\epsilon}{9k}\|\mathbf{R}-\mathbf{R}_{k/\epsilon}\|^{2}_{F}&\leq\frac{10\epsilon}{k}\|\mathbf{R}-\mathbf{R}_{k}\|^{2}_{F}+\frac{\epsilon}{9k}\|\mathbf{R}-\mathbf{R}_{k/\epsilon}\|^{2}_{F}\\\
&\leq
O\left(\frac{\epsilon}{k}\right)\|\mathbf{C}-\mathbf{C}_{k}\|^{2}_{F}\end{split}$
(2.8)

where the second inequality follows from Equation 2.7 and the third follows
from the fact that PCPs preserve Frobenius low-rank approximation cost up to a
constant factor, i.e.,
$\|\mathbf{R}-\mathbf{R}_{k/\epsilon}\|^{2}_{F}=\Theta(\|\mathbf{C}-\mathbf{C}_{k/\epsilon}\|^{2}_{F})=\Theta(\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F})$.
Therefore, $\mathbf{Z}\mathbf{Z}^{\top}$ is also a $(\epsilon,k)$\- Spectral-
Frobenius projection for $\mathbf{C}$. Here, we are faced with a few
challenges. First, the relative-error approximation spanned by the subspace
has rank $k/\epsilon$. Second, it is unclear how to obtain any reasonable
result for $\mathbf{A}$ from the above structural property, given that even
the dimensions of $\mathbf{Z}\mathbf{Z}^{\top}$ do not match $\mathbf{A}$.

Here, we observe that if we had a such a projection for the column-space of
$\mathbf{C}$, by Equation 2.5, it would also work for $\mathbf{A}$ and we
would be done. To this end, we consider the following optimization problem:

$\min_{\mathbf{W}\in\mathbb{R}^{n\times
k^{\prime}}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}$ (2.9)

We show that an orthonormal basis $\mathbf{Q}$ for an approximate minimizer to
Equation 2.9 is an SF projection for $\mathbf{C}$ and in turn $\mathbf{A}$
(since $\mathbf{C}$ is a column PCP for $\mathbf{A}$). Therefore, we focus on
optimizing Equation 2.9 and refer to this problem as Spectral Regression. We
note that unlike standard regression, here we minimize the Spectral (Operator)
norm. While the corresponding problem for minimizing Frobenius norm is
extensively studied and well understood, to the best of our knowledge the only
relevant related work on Spectral Regression is in the streaming model, by
Clarkson and Woodruff [CW09]. They construct an oblivious sketch, consisting
of random entries in $\\{-1,1\\}$, for Equation 2.9 that preserves the optimal
solution up to a $(1\pm\epsilon)$ factor. Unfortunately, we cannot use
oblivious sketching here, since $\mathbf{C}$ may be a dense matrix and we
cannot afford to read all of it.

Here, we emphasize that obtaining a sample-optimal algorithm for the
aforementioned Spectral Regression problem is crucial for our main algorithmic
result. Given that $\mathbf{C}$ is an $n\times\sqrt{nk/\epsilon}$ matrix, we
cannot query most of it and thus approximating its leverage scores is
infeasible. A natural approach here would be to follow the Affine Embedding
idea for Frobenius norm (refer to Lemma 4.5) and hope a similar guarantee
holds for the spectral norm as well. Here, one might hope to obtain a small
sketch that preserves the spectral norm cost of all $\mathbf{W}$ up to a
$(1\pm\epsilon)$ factor. While such a guarantee would suffice, we note that
$\mathbf{Z}$ could have rank as large as $k/\epsilon$ and we can no longer
afford a $(1+\epsilon)$-approximate affine embedding even for Frobenius norm,
without incurring a larger dependence on $\epsilon$. This precludes all known
approaches for sketching Equation 2.9 to preserve the optimal cost.

Instead, we relax the notion of approximation for our sketch. We observe that
it suffices to construct a sketch $\SS$ such that if
$\widehat{\mathbf{W}}=\arg\min_{\mathbf{W}}\|\mathbf{C}\SS-\mathbf{W}\mathbf{Z}^{\top}\SS\|^{2}_{2}$,
then

$\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}\leq
O(1)\left(\min_{\mathbf{W}\in\mathbb{R}^{n\times
k/\epsilon}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{C}-\mathbf{C}_{k/\epsilon}\|^{2}_{F}\right)$
(2.10)

Note, this is a weaker guarantee for the sketch $\SS$, since we only need to
preserve the cost of the optimal solution up to a mixed relative and additive
error. First, we observe such a guarantee suffices, since we can upper bound
the cost from Equation 2.10 by
$O(\epsilon/k)\cdot\|\mathbf{C}-\mathbf{C}_{k/\epsilon}\|^{2}_{F}$ and the
Spectral-Frobenius PCP from Equation 2.5 incurs this term anyway. In Theorem
4.18, we show that we can construct such a sketch $\SS$ satisfying Equation
2.10 by sampling $k/\epsilon$ columns of $\mathbf{C}$ proportional to the
leverage scores of $\mathbf{Z}^{\top}$.This is surprising since we completely
ignore all information about $\mathbf{C}$ and our sketch is not an oblivious
sketch.

The key technical lemma (Lemma 4.25) we prove here is a weak approximate
matrix product for $\mathbf{C}^{*}$ and $\mathbf{Z}^{\top}$ where
$\mathbf{C}^{*}=\mathbf{C}(\mathbf{I}-\mathbf{P}_{\mathbf{Z}^{\top}})$ is the
projection onto the orthogonal complement of $\mathbf{Z}^{\top}$. While
approximate matrix product has been extensively studied [DFK+04, Sar06, CW13],
even for spectral norm [CNW15], it is important to emphasize here that all
known constructions are either oblivious sketches or require sketches that are
sampled proportional to both $\mathbf{C}^{*}$ and $\mathbf{Z}^{\top}$. Since
$\mathbf{Z}^{\top}$ has no information about the spectrum of $\mathbf{C}^{*}$,
the main challenge here is to control the spectrum of
$\mathbf{C}^{*}\SS\SS^{\top}\mathbf{Z}^{\top}$.

In order to bound $\|\mathbf{C}^{*}\SS\SS^{\top}\mathbf{Z}^{\top}\|_{2}$, we
analyze how sampling columns of $\mathbf{C}^{*}$ proportional to the leverage
scores of $\mathbf{Z}^{\top}$ affects the spectrum of $\mathbf{C}^{*}$. An
important tool in our analysis is the following result by Rudelson and
Vershynin on how the spectral norm of a matrix degrades when we sample a
uniformly random subset of rows [RV07]. They show that sampling $q$ rows of a
matrix $\mathbf{M}$ uniformly at random, indexed by the set $\mathcal{Q}$,
results in a matrix $\mathbf{M}_{|\mathcal{Q}}$ such that

${\bf\mathbb{E}}\left[\big{\|}\mathbf{M}_{|\mathcal{Q}}\big{\|}_{2}\right]=O\left(\sqrt{\frac{q}{n}}\|\mathbf{M}\|_{2}+\sqrt{\log(q)}\|\mathbf{M}\|_{(n/q)}\right)$

where $\|\mathbf{A}\|_{(n/q)}$ is the average of the largest $n/q$
$\ell_{2}$-norms of columns of $\mathbf{A}$. Here, we prove that expected
spectral norm of $\mathbf{C}$ restricted to the columns sampled by $\SS$
proportional to the leverage scores of $\mathbf{Z}^{T}$ only exceeds that of a
random subset by a polylogarithmic factor. This result may be of independent
interest in applications where we would want to bound the spectrum of random
submatrices, where the rows or columns are not sampled uniformly.

Intuitively, there are two technical challenges we overcome in order to apply
the Rudelson-Vershynin result in our setting. First, a leverage score sampling
matrix $\SS$ need not sample columns uniformly at random, since we have no
control over the squared column norms of $\mathbf{Z}^{\top}$. Given that the
squared column norms of $\mathbf{Z}^{\top}$ may be lopsided, the subset of
columns we select could be far from a uniform sample in the worst case.
Second, the matrix we apply it to is not square and $\|\cdot\|_{(n/q)}$ norm
only shrinks substantially when the columns of $\mathbf{A}$ have the same
$\ell^{2}_{2}$ norm, up to a constant.

We therefore obtain a variant of Spectral norm decay for rectangular matrices,
i.e. for any $n\times m$ matrix $\mathbf{M}$ with roughly the same squared
column norms, we show that

${\bf\mathbb{E}}\left[\big{\|}\mathbf{M}_{|\mathcal{Q}}\big{\|}_{2}\right]=O\left(\sqrt{\frac{q}{n}}\|\mathbf{M}\|_{2}+\sqrt{\log(q)/b}\|\mathbf{M}\|_{(n/q)}\right)$
(2.11)

where $b=\left\lceil n/m\right\rceil$. To apply the above result, we then
partition the rows of $\mathbf{C}$ (since $\SS$ samples columns of
$\mathbf{C}$ as opposed to rows) into $\log(n)$ groups such that within each
group, all rows have roughly the same squared norm. We then analyze leverage
score sampling proportional to the column norms of $\mathbf{Z}^{\top}$ on each
group independently. We show that we can obtain a coupling between the two
random processes, namely uniform sampling and leverage score sampling, such
that we obtain a decay bound similar to Equation 2.11, up to log factors. We
describe our solution in more detail in Section 4. We note that our results
extend to outputting a low-rank PSD matrix as well.

#### Negative-Type Matrices.

We then use the techniques developed above to obtain an optimal relative-error
low-rank approximation for Negative-Type distance matrices. While arbitrary
metrics do not admit sublinear time algorithms for relative-error low-rank
approximation (see Theorem 7.1 in [BW18]) Bakshi and Woodruff provided a
sublinear time algorithm for metrics that satisfy negative-type inequalities.
They obtain a $(1+\epsilon)$-relative-error approximation, that queries
$\widetilde{O}(nk/\epsilon^{2.5})$ entries in the input. However, this
algorithm outputs a bi-criteria solution, i.e., given a negative-type matrix
$\mathbf{A}$, it outputs a rank-$(k+4)$ matrix $\mathbf{M}$ such that
$\|\mathbf{A}-\mathbf{B}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$.

The key observation they make is that negative-type metrics can be realized as
the distances corresponding to a point set $\mathcal{P}=\\{x_{1},x_{2},\ldots
x_{n}\\}$ such that
$\mathbf{A}_{i,j}=\|x_{i}-x_{j}\|^{2}_{2}=\|x_{i}\|^{2}_{2}+\|x_{j}\|^{2}_{2}-2\langle
x_{i},x_{j}\rangle$. Therefore, $\mathbf{A}$ admits the following
decomposition: $\mathbf{A}=\mathbf{R}_{1}+\mathbf{R}_{2}-2\mathbf{B}$, where
for all $j\in[n]$, $(\mathbf{R}_{1})_{i,j}=\|x_{i}\|^{2}_{2}$,
$\mathbf{R}_{2}=\mathbf{R}_{1}^{\top}$ and $\mathbf{B}$ is PSD. Observe that
query access to $\mathbf{A}$ suffices to obtain query access to $\mathbf{B}$
by simply assuming w.l.o.g. that $x_{1}$ is centered at the origin and the
$i$-th entry in the first row corresponds to $\|x_{i}\|^{2}_{2}$. Therefore,
any PSD low-rank approximation algorithm can be simulated on the matrix
$\mathbf{B}$ by only having query access to $\mathbf{A}$. Bakshi and Woodruff
show that obtaining the low-rank approximation for $\mathbf{B}$ and appending
the column span of $\mathbf{R}_{1}$ and $\mathbf{R}_{2}$ to it results in a
rank-$(k+4)$ bi-criteria approximation to $\mathbf{A}$. The bi-criteria
algorithm can be improved to $k+2$ using Cauchy’s Interlacing Theorem [Fis05]
and observing $\mathbf{R}_{1},\mathbf{R}_{2}$ are rank-$1$ updates to
$\mathbf{B}$, but this seems to be the limit of such approaches.

We show here that our SF projection framework can be used to obtain a sample-
optimal algorithm for negative-type metrics, and the bi-criteria approximation
is not necessary. Recall, from our discussion above, that an SF projection for
$\mathbf{A}$ suffices to obtain a low-rank approximation for $\mathbf{A}$. Our
key observation is that we can use the techniques we developed for PSD
matrices to obtain an SF projection, $\mathbf{Q}\mathbf{Q}^{T}$ , for
$\mathbf{B}$ (in the decomposition above), to which we append the column span
of $\mathbf{R}_{1},\mathbf{R}_{2}$ to $\mathbf{Q}$, and the resulting
projection (denoted by $\mathbf{\Omega}$) is an SF projection for
$\mathbf{A}$. To see this, observe,
$\|\mathbf{A}-\mathbf{\Omega}\mathbf{A}_{k}\mathbf{\Omega}\|^{2}_{F}=\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}+\|\mathbf{A}_{k}-\mathbf{\Omega}\mathbf{A}_{k}\mathbf{\Omega}\|^{2}_{F}+2\mathrm{Tr}\left((\mathbf{A}-\mathbf{A}_{k})(\mathbf{I}-\mathbf{\Omega})\mathbf{A}_{k}\mathbf{\Omega}\right)$.
A simple calculation using Von-Neumann’s trace inequality bounds
$\|\mathbf{A}_{k}-\mathbf{\Omega}\mathbf{A}_{k}\mathbf{\Omega}\|^{2}_{F}+2\mathrm{Tr}\left((\mathbf{A}-\mathbf{A}_{k})(\mathbf{I}-\mathbf{\Omega})\mathbf{A}_{k}\mathbf{\Omega}\right)$
by $O(k\|\mathbf{A}(\mathbf{I}-\mathbf{\Omega})\|^{2}_{2})$. Since
$\mathbf{\Omega}$ spans $\mathbf{R}_{1}$ and $\mathbf{R}_{2}$, and is an SF
projection for $\mathbf{B}$, we can bound the above cost by
$O(\epsilon/k)\|\mathbf{B}-\mathbf{B}_{k+2}\|_{F}$. It is easy to see that
$\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}=O(\|\mathbf{B}-\mathbf{B}_{k+2}\|^{2}_{F})$
and therefore, we conclude
$\|\mathbf{A}-\mathbf{\Omega}\mathbf{A}_{k}\mathbf{\Omega}\|^{2}_{F}\leq(1+O(\epsilon))\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
(see Lemma 4.29 for details). Subsequently, we use the sublinear algorithm we
developed for PSD matrices to obtain a low-rank approximation for
$\mathbf{A}$.

#### Ridge Regression.

Our techniques also naturally extend to ridge regression, when the design
matrix is PSD. This connection was originally outlined by Musco and Woodruff
and they obtain sublinear time algorithms for solving ridge regression,
parametrized by the statistical dimension $s_{\lambda}$. At a high level, we
compute a rank-$(s_{\lambda}/\epsilon^{2})$ spectral approximation to the
input and solve ridge regression on the resulting matrix, i.e., given a PSD
matrix $\mathbf{A}$, we compute a low-rank matrix $\mathbf{B}$ such that
$\|\mathbf{A}-\mathbf{B}\|^{2}_{2}\leq
O(\epsilon/k)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$. Further, we observe that
the low-rank matrix is in fact a coreset for the input as it simultaneously
preserves the cost of all $x$ and $y$.

We then obtain a matching query lower bound for constructing coresets for
ridge regression. Our lower bound proceeds by showing that a coreset can
output a low-rank approximation on the instance of Musco and Woodruff with a
stronger quadratic, rather than a linear dependence on $\epsilon$.
Intuitively, the hard instance has multiple principle submatrices of all $1$s
placed randomly over the matrix. Since a coreset simultaneously preserves the
ridge regression cost for all $x,y$, it suffices to query the coreset on
tuples of (scaled) eigenvectors and learn the positions of the blocks.
However, a priori, we do not know what the eingenvectors of $\mathbf{A}$ are.
Instead, we query the coreset on every vector with a bounded support, and pick
all vectors with small regression cost. We show that our resulting set only
contains vectors which do not overlap much on the locations of the hidden
blocks and we show this suffices.

###  2.2 Robust Low-Rank Approximation

The robustness model we consider is as follows: we begin with an $n\times n$
PSD matrix $\mathbf{A}$. An adversary is then allowed to arbitrarily corrupt
$\mathbf{A}$ by adding a perturbation matrix $\mathbf{N}$ such that
$\|\mathbf{N}\|^{2}_{F}\leq\eta\|\mathbf{A}\|^{2}_{F}$ and for all $i\in[j]$,
$\|\mathbf{N}_{i,*}\|^{2}_{2}\leq c\|\mathbf{A}_{i,*}\|^{2}_{2}$, for a fixed
constant $c$. Note, while the adversary is unrestricted in the entries of
$\mathbf{A}$ that it corrupts, the Frobenius norm of the corruption is bounded
in terms of the Frobenius norm of $\mathbf{A}$ and the corruption is well-
spread. The motivation for considering such a model is that many matrices that
we observe in practice might be close but not exactly PSD, for instance, small
perturbations to PSD matrices.

It is impossible to obtain a relative-error low-rank approximation in this
setting, since we cannot even identify $\mathbf{A}$ after querying all $n^{2}$
entries of $\mathbf{A}+\mathbf{N}$. To see this, consider the case where
$\mathbf{A}$ is rank-$k$, and observe that a relative-error algorithm requires
identifying $\mathbf{A}$ exactly. However, by querying all entries of
$\mathbf{A}+\mathbf{N}$, we can determine the row norms exactly. Therefore, we
can run the algorithm of Frieze-Kannan-Vempala [FKV04] to obtain a rank-$k$
matrix $\mathbf{X}\mathbf{Y}^{\top}$ (in factored form) such that with
probability at least $99/100$,

$\begin{split}\|\mathbf{A}+\mathbf{N}-\mathbf{X}\mathbf{Y}^{\top}\|^{2}_{F}&\leq\|\mathbf{A}+\mathbf{N}-(\mathbf{A}+\mathbf{N})_{k}\|^{2}_{F}+\epsilon\|\mathbf{A}+\mathbf{N}\|^{2}_{F}\\\
&\leq\|\mathbf{A}+\mathbf{N}-\mathbf{A}_{k}\|^{2}_{F}+\epsilon\|\mathbf{A}+\mathbf{N}\|^{2}_{F}\\\
&\leq\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}+\|\mathbf{N}\|^{2}_{F}+2\langle\mathbf{A}-\mathbf{A}_{k},\mathbf{N}\rangle+(3+\eta)\epsilon\|\mathbf{A}\|^{2}_{F}\\\
&\leq\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}+O(\epsilon+\sqrt{\eta})\|\mathbf{A}\|^{2}_{F}\end{split}$
(2.12)

where the second inequality follows from $(\mathbf{A}+\mathbf{N})_{k}$ being
the best rank-$k$ approximation to $\mathbf{A}+\mathbf{N}$ and
$\mathbf{A}_{k}$ is any other rank-$k$ matrix. The third inequality uses
$\|\mathbf{A}+\mathbf{N}\|^{2}_{F}\leq
2(\|\mathbf{A}\|^{2}_{F}+\|\mathbf{N}\|^{2}_{F})$, which follows from
$\ell^{2}_{2}$ distance satisfying triangle-inequality up to a factor of $2$.
The last inequaliity uses Cauchy-Schwarz on
$2|\langle\mathbf{A}-\mathbf{A}_{k},\mathbf{N}\rangle|\leq
2\|\mathbf{A}\|_{F}\cdot\|\mathbf{N}\|_{F}\leq
2\sqrt{\eta}\|\mathbf{A}\|^{2}_{F}$, which follows from the assumption on
$\mathbf{N}$. Additionally

$\begin{split}\|\mathbf{A}+\mathbf{N}-\mathbf{X}\mathbf{Y}^{\top}\|^{2}_{F}&=\|\mathbf{A}-\mathbf{X}\mathbf{Y}^{\top}\|^{2}_{F}+\|\mathbf{N}\|^{2}_{F}+2\langle\mathbf{A}-\mathbf{X}\mathbf{Y}^{\top},\mathbf{N}\rangle\\\
&\geq\|\mathbf{A}-\mathbf{X}\mathbf{Y}^{\top}\|^{2}_{F}-2\sqrt{\eta}\|\mathbf{A}\|^{2}_{F}\end{split}$
(2.13)

Combining Equations 2.12 and 2.13, we have
$\|\mathbf{A}-\mathbf{X}\mathbf{Y}^{\top}\|^{2}_{F}\leq\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}+O(\epsilon+\sqrt{\eta})\|\mathbf{A}\|^{2}_{F}$.
While this algorithm is far from optimal in terms of sample complexity, it
indicates that relaxing our guarantees to additive-error is amenable to robust
algorithms and indicates why we pick up a $\sqrt{\eta}$ term. The central
question we focus on in this section is whether there exists a robust
sublinear time and query algorithm to obtain an additive-error low-rank
approximation for PSD matrices.

We begin by showing a sample complexity lower bound if $\mathbf{A}$ is an
arbitrary PSD matrix. The intuition from the relative-error setting still
applies and the diagonal entries are crucial for sublinear algorithms. In tune
with this intuition, the adversary corrupts large diagonal entries to decrease
their magnitude and thus obfuscate rows that contain large off-diagonal
entries. We therefore parameterize our lower bound and algorithms by the
largest ratio between a diagonal entry of $\mathbf{A}$ and
$\mathbf{A}+\mathbf{N}$, denoted by
$\phi_{\max}=\max_{j\in[n]}\mathbf{A}_{j,j}/|(\mathbf{A}+\mathbf{N})_{j,j}|$.
Recall, we obtain the following lower bound:

Theorem 5.4. (Informal lower bound.) Let $\epsilon>\eta>0$. Given
$\mathbf{A}+\mathbf{N}$ such that $\mathbf{A}$ is PSD and $\mathbf{N}$ is a
corruption matrix as defined above, any randomized algorithm, that with
probability at least $2/3$, outputs a rank-$k$ approximation up to additive
error $(\epsilon+\eta)\|\mathbf{A}\|^{2}_{F}$ must read
$\Omega\left(\phi_{\max}^{2}nk/\epsilon\right)$ entries of
$\mathbf{A}+\mathbf{N}$.

In our hard instance, we have a block matrix $\mathbf{A}$, where we place a
random $\epsilon/\eta\times\epsilon/\eta$, rank-$1$, non-contiguous block
$\mathbf{B}_{1}$ such that each entry in the block is
$\sqrt{\eta^{2}n/\epsilon}$ and the remaining matrix has $1$s on the diagonals
and $0$s everywhere else. It is easy to see this matrix is PSD. We observe
that the block $\mathbf{B}_{1}$ contributes an $\epsilon$-fraction of the
Frobenius norm of $\mathbf{A}$, and the $\ell^{2}_{2}$ norm of the diagonals
is an $\eta$-fraction of the Frobenius norm of $\mathbf{A}$. Therefore, the
adversary can afford to corrupt all the diagonal entries in $\mathbf{B}_{1}$
and set them to be $1$. Such a perturbation is feasible in our model and
successfully obfuscates the large off-diagonal entries. Note, for this
perturbation $\phi^{2}_{max}=\eta^{2}n/\epsilon$.

Let the resulting matrix be denoted by $\mathbf{A}+\mathbf{N}$. Here, we
observe any $\epsilon$-additive-error low-rank approximation cannot ignore the
block $\mathbf{B}_{1}$. Since the diagonals of $\mathbf{A}+\mathbf{N}$ now
provide no information about the off-diagonal entries, any algorithm that
correctly outputs a low-rank approximation for both $\mathbf{A}+\mathbf{N}$
and $\mathbf{I}$ must detect at least one entry in $\mathbf{B}_{1}$. Since
$\mathbf{B}_{1}$ has $\epsilon^{2}/\eta^{2}$ non-zeros, any algorithm must
query $\Omega(\eta^{2}n^{2}/\epsilon^{2})=\Omega(\phi^{2}_{\max}n/\epsilon)$
entries to detect one entry. To obtain a linear dependence on $k$, we simply
create $k$ independent copies of $\mathbf{B}_{1}$.

Robust Algorithm. Next, we focus on a robust, additive-error low-rank
approximation algorithm, where the sample complexity is parameterized by
$\phi_{\max}$. We begin by introducing a new sampling procedure to construct
projection-cost preserving sketches. Our construction is simple to state: we
sample each column proportional to the corresponding diagonal entry. Computing
these sampling probabilities exactly requires reading only $n$ entries in
$\mathbf{A}+\mathbf{N}$. We show that sampling
$\widetilde{O}\left(\phi_{\max}^{2}\sqrt{n}k^{2}/\epsilon^{2}\right)$ columns
proportional to this distribution preserves the projection of the columns of
$\mathbf{A}$ onto the orthogonal complement of any rank-$k$ subspace, up to
additive error $(\epsilon+\sqrt{\eta})\|\mathbf{A}\|^{2}_{F}$.

Theorem 5.7. (Informal Robust Column PCP.) Let $\mathbf{A}+\mathbf{N}$ be an
$n\times n$ matrix following the assumptions of our noise model. Let $k\in[n]$
and $\epsilon,\sqrt{\eta}>0$. Let $q=\\{q_{1},q_{2}\ldots q_{n}\\}$ be a
probability distribution over the columns of $\mathbf{A}$ such that
$q_{j}=(\mathbf{A}+\mathbf{N})_{j,j}/\mathrm{Tr}\left(\mathbf{A}+\mathbf{N}\right)$.
Construct $\mathbf{C}$ by sampling
$\widetilde{O}(\phi^{2}_{\max}\sqrt{n}k^{2}/\epsilon^{2})$ columns of
$\mathbf{A}+\mathbf{N}$ proportional to $q$ and rescaling appropriately. Then,
with probability at least $1-c$, for any rank-$k$ orthogonal projection
$\mathbf{X}$,

$\|\mathbf{C}-\mathbf{X}\mathbf{C}\|^{2}_{F}=\|\mathbf{A}-\mathbf{X}\mathbf{A}\|^{2}_{F}\pm(\epsilon+\sqrt{\eta})\|\mathbf{A}\|^{2}_{F}$

We note that all prior PCP constructions work in the noiseless setting. As a
comparison, the construction of Cohen et. al. [CMM17] works for arbitrary
$\mathbf{A}$, but requires $nnz(\mathbf{A})$ time and queries to compute the
approximate ridge-leverage scores of $\mathbf{A}$. Musco and Musco [MM17]
describe how to approximately compute the ridge leverage scores of
$\mathbf{A}^{1/2}$ (if $\mathbf{A}$ is PSD) using the Nystrom approximation,
where $\mathbf{A}=\mathbf{A}^{1/2}\cdot\mathbf{A}^{1/2}$. Musco and Woodruff
[MW17] use this method to compute the ridge leverage scores of
$\mathbf{A}^{1/2}$ with $\Theta(nk)$ queries and show that the this provides a
$(\sqrt{n/k})$-approximation to the ridge leverage scores of $\mathbf{A}$. We
note that the guarantees obtained by [CMM17, MW17] are relative error, as
opposed to the additive error guarantee in the theorem above. Finally, Bakshi
and Woodruff [BW18] provide an additive-error sublinear time construction for
distance matrices by sampling proportional to column norms. In all the
aforementioned constructions, computing the sampling distribution is a non-
trivial task, whereas we simply sample proportional to the diagonal entries.

We observe that we sample columns of $\mathbf{A}+\mathbf{N}$, to obtain
$\mathbf{C}$ which is an unbiased estimator for
$\|\mathbf{A}+\mathbf{N}\|^{2}_{F}$. The main technical challenge in our
construction is to relate the cost of rank-$k$ projections for the column
space of $\mathbf{A}$ to that of $\mathbf{C}$, while obtaining an optimal
dependence on $n$ and $k$. Note, while we do not obtain the correct dependence
on $\epsilon$, we do not have to explicitly compute all of $\mathbf{C}$, only
a subset of it.

We then extend the diagonal sampling algorithm to construct a robust row PCP
for the matrix $\mathbf{C}$. We note that the construction for $\mathbf{A}$
does not immediately give a row PCP for $\mathbf{C}$ since $\mathbf{C}$ is no
longer a corrupted PSD matrix or even a square matrix, and thus there is no
notion of a diagonal. Here, all previous approaches to construct a PCP with a
sublinear number of queries hit a roadblock, since the matrix $\mathbf{C}$
need not have any well-defined structure apart from being a scaled subset of
the columns of the original corrupted PSD matrix $\mathbf{A}+\mathbf{N}$.
However, we show that sampling rows of $\mathbf{C}$ proportional to the
diagonal entries of $\mathbf{A}+\mathbf{N}$ results in a row PCP for
$\mathbf{C}$.

Theorem 5.10. (Informal Robust Row PCP.) Let $\mathbf{A}+\mathbf{N}$ be an
$n\times n$ matrix corresponding to our noise model and let $\mathbf{C}$ be a
column PCP for $\mathbf{A}$ as defined above. Let $p=\\{p_{1},p_{2}\ldots
p_{n}\\}$ be a probability distribution over the rows of $\mathbf{C}$ such
that
$p_{j}=(\mathbf{A}+\mathbf{N})_{j,j}/\mathrm{Tr}\left(\mathbf{A}+\mathbf{N}\right)$.
Construct $\mathbf{R}$ by sampling
$\widetilde{O}(\phi_{\max}\sqrt{n}k^{2}/\epsilon^{2})$ rows of $\mathbf{C}$
proportional to $p$ and scaling appropriately. With probability at least
$1-c$, for any rank-$k$ orthogonal projection $\mathbf{X}$,

$\|\mathbf{R}-\mathbf{R}\mathbf{X}\|^{2}_{F}=\|\mathbf{C}-\mathbf{C}\mathbf{X}\|^{2}_{F}\pm(\epsilon+\sqrt{\eta})\|\mathbf{A}\|^{2}_{F}$

For our algorithm, we begin by constructing column and row PCPs of
$\mathbf{A}+\mathbf{N}$, to obtain a $t\times t$ matrix $\mathbf{R}$, where
$t=\widetilde{O}(\phi_{\max}\sqrt{n}k^{2}/\epsilon^{2})$. Instead of reading
the entire matrix, we uniformly sample $\epsilon^{3}t/k^{3}$ entries in each
row of $\mathbf{R}$, and query these entries. Note, this corresponds to
reading $\epsilon^{3}t^{2}/k^{3}=\widetilde{O}(\phi^{2}_{\max}nk/\epsilon)$
entries in $\mathbf{A}+\mathbf{N}$. Ideally we would want to estimate the
$\ell^{2}_{2}$ norms of each row of $\mathbf{R}$ to then use a result of
Frieze-Kannan-Vempala to obtain a low-rank approximation for $\mathbf{R}$
[FKV04]. It is well known that to recover a low-rank approximation for
$\mathbf{R}$, one can sample rows of $\mathbf{R}$ proportional to row norm
estimates, denoted by $\mathcal{Y}_{i}$ [FKV04]. As shown in [IVWW19] the
following two conditions are a relaxation of those required in [FKV04], and
suffice to obtain an additive error low-rank approximation :

  1. 1.

For all $i\in[t]$, the corresponding estimate over-estimates the row norm of
$\mathbf{R}_{i,*}$, i.e., $\mathcal{Y}_{i}\geq\|\mathbf{R}_{i,*}\|^{2}_{2}$.

  2. 2.

The sum of the over-estimates is not too much larger than the Frobenius norm
of the matrix, i.e.,
$\sum_{i\in[t]}\mathcal{Y}_{i}\leq\phi^{2}_{\max}n/t\|\mathbf{R}\|^{2}_{F}$

If the two conditions are satisfied, Frieze-Kannan-Vempala implies sampling
$s$ rows of $\mathbf{R}$ proportional to $\mathcal{Y}_{i}$ results in an
$s\times t$ matrix $\SS$ such that the row space of $\SS$ contains a good
rank-$k$ approximation, where $s=O(\phi^{2}_{\max}nk/\epsilon t)$, which
matches our desired sample complexity. We show that, unfortunately, such a
guarantee is not possible even in the uncorrupted case, where $\mathbf{N}=0$
and $\phi_{\max}=1$. Intuitively, $\mathbf{R}$ may be a sparse matrix such
that many rows have large norm, and uniform sampling cannot obtain
concentration for all such rows, as required by the aforementioned conditions.

Instead, we settle for a weaker statement, where we show that the estimator
obtained by uniform sampling in each row is accurate with $o(1)$ probability.
At a high level, we show that we can design a sampling process that is
statistically close to $\ell^{2}_{2}$ sampling described by Frieze-Kannan-
Vempala [FKV04]. We then open up the analyzes of Frieze-Kannan-Vempala and
show that our sampling process suffices to recover the low-rank approximation
guarantee. Given the flurry of recent work in quantum computing [KP16, CLW18,
Tan19, RSML18, GLT18] that uses Frieze-Kannan-Vempala $\ell_{2}^{2}$ sampling
as a key algorithmic primitive, our analysis may be of independent interest.

Lemma 5.14. (Informal Estimation of Row Norms.) Let
$\mathbf{R}\in\mathbb{R}^{t\times t}$ be the row PCP as defined above. For all
$i\in[t]$ let
$\mathcal{X}_{i}=\sum_{j\in[\epsilon^{3}t/k^{3}]}\mathcal{X}_{i,j}$ such that
$\mathcal{X}_{i,j}=k^{3}\mathbf{R}^{2}_{i,j^{\prime}}/\epsilon^{3}$ with
probability $1/t$, for all $j^{\prime}\in[t]$. Then, for all $i\in[t]$,
$\mathcal{X}_{i}=\left(1\pm 0.1\right)\|\mathbf{R}_{i,*}\|^{2}_{2}$ with
probability at least $\min(\|\mathbf{R}_{i,*}\|^{2}_{2}k/\epsilon n,1)$.

We now face two major challenges: first, the probability with which the
estimators are accurate is too small to even detect all rows with norm larger
than $\phi_{\max}^{2}n\|\mathbf{R}\|^{2}_{F}/t^{2}$, and second, there is no
small query certificate for when an estimator is accurate in estimating the
row norms. Therefore, we cannot even identify the rows where we obtain an
accurate estimate of norm.

To address the first issue, we make the crucial observation that while we
cannot estimate the norm of each row accurately, we can hope to sample the row
with the same probability as Frieze-Kannan-Vempala [FKV04]. Recall, their
algorithm requires sampling row $\mathbf{R}_{i,*}$ with probability at least
$\|\mathbf{R}_{i,*}\|^{2}_{2}/\|\mathbf{R}\|^{2}_{F}$, which matches the
probability in Lemma 5.14. Therefore, we can focus on designing a weaker
notion of identifiability, that may potentially include extra rows.

We begin by partitioning the rows of $\mathbf{R}$ into two sets. Let
$\mathcal{H}=\left\\{i\textrm{ }\big{|}\textrm{
}\|\mathbf{R}_{i,*}\|^{2}_{2}\geq\phi_{\max}^{2}n/t^{2}\|\mathbf{R}\|^{2}_{F}\right\\}$
be the set of heavy rows and $[t]\setminus\mathcal{H}$ be the remaining rows.
Note,
$|\mathcal{H}|=O\left(t^{2}/\phi_{\max}^{2}n\right)=O\left(k^{4}\log^{4}(n)/\epsilon^{4}\right)$.
We then condition on our estimator having norm at least
$\phi_{\max}^{2}n\|\mathbf{R}\|^{2}_{F}/t^{2}$. Conditioned on this event, we
sample the corresponding row of $\mathbf{R}$ with probability $1$. As before,
we want to prevent sampling too many spurious rows, but we show only a subset
of the rows in $\mathcal{H}$ satisfy this condition. This ensures we identify
rows in $\mathcal{H}$ with the right probability. For all the remaining rows,
we know the norm is at most $\phi_{\max}^{2}n/t^{2}\|\mathbf{R}\|^{2}_{F}$. We
show that uniformly sampling $\phi_{\max}^{2}n/t$ such rows suffices to
simulate row norm sampling.

We then open up the analysis of Frieze-Kannan-Vemapala to show that the above
sampling procedure suffices to bound the overall variance, resulting in a
relaxation of the conditions required to obtain an additive error low-rank
approximation to $\mathbf{R}$. Once we compute a good low-rank approximation
for $\mathbf{R}$ we can follow the approach of [CMM17, MW17, BW18], where we
set up two regression problems, and use the sketch and solve paradigm to
compute a low-rank approximation for $\mathbf{A}$, culminating in Theorem
5.11.

For corrupted correlation matrices, we observe that the true uncorrupted
matrix has all diagonal entries equal to $1$. Therefore, we can discard the
diagonal entries of $\mathbf{A}+\mathbf{N}$ and assume they are $1$. In this
case, no matter what the adversary does to the diagonal, $\phi_{\max}=1$ and
we obtain an $\widetilde{O}(nk/\epsilon)$ query algorithm that satisfies the
above guarantee. Further, we show a matching sample complexity lower bound of
$\Omega(nk/\epsilon)$, to obtain $\epsilon$-additive-error, even in the
presence of no noise.

##  3 Preliminaries and Notation

Given an $m\times n$ matrix $\mathbf{A}$ with rank $r$, we can compute its
singular value decomposition, denoted by
$\texttt{SVD}(\mathbf{A})=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, such
that $\mathbf{U}$ is an $m\times r$ matrix with orthonormal columns,
$\mathbf{V}^{\top}$ is an $r\times n$ matrix with orthonormal rows and
$\mathbf{\Sigma}$ is an $r\times r$ diagonal matrix. The entries along the
diagonal are the singular values of $\mathbf{A}$, denoted by
$\sigma_{1},\sigma_{2}\ldots\sigma_{r}$. Given an integer $k\leq r$, we define
the truncated singular value decomposition of $\mathbf{A}$ that zeros out all
but the top $k$ singular values of $\mathbf{A}$, i.e.,
$\mathbf{A}_{k}=\mathbf{U}\mathbf{\Sigma}_{k}\mathbf{V}^{\top}$, where
$\mathbf{\Sigma}_{k}$ has only $k$ non-zero entries along the diagonal. It is
well known that the truncated SVD computes the best rank-$k$ approximation to
$\mathbf{A}$ under the Frobenius norm, i.e.,
$\mathbf{A}_{k}=\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\mathbf{X}\|_{F}$. More generally, for any matrix $\mathbf{M}$,
we use the notation $\mathbf{M}_{k}$ and $\mathbf{M}_{\setminus k}$ to denote
the first $k$ components and all but the first $k$ components respectively. We
use $\mathbf{M}_{i,*}$ and $\mathbf{M}_{*,j}$ to refer to the $i^{th}$ row and
$j^{th}$ column of $\mathbf{M}$ respectively. For an $n\times n$ PSD matrix
$\mathbf{A}$, we denote the singular (eigenvalue) decomposition by
$\mathbf{U}\mathbf{\Sigma}\mathbf{U}^{\top}$. Further, since
$\mathbf{\Sigma}_{i,i}\geq 0$, let
$\mathbf{A}^{1/2}=\mathbf{U}\mathbf{\Sigma}^{\frac{1}{2}}\mathbf{U}^{\top}$ be
the square root of $\mathbf{A}$. Note that
$\mathbf{A}_{i,j}=\langle\mathbf{A}^{1/2}_{i,*},\mathbf{A}^{1/2}_{j,*}\rangle$.
By Cauchy-Schwarz, for all $i,j\in[n]$,
$\mathbf{A}^{2}_{i,j}=\langle\mathbf{A}^{1/2}_{i,*},\mathbf{A}^{1/2}_{j,*}\rangle^{2}\leq\|\mathbf{A}^{1/2}_{i,*}\|^{2}_{2}\cdot\|\mathbf{A}^{1/2}_{j,*}\|^{2}_{2}=\mathbf{A}_{i,i}\cdot\mathbf{A}_{j,j}$.
We use $\textsf{nnz}(\mathbf{A})$ to denote the number of non-zero entries
(sparsity) of $\mathbf{A}$. We use operator and spectral norm interchangeably
to denote $\|\mathbf{M}\|_{2}=\max_{\|y\|_{2}=1}\|\mathbf{M}y\|_{2}$. We also
use the notation $\mathbf{M}^{\dagger}$ to denote the Moore-Penrose
pseudoinverse.

##  4 Relative Error PSD Low-Rank Approximation

In this section, we describe our main algorithm for relative-error PSD Low-
Rank Approximation, where we query only $\widetilde{O}(nk/\epsilon)$ of the
input matrix $\mathbf{A}$. This improves the best known algorithm by Musco and
Woodruff that queries $\widetilde{O}(nk/\epsilon^{2.5})$ and matches their
query lower bound of $\Omega(nk/\epsilon)$ up to polylogarithmic factors
[MW17]. Formally, we prove the following:

######  Theorem 4.1.

(Sample-Optimal PSD Low-Rank Approximation.) Given an $n\times n$ PSD matrix
$\mathbf{A}$, an integer $k$, and $1>\epsilon>0$, Algorithm 4.3 samples
$\widetilde{O}(nk/\epsilon)$ entries in $\mathbf{A}$ and outputs matrices
$\mathbf{M},\mathbf{N}^{\top}\in\mathbb{R}^{n\times k}$ such that with
probability at least $9/10$,

$\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

Further, the algorithm runs in
$\widetilde{O}(n(k/\epsilon)^{\omega-1}+(k/\epsilon^{3})^{\omega})$ time.

We begin by defining various statistical quantities associated with a given
matrix, such as the leverage and ridge-leverage scores. The leverage score of
a given row measures the importance of this row in composing the row span.
Leverage scores have found numerous applications in regression,
preconditioning, linear programming and graph sparsification [Sar06, SS11,
LS15, CLM+15]. In the special case of graphs, they are referred to as
effective resistances.

######  Definition 4.2.

(Leverage Scores.) Given a matrix $\mathbf{M}\in\mathbb{R}^{n\times m}$, let
$\textbf{m}_{i}=\mathbf{M}_{i,*}$ be the $i$-th row of $\mathbf{M}$. Then, for
all $i\in[n]$ the $i$-th row leverage score of $\mathbf{M}$ is given by

$\tau_{i}(\mathbf{M})=\textbf{m}_{i}(\mathbf{M}^{\top}\mathbf{M})^{\dagger}\textbf{m}_{i}^{\top}$

The column leverage scores can be defined analogously. Note, in the special
case where $\mathbf{M}$ has orthonormal columns, the row leverage scores of
$\mathbf{M}$ are simply the $\ell^{2}_{2}$ norms of the rows i.e.,
$\tau_{i}(\mathbf{M})=\|\textbf{m}_{i}\|^{2}_{2}$. It is well-known that
sampling rows of a matrix proportional to the leverage scores satisfies the
subspace embedding property (Spectral Sparsification for Graphs) and leads to
faster algorithms for $\ell_{2}$-norm Regression. Recall, for an $n\times m$
matrix $\mathbf{A}$, a leverage score sampling matrix
$\SS=\mathbf{D}\mathbf{\Omega}^{\top}$, where $\mathbf{D}$ is a $t\times t$
diagonal matrix and $\mathbf{\Omega}$ is an $n\times t$ sampling matrix. For
all $j\in[t]$, select row index $i\in[n]$ with probability
$p_{i}=\tau_{i}(\mathbf{A})/\sum_{i}\tau_{i}(\mathbf{A})$ and set
$\mathbf{\Omega}_{i,j}=1$ and $\mathbf{D}_{j,j}=1/\sqrt{tp_{i}}$.

######  Lemma 4.3.

(Subspace Embedding.) Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times m}$,
$\epsilon>0$, and a leverage score sampling matrix $\SS$ with
$t=O(mlog(m)/\epsilon^{2})$ rows, with probability at least $99/100$, for all
$x\in\mathbb{R}^{m}$

$\|\SS\mathbf{A}x\|^{2}_{2}=(1\pm\epsilon)\|\mathbf{A}x\|^{2}_{2}$

This simply follows from an application of the Matrix Chernoff bound. Observe
that the sketch preserves all the singular values of $\mathbf{A}$ up to a
factor of $1\pm\epsilon$. We refer the reader to a recent survey for more
details [Woo14]. Next, we recall that leverage score sampling results in a
fast algorithm for regression.

######  Lemma 4.4.

(Fast Regression, Theorem 38 [CW13].) Given matrices
$\mathbf{A}\in\mathbb{R}^{n\times m},\mathbf{B}\in\mathbb{R}^{n\times d}$ such
that rank$(\mathbf{A})\leq r$ and $\epsilon>0$, sample
$O(r\log(r)+r/\epsilon)$ rows of $\mathbf{A},\mathbf{B}$ proportional to the
leverage scores of $\mathbf{A}$ to obtain a sketch $\SS$ such that
$\mathbf{Y}^{*}=\arg\min_{\mathbf{Y}}\|\SS\mathbf{A}\mathbf{Y}-\SS\mathbf{B}\|^{2}_{F}$.
Then, with probability at least $1-c$,

$\|\mathbf{A}\mathbf{Y}^{*}-\mathbf{B}\|^{2}_{F}\leq(1+\epsilon)\min_{\mathbf{Y}}\|\mathbf{A}\mathbf{Y}-\mathbf{B}\|^{2}_{F}$

for a fixed small constant $c$. Further, the time to compute $\mathbf{Y}^{*}$
is
$O(\textsf{nnz}(\mathbf{A})\log(r/\epsilon)+(n+d)(r/\epsilon)^{\omega-1}+\textrm{poly}(r/\epsilon))$.

Note, the terms in the running time follow from using Cohen’s construction for
OSNAP [Coh16]. Leverage score sampling matrices also approximately preserve
norms in affine spaces, which leads to faster algorithms for multi-response
regression, i.e.,
$\min_{\mathbf{X}}\|\mathbf{A}\mathbf{X}-\mathbf{B}\|^{2}_{F}$, where
$\mathbf{B}$ now has a large number of columns.

######  Lemma 4.5.

(Affine Embeddings, Theorem 39 [CW13].) Given matrices
$\mathbf{A}\in\mathbb{R}^{n\times m}$ , such that
$\textrm{rank}(\mathbf{A})=r$, and $\mathbf{B}\in\mathbb{R}^{n\times d}$, let
$\SS$ be a leverage score sampling matrix with $t=O(r/\epsilon^{2})$ rows.
Further, let $\mathbf{X}^{*}$ be the optimizer for
$\min_{\mathbf{X}}\|\mathbf{A}\mathbf{X}-\mathbf{B}\|^{2}_{F}$ and let
$\mathbf{B}^{*}=\mathbf{A}\mathbf{X}^{*}-\mathbf{B}$. Then, with probability
at least $1-c$, for all $\mathbf{X}\in\mathbb{R}^{m\times d}$

$\|\SS\mathbf{A}\mathbf{X}-\SS\mathbf{B}\|^{2}_{F}-\|\SS\mathbf{B}^{*}\|^{2}_{F}=(1\pm\epsilon)\|\mathbf{A}\mathbf{X}-\mathbf{B}\|^{2}_{F}-\|\mathbf{B}^{*}\|^{2}_{F}$

for a fixed small constant $c$.

An important application of the above lemma (which we use extensively) is to
sketch constrained regression problems, for example, when the matrix
$\mathbf{X}$ has a fixed small rank. Since affine embeddings approximately
preserve the cost of all affine spaces up to a fixed shift, this guarantee in
particular holds for $\mathbf{X}$ with small rank. Recall, an important caveat
here is that the cost of the sketched problem is not a relative-error
approximation to the cost of the original problem since we cannot estimate
$\|\mathbf{B}^{*}\|^{2}_{F}$ in general. However, the upshot here is that the
aforementioned guarantee still suffices for optimization since the fixed shift
does not change the optimizer.

The next tool we use is input-sparsity time low-rank approximation. This was
achieved by Clarkson and Woodruff [CW13] and the exact dependence on
$k,\epsilon$ was improved in subsequent works [MM13, NN13, BDN15, Coh16].
While the standard low-rank approximation guarantee achieves relative-error
under Frobenius norm, here we will require a spectral norm bound, which
follows from results of [CEM+15, CMM17].

######  Lemma 4.6.

(Input-Sparsity Spectral LRA [CEM+15, CMM17].) Given a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$, $\epsilon,\delta>0$ and
$k\in\mathbb{N}$, let $k^{\prime}=O(k/\epsilon)$. Then, there exists an
algorithm that outputs a matrix
$\mathbf{Z}^{\top}\in\mathbb{R}^{k^{\prime}\times n}$ such that with
probability at least $1-\delta$,

$\|\mathbf{A}-\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\|^{2}_{2}\leq
O\left(\frac{\epsilon}{k}\right)\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}$

in time and query complexity
$\widetilde{O}\left(\textsf{nnz}(\mathbf{A})+(n+d)\textrm{poly}(k/\epsilon\delta)\right)$.

###### Proof.

By Lemma 18 from [CEM+15] it suffices to use any obvious subspace embedding
matrix with $\epsilon=O(1)$ and $k=k/\epsilon$. Here, we use OSNAP in the
regime that requires $\widetilde{O}(k/\epsilon^{2})$ rows and sparsity
$\textrm{polylog}(k)/\epsilon$ [NN13]. Instantiating this OSNAP construction
with $\epsilon=O(1)$ and $k=k/\epsilon$ results in $\mathbf{Z}^{\top}$ with
$k/\epsilon$ rows in the desired running time. ∎

Next we define the ridge leverage scores of a matrix. The ridge leverage
scores were used as sampling probabilities in the context of linear regression
and spectral approximation [LMP13, KLM+17, AM15], and low-rank approximation
[CMM17, MW17]. Intuitively, the ridge leverage scores can be thought of as
adding a regularization term that attenuates the smaller singular directions
such that they are sampled with proportionately lower probability.

######  Definition 4.7.

(Ridge Leverage Scores.) Given a matrix $\mathbf{M}\in\mathbb{R}^{n\times m}$
and an integer $k$, let $\textbf{m}_{i}=\mathbf{M}_{i,*}$ be the $i$-th row of
$\mathbf{M}$. Then, for all $i\in[n]$, the $i$-th rank-$k$ ridge leverage
score of $\mathbf{M}$ is

$\rho^{k}_{i}(\mathbf{M})=\textbf{m}_{i}\left(\mathbf{M}^{\top}\mathbf{M}+\frac{\|\mathbf{M}-\mathbf{M}_{k}\|^{2}_{F}}{k}\mathbf{I}\right)^{\dagger}\textbf{m}^{\top}_{i}$

Since we typically use the row ridge leverage scores to define a probability
distribution over the rows and sample according to this distribution, it is
crucial that their sum is small as this controls the number of rows we would
need to sample. This follows from a straightforward calculation:

######  Lemma 4.8.

(Lemma 4 from [CMM17].) Let $\rho^{k}_{i}(\mathbf{M})$ be the $i$-th ridge
leverage score of $\mathbf{M}$. Then,

$\sum_{i\in[n]}\rho^{k}_{i}(\mathbf{M})\leq 2k$

Cohen et. al. [CMM17] show that the ridge leverage scores of a matrix can be
approximated up to a small constant in $O(\textsf{nnz}(\mathbf{A}))$ time,
however this involves reading the entire matrix $\mathbf{A}$. For the special
case of $\mathbf{A}$ being PSD, Musco and Musco [MM17] show that the ridge
leverage scores of $\mathbf{A}^{\frac{1}{2}}$ can be approximated up to a
small constant using a so-called Nyström approximation.

######  Lemma 4.9.

(Lemma 4 of [MW17].) Given a PSD matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$
and integer $k$, there exists an algorithm that accesses $O(nk\log(k/\delta))$
entries in $\mathbf{A}$ and computes
$\hat{\rho}^{k}_{i}(\mathbf{A}^{\frac{1}{2}})$ for all $i\in[n]$, such that
with probability $1-\delta$,

$\rho^{k}_{i}(\mathbf{A}^{1/2})\leq\hat{\rho}^{k}_{i}(\mathbf{A}^{1/2})\leq
3\rho^{k}_{i}(\mathbf{A}^{1/2})$

and runs in time $O(n(k\log(k/\delta))^{\omega-1})$, where $\omega$ is the
matrix multiplication exponent.

Note, while it is not known how to compute ridge leverage scores of a PSD
matrix in sublinear time, Musco and Woodruff [MW17] show that the ridge
leverage scores of $\mathbf{A}^{\frac{1}{2}}$ are a coarse approximation to
the ridge leverage scores of $\mathbf{A}$.

######  Lemma 4.10.

(Lemma 5 in [MW17].) Given a PSD matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$,
for all $i\in[n]$,

$\rho^{k}_{i}(\mathbf{A})\leq
2\sqrt{\frac{n}{k}}\rho^{k}_{i}(\mathbf{A}^{\frac{1}{2}})$

Musco and Woodruff then show that sampling columns of $\mathbf{A}$, according
to the corresponding ridge leverage scores of $\mathbf{A}^{\frac{1}{2}}$,
suffices to obtain a column projection-cost preserving sketch (PCP), if we
oversample by a $\sqrt{n/k}$ factor. Projection-cost preserving sketches were
introduced by Feldmen et. al. [FSS13] and Cohen et. al. [CEM+15] and studied
in the context of low-rank approximation in [CMM17, MW17, BW18].

######  Lemma 4.11.

(Column PCP from [MW17].) Given a PSD matrix $\mathbf{A}\in\mathbb{R}^{n\times
n}$, integer $k$ and $\epsilon>0$, for all $j\in[n]$ let
$\bar{\rho}^{k}_{j}(\mathbf{A}^{\frac{1}{2}})$ be a constant approximation to
the column-ridge leverage scores of $\mathbf{A}^{\frac{1}{2}}$. Let
$q_{j}=\bar{\rho}^{k}_{j}(\mathbf{A}^{\frac{1}{2}})/\sum_{j}\bar{\rho}^{k}_{j}(\mathbf{A}^{\frac{1}{2}})$
and let
$t=O\left(\sqrt{\frac{n}{k}}\sum_{j}\bar{\rho}^{k}_{j}(\mathbf{A}^{\frac{1}{2}})\log(k/\delta)/\epsilon^{2}\right)=O\left(\sqrt{nk}\log(k/\delta)/\epsilon^{2}\right)$.
Construct $\mathbf{C}\in\mathbb{R}^{n\times t}$ by sampling $t$ columns of
$\mathbf{A}$ and setting each one to be
$\frac{1}{\sqrt{tq_{j}}}\mathbf{A}_{*,j}$ with probability $q_{j}$. Then, with
probability $1-\delta$, for any rank-$k$ projection matrix
$\mathbf{X}\in\mathbb{R}^{n\times n}$,

$(1-\epsilon)\|\mathbf{A}-\mathbf{X}\mathbf{A}\|^{2}_{F}\leq\|\mathbf{C}-\mathbf{X}\mathbf{C}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{X}\mathbf{A}\|^{2}_{F}$

Further, such a $\mathbf{C}$ can be computed by accessing $\widetilde{O}(nk)$
entries in $\mathbf{A}$ and in time $O(nk^{\omega-1})$.

This result also implies that the resulting matrix $\mathbf{C}$ is a Spectral-
Frobenius PCP for $\mathbf{A}$ (Lemma 24 in [MW17]), i.e., for any rank-$k$
projection matrix $\mathbf{X}$,

$(1-\epsilon)\|\mathbf{A}-\mathbf{X}\mathbf{A}\|^{2}_{2}-\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}\leq\|\mathbf{C}-\mathbf{X}\mathbf{C}\|^{2}_{2}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{X}\mathbf{A}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
(4.1)

As noted by Musco and Woodruff, the resulting matrix $\mathbf{C}$ is not even
square and thus is it unclear how to sample rows of $\mathbf{C}$ to obtain a
row-PCP in sublinear time and queries. In particular, the ridge-leverage
scores of rows of $\mathbf{C}$ can be an $n/k$-factor larger than the
corresponding ridge-leverage scores of $\mathbf{A}^{\frac{1}{2}}$. Instead,
Musco and Woodruff sample rows of $\mathbf{C}$ proportional to the
rank-$k/\epsilon^{2}$ ridge leverage scores of $\mathbf{A}^{\frac{1}{2}}$. In
addition, they show the stronger guarantee that a Spectral-Frobenius PCP holds
(by Lemma 8 of [MW17]) for PSD Matrices.

######  Lemma 4.12.

(Spectral-Frobenius PCP.) Given a PSD matrix $\mathbf{A}\in\mathbb{R}^{n\times
n}$, an integer $k$ and $\epsilon>0$, let $\mathbf{C}\in\mathbb{R}^{n\times
t}$ be a column PCP for $\mathbf{A}$, following Lemma 4.11. Let
$k^{\prime}=k/\epsilon^{2}$. For all $i\in[n]$, let
$\bar{\rho}^{k^{\prime}}_{i}(\mathbf{A}^{\frac{1}{2}})$ be a constant
approximation to the rank-$k^{\prime}$ row ridge leverage scores of
$\mathbf{A}^{\frac{1}{2}}$. Let
$p_{i}=\bar{\rho}^{k^{\prime}}_{i}(\mathbf{A}^{\frac{1}{2}})/\sum_{i}\bar{\rho}^{k^{\prime}}_{i}(\mathbf{A}^{\frac{1}{2}})$
and let
$t=O\left(\sqrt{\frac{n}{k}}\sum_{i}\bar{\rho}^{k^{\prime}}_{i}(\mathbf{A}^{\frac{1}{2}})\log(n)/\epsilon\right)=O\left(\sqrt{nk}\log(n)/\epsilon^{3}\right)$.
Then, with probability $1-c$, for all rank-$k^{\prime}$ projection matrices
$\mathbf{X}$,

$(1-\epsilon)\|\mathbf{C}-\mathbf{C}\mathbf{X}\|^{2}_{2}-\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}\leq\|\mathbf{R}-\mathbf{R}\mathbf{X}\|^{2}_{2}\leq(1+\epsilon)\|\mathbf{C}-\mathbf{C}\mathbf{X}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

We observe that we could compute a low-rank approximation to $\mathbf{R}$ in
input sparsity time, which already requires querying
$\Omega(\textsf{nnz}(\mathbf{R}))=\Omega(nk/\epsilon^{4})$ entries in
$\mathbf{A}$ and is far from optimal in terms of the dependence on $\epsilon$.
It is here that we digress from the approach of Musco and Woodruff. We observe
that the dependence on $n$ and $k$ is optimal and thus we instantiate the
aforementioned column and row PCPs with $\epsilon=O(1)$ and $k=k/\epsilon$.
While this results in weaker PCP guaranteees, the resulting matrix
$\mathbf{R}$ is a $\sqrt{nk/\epsilon}\times\sqrt{nk/\epsilon}$ matrix and we
can now afford to read all of it and thus we can compute a rank-$k$ low-rank
approximation to $\mathbf{R}$ using the input sparsity time algorithm of
Clarkson and Woodruff [CW13].

However, the main technical challenge here is that we can no longer use the
approach of [CMM17, MW17, BW18] to use the low-rank approximation for
$\mathbf{R}$ and solve regression problems to recover an
$\epsilon$-approximate low-rank matrix for $\mathbf{A}$. In particular, we can
now only hope for an $O(1)$ approximation if we use the standard technique of
iteratively solving regression problems. Our first insight is that computing a
Spectral Low-Rank Approximation to $\mathbf{R}$ results in a structured
projection matrix for $\mathbf{C}$, from which we can compute a structured
projection matrix for $\mathbf{A}$. Further, this structured projection can be
computed with only $\widetilde{O}(nk/\epsilon)$ queries. We first describe how
this structured projection matrix for $\mathbf{A}$ results in an efficient
low-rank approximation algorithm.

###  4.1 Structured Projections to Low-Rank Approximation

Our starting point is a structural result based on the Spectral-Frobenius
projection (SF) property introduced by Clarkson and Woodruff in the context of
approximating arbitrary matrices with low-rank PSD matrices [CW17]. In this
subsection, we show that if we are given a projection matrix that satisfies
the SF property, we can obtain a query-optimal algorithm for Low-Rank
Approximation. We begin by defining this property:

######  Definition 4.13.

($(\epsilon,k)$-SF Projection.) Given any matrix
$\mathbf{A}\in\mathbb{R}^{n\times n}$, integer $k$, and $\epsilon>0$, a
projection matrix $\mathbf{P}\in\mathbb{R}^{n\times n}$ is $(\epsilon,k)$-SF
w.r.t. $\mathbf{A}$ if

$\|\mathbf{A}-\mathbf{A}\mathbf{P}\|^{2}_{2}\leq\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

or

$\|\mathbf{A}-\mathbf{P}\mathbf{A}\|^{2}_{2}\leq\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

Intuitively, the following structural result of Clarkson and Woodruff relates
an $(\epsilon,k)$-SF projection to a relative-error low-rank approximation. We
leverage this connection heavily in subsequent sections.

Input: A PSD Matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, integer $k$,
$\epsilon>0$, an orthonormal matrix $\mathbf{Q}\in\mathbb{R}^{n\times
k^{\prime}}$ such that the projection matrix
$\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$ satisfies
$\|\mathbf{A}-\mathbf{P}\mathbf{A}\|^{2}_{2}\leq\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
1. Consider the optimization problem: $\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}\vspace{-0.1in}$
2. For all $i\in[n]$, compute the leverage scores, $\tau_{i}(\mathbf{Q})$.
Since $\mathbf{Q}$ has orthonormal columns,
$\tau_{i}(\mathbf{Q})=\|\mathbf{Q}_{i,*}\|^{2}_{2}$ and can be computed
exactly. Let $p=\\{p_{1},p_{2}\ldots p_{n}\\}$ denote a distribution over rows
of $\mathbf{A}$ for which
$p_{i}=\tau_{i}(\mathbf{Q})/\sum_{i^{\prime}}\tau_{i^{\prime}}(\mathbf{Q})$.
3. Let $t=k^{\prime}/\epsilon^{2}$. Construct a leverage score sampling matrix
$\SS$ by sampling $t$ rows of $\mathbf{A}$, such that
$\SS=\mathbf{D}\mathbf{\Omega}^{\top}$, where $\mathbf{D}$ is a $t\times t$
diagonal matrix and $\mathbf{\Omega}$ is an $n\times t$ sampling matrix. For
all $j\in[t]$, select row index $i\in[n]$ with probability $p_{i}$ and set
$\mathbf{\Omega}_{i,j}=1$ and $\mathbf{D}_{j,j}=1/\sqrt{tp_{i}}$. Repeat this
sampling process to construct another leverage score sampling matrix
$\mathbf{T}$. 4. Consider the sketched optimization problem :
$\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\SS\mathbf{A}\mathbf{T}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}\vspace{-0.1in}$
Compute $\SS\mathbf{A}\mathbf{T}$, $\mathbf{P}_{\SS\mathbf{Q}}$,
$\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}}$,$(\SS\mathbf{Q})^{\dagger}$ and
$(\mathbf{Q}^{\top}\mathbf{T})^{\dagger}$, where $\mathbf{P}_{\SS\mathbf{Q}}$
and $\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}}$ are the projections onto
$\SS\mathbf{Q}$ and $\mathbf{Q}^{\top}\mathbf{T}$ respectively. Compute
$\textsf{SVD}(\mathbf{P}_{\SS\mathbf{Q}}\SS\mathbf{A}\mathbf{T}\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}})$.
By Theorem 4.15 the sketched problem is minimized by
$\mathbf{X}^{*}=(\SS\mathbf{Q})^{\dagger}[\mathbf{P}_{\SS\mathbf{Q}}\SS\mathbf{A}\mathbf{T}\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}}]_{k}(\mathbf{Q}^{\top}\mathbf{T})^{\dagger}$.
5. Let $\mathbf{U}^{*}\in\mathbb{R}^{k^{\prime}\times k}$ be an orthonormal
basis for the columns of $\mathbf{X}^{*}$. Compute an orthonormal basis
$\mathbf{M}$ for $\mathbf{Q}\mathbf{U}^{*}$. Consider the following regression
problem: $\min_{\mathbf{Y}\in\mathbb{R}^{k\times
n}}\|\mathbf{A}-\mathbf{M}\mathbf{Y}\|^{2}_{F}$. For all $i\in[n]$, compute
$\tau_{i}(\mathbf{M})=\|\mathbf{M}_{i,*}\|^{2}_{2}$. Let
$q=\\{q_{1},q_{2},\ldots,q_{n}\\}$ be a distribution over the rows of
$\mathbf{A}$ such that
$q_{i}=\tau_{i}(\mathbf{M})/\sum_{i^{\prime}\in[n]}\tau_{i^{\prime}}(\mathbf{M})$.
Let $\mathbf{W}$ be a leverage score sampling matrix with $k/\epsilon$ rows
sampled proportional to $q$. 6. Consider the sketched regression problem:
$\min_{\mathbf{Y}\in\mathbb{R}^{k\times
n}}\|\mathbf{W}\mathbf{A}-\mathbf{W}\mathbf{M}\mathbf{Y}\|^{2}_{F}$. Let
$\mathbf{N}$ be the minimizer to this regression problem computed using the
algorithm from Lemma 4.4. Output:
$\mathbf{M},\mathbf{N}^{\top}\in\mathbb{R}^{n\times k}$ such that
$\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

######  Lemma 4.14.

(Structured Projections and Low-Rank Approximation [CW17].) Let
$\mathbf{P}\in\mathbb{R}^{n\times n}$ be an $(\epsilon,k)$-SF projection w.r.t
$\mathbf{A}$, then

$\|\mathbf{A}-\mathbf{P}\mathbf{A}_{k}\mathbf{P}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

Ignoring computational and query complexity constraints, suppose we were given
a matrix $\mathbf{Q}\in\mathbb{R}^{n\times k^{\prime}}$ with orthonormal
columns such that $\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$ is an
$(\epsilon,k)$-SF Projection, where $k^{\prime}$ is the dimension of the space
$\mathbf{P}$ projects onto. Note, for now it suffices to set
$k^{\prime}=\textrm{poly}(k/\epsilon)$. As a consequence of Lemma 4.14, we
observe that solving the following constrained regression problem suffices to
obtain a $(1+\epsilon)$-relative error solution to the Low-Rank Approximation
problem:

$\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}$ (4.2)

However, there are several challenges pertaining to this approach. As noted
above, it is not immediately clear how to obtain such a $\mathbf{Q}$ with
$nk/\epsilon$ queries to $\mathbf{A}$. Further, it is not immediately clear
how to solve Equation 4.2 efficiently. While we have reduced to optimizing
over $k^{\prime}\times k^{\prime}$ sized matrices $\mathbf{X}$ with rank at
most $k$, the problem still seems intractable in sublinear time and queries.

We begin by describing how to solve the optimization problem in Equation 4.2
with $\widetilde{O}(nk/\epsilon)$ queries given that we have access to
$\mathbf{Q}$ and $\mathbf{Q}\mathbf{Q}^{\top}$ is an $(\epsilon,k)$-SF
Projection. At a high level, our approach is to sketch the problem by sampling
rows and columns proportional to the row leverage scores of $\mathbf{Q}$. We
observe that since $\mathbf{Q}$ has orthonormal columns, the row leverage
scores of $\mathbf{Q}$ are simply the $\ell^{2}_{2}$ norms of corresponding
rows. Therefore, we create sampling matrices $\SS$ and $\mathbf{T}$ that
sample $\textrm{poly}(k^{\prime})$ rows proportional to the leverage scores of
$\mathbf{Q}$ and consider the resulting optimization problem:

$\min_{\textrm{rank}(X)\leq
k}\|\SS\mathbf{A}\mathbf{T}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}$
(4.3)

We then show that the minimizer for Equation 4.3 is an approximate minimizer
for Equation 4.2. Further, the optimization problem in Equation 4.3 is
referred to as Generalized Low-Rank Approximation and admits a closed form
solution:

######  Theorem 4.15.

(Generalized Low-Rank Approximation [FT07].) Let
$\mathbf{A}\in\mathbb{R}^{n\times n}$, $\mathbf{B}\in\mathbb{R}^{n\times
k^{\prime}}$ and $\mathbf{C}\in\mathbb{R}^{k^{\prime}\times n}$ and
$k\in\mathbb{N}$. Then, the Generalized Low-Rank Approximation problem

$\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\mathbf{B}\mathbf{X}\mathbf{C}\|^{2}_{F}$

is minimized by
$\mathbf{X}=\mathbf{B}^{\dagger}[\mathbf{P}_{B}\mathbf{A}\mathbf{P}_{C}]_{k}\mathbf{C}^{\dagger}$,
where $\mathbf{P}_{B},\mathbf{P}_{C}$ are the projection matrices onto
$\mathbf{B}$ and $\mathbf{C}$ respectively.

We apply the above theorem to Equation 4.3. Both the query complexity and
running time here contribute a lower-order term and we can afford to compute
the SVD for each term. Let $\mathbf{X}^{*}$ be the solution to the sketched
optimization problem in Equation 4.3. Then, we can compute $\mathbf{U}^{*}$,
an orthonormal column basis for $\mathbf{X}^{*}$ and consider $\mathbf{M}$, an
orthonormal basis for $\mathbf{Q}\mathbf{U}^{*}\in\mathbb{R}^{n\times k}$ to
be one of the low-rank factors for $\mathbf{A}$. To find the second factor, we
set up the following regression problem:

$\min_{\mathbf{N}\in\mathbb{R}^{k\times
n}}\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{F}$ (4.4)

Again, $\mathbf{M}$ has orthonormal columns and thus we can efficiently
compute the corresponding row leverage scores and sample $k/\epsilon$ rows. By
Lemma 4.5 this achieves a $(1+\epsilon)$-approximation to the optimal cost in
Equation 4.4 and obtains an $\mathbf{N}^{*}$ with $\widetilde{O}(nk/\epsilon)$
queries to $\mathbf{A}$. At this stage, we have obtained a
$(1+\epsilon)$-approximate rank-$k$ solution to Equation 4.2 and Lemma 4.14
implies that we are done. We now formalize this argument:

######  Theorem 4.16.

(Structured Projection to Low-Rank Approximation.) Given a rank-$k^{\prime}$
projection matrix $\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$, such that
$\mathbf{P}$ is an $(\epsilon,k)$-SF projection, Algorithm 4.1 queries
$\widetilde{O}(nk/\epsilon+k^{\prime 2}/\epsilon^{4})$ entries in $\mathbf{A}$
and with probability $99/100$ outputs
$\mathbf{M},\mathbf{N}^{\top}\in\mathbb{R}^{n\times k}$ such that

$\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

Further, Algorithm 4.1 runs in time
$\widetilde{O}(n(k/\epsilon)^{\omega-1}+nk^{\prime\omega-1}+(k^{\prime}/\epsilon^{2})^{\omega})$.

###### Proof.

Since $\mathbf{P}$ is an $(\epsilon,k)$-SF projection and $\mathbf{A}_{k}$ is
a feasible solution to $\min_{\textrm{rank}(\mathbf{Y})\leq
k}\|\mathbf{A}-\mathbf{P}\mathbf{Y}\mathbf{P}\|^{2}_{F}$, from Lemma 4.14 we
have

$\min_{\textrm{rank}(\mathbf{Y})\leq
k}\|\mathbf{A}-\mathbf{P}\mathbf{Y}\mathbf{P}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
(4.5)

Since $\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$, we can substitute it in
Equation 4.5 to get $\min_{\textrm{rank}(\mathbf{Y})\leq
k}\|\mathbf{A}-\mathbf{Q}\mathbf{Q}^{\top}\mathbf{Y}\mathbf{Q}\mathbf{Q}^{\top}\|^{2}_{F}$.
We further relax this by optimizing over all rank-$k$ matrices
$\mathbf{X}\in\mathbb{R}^{k^{\prime}\times k^{\prime}}$ instead of matrices of
the form $\mathbf{Q}\mathbf{Y}\mathbf{Q}^{\top}$. Therefore,

$\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
(4.6)

where we are now optimizing over a $k^{\prime}\times k^{\prime}$ matrix
$\mathbf{X}$, which is considerably smaller than $\mathbf{Y}$. Let
$\SS,\mathbf{T}^{\top}\in\mathbb{R}^{k^{\prime}/\epsilon^{2}\times n}$ be the
leverage score sampling matrices as defined in Algorithm 4.1. Observe, from
Lemma 4.5 we know $\SS$ has a sufficient number of rows to be an affine
embedding for Equation 4.6. However, we cannot directly apply the affine
embedding guarantee since $\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}$
is not an affine space. Let $\mathbf{H}$ be $k^{\prime}\times n$ matrix, let
$\mathbf{H}^{*}=\arg\min_{\mathbf{H}}\|\mathbf{A}-\mathbf{Q}\mathbf{H}\|^{2}_{F}$
and let $\mathbf{A}^{*}$ be $\mathbf{A}-\mathbf{Q}\mathbf{H}^{*}$. Then, with
probability at least $1-c_{1}$, for all $\mathbf{H}$,

$\|\SS\mathbf{A}-\SS\mathbf{Q}\mathbf{H}\|^{2}_{F}-\|\SS\mathbf{A}^{*}\|^{2}_{F}=(1\pm\epsilon)\|\mathbf{A}-\mathbf{Q}\mathbf{H}\|^{2}_{F}+\|\mathbf{A}^{*}\|^{2}_{F}$
(4.7)

Since Equation 4.7 holds for all $\mathbf{H}$, in particular it holds for all
rank-$k$ matrices $\mathbf{X}$ such that
$\mathbf{H}=\mathbf{X}\mathbf{Q}^{\top}$. Therefore, with probability at least
$1-c_{1}$, for all rank $k$ matrices $\mathbf{X}$,

$\|\SS\mathbf{A}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}-\|\SS\mathbf{A}^{*}\|^{2}_{F}=(1\pm\epsilon)\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}+\|\mathbf{A}^{*}\|^{2}_{F}$
(4.8)

Here, we observe that while we cannot estimate $\|\mathbf{A}^{*}\|^{2}_{F}$
accurately, it is a fixed matrix independent of $\mathbf{X}$ and thus we can
still approximately optimize. Let $\zeta_{1}$ be the event that Equation 4.8
holds. We now use the sampling matrix $\mathbf{T}$ to sketch
$\|\SS\mathbf{A}-\mathbf{Z}\mathbf{Q}^{\top}\|^{2}_{F}$. Let
$\mathbf{Z}^{\prime}=\arg\min_{\mathbf{Z}}\|\SS\mathbf{A}-\mathbf{Z}\mathbf{Q}^{\top}\|^{2}_{F}$
and let
$\SS\mathbf{A}^{\prime}=\SS\mathbf{A}-\mathbf{Z}^{\prime}\mathbf{Q}^{\top}$.
Then, with probability at least $1-c_{2}$, for all $\mathbf{Z}$,

$\|\SS\mathbf{A}\mathbf{T}-\mathbf{Z}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}-\|\SS\mathbf{A}^{\prime}\mathbf{T}\|^{2}_{F}=(1\pm\epsilon)\|\SS\mathbf{A}-\mathbf{Z}\mathbf{Q}^{\top}\|^{2}_{F}+\|\SS\mathbf{A}^{\prime}\|^{2}_{F}$
(4.9)

In particular, the above equation holds for all rank-$k$ matrices $\mathbf{X}$
such that $\mathbf{Z}=\SS\mathbf{Q}^{\top}\mathbf{X}$. Let $\zeta_{2}$ be the
event that the aforementioned equation holds. Combining equations 4.8 and 4.9
and conditioning on $\zeta_{1}$ and $\zeta_{2}$, for all rank-$k$ matrices
$\mathbf{X}$,

$\|\SS\mathbf{A}\mathbf{T}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}-\|\SS\mathbf{A}^{\prime}\mathbf{T}\|^{2}_{F}=(1\pm\epsilon)^{2}\left(\|\mathbf{A}-\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}+\|\SS\mathbf{A}^{*}\|^{2}_{F}+\|\mathbf{A}^{*}\|^{2}_{F}\right)+\|\SS\mathbf{A}^{\prime}\|^{2}_{F}$
(4.10)

Here, we observe that while the sketch does not preserve the cost of all
$\mathbf{X}$ up to relative error $(1+\epsilon)$, the additive error
$\Delta\leq(1+\epsilon)\left(\|\SS\mathbf{A}^{*}\|^{2}_{F}+\|\mathbf{A}^{*}\|^{2}_{F}+\|\SS\mathbf{A}^{\prime}\|^{2}_{F}+\|\SS\mathbf{A}^{\prime}\mathbf{T}\|^{2}_{F}\right)$
is fixed and is independent of $\mathbf{X}$. Let
$\mathbf{X}^{*}=\arg\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\SS\mathbf{A}\mathbf{T}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\mathbf{T}\|^{2}_{F}$.
Then, union bounding over $\zeta_{1}$ and $\zeta_{2}$, with probability
$1-c_{1}-c_{2}$,

$\|\mathbf{A}-\mathbf{Q}\mathbf{X}^{*}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+\epsilon)\min_{\textrm{rank}(\mathbf{X})\leq
k}\|\mathbf{A}-\SS\mathbf{Q}\mathbf{X}\mathbf{Q}^{\top}\|^{2}_{F}$ (4.11)

Therefore, it suffices to efficiently compute $\mathbf{X}^{*}$. By Theorem
4.15, we know that the sketched optimization problem above is minimized by
$\mathbf{X}^{*}=(\SS\mathbf{Q})^{\dagger}[\mathbf{P}_{\SS\mathbf{Q}}\SS\mathbf{A}\mathbf{T}\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}}]_{k}(\mathbf{Q}^{\top}\mathbf{T})^{\dagger}$,
which can be computed exactly as shown in Step 4 of Algorithm 4.1. We note
that we can now explicitly compute $\SS\mathbf{A}\mathbf{T}$ by querying the
relevant entries in $\mathbf{A}$. Further, we can compute $\SS\mathbf{Q}$ and
$\mathbf{Q}^{\top}\mathbf{T}$ without querying $\mathbf{A}$ at all. Recalling
equation 4.11 we can approximate the optimal low rank approximation cost:

$\|\mathbf{A}-\mathbf{Q}\mathbf{X}^{*}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+O(\epsilon))\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

While we have now approximately minimized the optimization problem from
Equation 4.2, recall our goal was to obtain a rank-$k$ approximation to
$\mathbf{A}$ in factored form i.e., outputting $n\times k$ matrices
$\mathbf{M},\mathbf{N}^{\top}$ such that the low rank approximation is given
by $\mathbf{M}\mathbf{N}$. Towards this end, we compute $\mathbf{U}^{*}$, an
orthonormal column basis for $\mathbf{X}^{*}$ such that
$\mathbf{X}^{*}=\mathbf{U}^{*}\mathbf{V}^{*}$. Substituting this in the above
equation we have

$\|\mathbf{A}-\mathbf{Q}\mathbf{U}^{*}\mathbf{V}^{*}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
(4.12)

Let $\mathbf{M}=\mathbf{Q}\mathbf{U}^{*}\in\mathbb{R}^{n\times k}$ be one of
the low-rank factors for $\mathbf{A}$. To find the second factor, we observe :

$\min_{\mathbf{Y}\in\mathbb{R}^{k\times
n}}\|\mathbf{A}-\mathbf{M}\mathbf{Y}\|^{2}_{F}\leq\|\mathbf{A}-\mathbf{M}\mathbf{V}^{*}\mathbf{Q}^{\top}\|^{2}_{F}$
(4.13)

and therefore, approximately optimizing
$\min_{\mathbf{Y}\in\mathbb{R}^{k\times
n}}\|\mathbf{A}-\mathbf{M}\mathbf{Y}\|^{2}_{F}$ suffices. Again, $\mathbf{M}$
has orthonormal columns and thus we can efficiently compute the corresponding
leverage scores to create a sketch $\mathbf{W}$ with $O(k/\epsilon)$ rows.
From Lemma 4.4, with probability at least $1-c_{3}$ for all $\mathbf{Y}$,

$\|\mathbf{W}\mathbf{A}-\mathbf{W}\mathbf{M}\mathbf{Y}\|^{2}_{F}=(1\pm\epsilon)\|\mathbf{A}-\mathbf{M}\mathbf{Y}\|^{2}_{F}$

Let $\mathbf{N}$ be the optimal solution for the sketched problem as defined
in Algorithm 4.1. Then, with probability at least $1-c_{3}$,

$\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{F}\leq\left(\frac{1+\epsilon}{1-\epsilon}\right)\min_{\mathbf{Y}\in\mathbb{R}^{k\times
n}}\|\mathbf{A}-\mathbf{M}\mathbf{Y}\|^{2}_{F}$ (4.14)

We conclude correctness by union bounding over the failure probabilities of
all the sketches and observing that with probability at least $99/100$,

$\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{F}\leq(1+O(\epsilon))\|\mathbf{A}-\mathbf{M}\mathbf{V}^{*}\mathbf{Q}^{\top}\|^{2}_{F}\leq(1+O(\epsilon))\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

where the inequalities follow from Equations 4.4, 4.13 and 4.12.

Finally, we analyze the query complexity and running time of our algorithm.
Since Algorithm 4.1 is given $\mathbf{Q}$ as input, computing the leverage
scores in Step 2 requires no queries to $\mathbf{A}$ and requires
$O(nk^{\prime})$ time. Next, observe we do not have to explicitly compute
$\SS\mathbf{A}$ or $\mathbf{A}\mathbf{T}$, since $\SS\mathbf{A}\mathbf{T}$ is
simply a submatrix of $\mathbf{A}$ with $(k^{\prime}/\epsilon^{2})^{2}$
entries appropriately scaled, it suffices to query them.
$\SS\mathbf{A}\mathbf{T}$ can be computed in $O(k^{\prime 2}/\epsilon^{4})$
time. Next, we compute $\textsf{SVD}(\SS\mathbf{Q})$ and
$\textsf{SVD}(\mathbf{Q}^{\top}\mathbf{T})$, which requires no queries to
$\mathbf{A}$ and time $O(k^{\prime\omega}/\epsilon^{2})$. We can then compute
$(\SS\mathbf{Q})^{\dagger},(\mathbf{Q}^{\top}\mathbf{T})^{\dagger},\mathbf{P}_{\SS\mathbf{Q}}$
and $\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}}$ from the aforementioned SVDs.
Next, we compute the matrix
$\mathbf{P}_{\SS\mathbf{Q}}\SS\mathbf{A}\mathbf{T}\mathbf{P}_{\mathbf{Q}^{\top}\mathbf{T}}$,
which requires no extra queries to $\mathbf{A}$ and time
$O((k^{\prime}/\epsilon^{2})^{\omega})$, which is also the time required to
compute its SVD. We can then compute $\mathbf{X}^{*}$ in Step 4 with a total
of $O(k^{\prime 2}/\epsilon^{4})$ queries to $\mathbf{A}$ in time
$O(nk^{\prime}+(k^{\prime}/\epsilon^{2})^{\omega})$.

In Step 5, we can compute $\mathbf{U}^{*}$ by computing the SVD of
$\mathbf{X}^{*}$ and compute $\mathbf{M}$ in time
$O(nk^{\prime\omega-1}+k^{\prime\omega})$ and do not require any queries to
$\mathbf{A}$. In Step 6, computing $\mathbf{W}\mathbf{A}$ requires
$O(nk/\epsilon)$ queries to $\mathbf{A}$, since $\mathbf{W}$ has
$\widetilde{O}(k/\epsilon)$ rows. Note, this step contributes the leading term
to the query complexity and it is crucial $\mathbf{W}$ does not have more
rows. By Lemma 4.4, $\mathbf{N}$ can be computed in time
$\widetilde{O}(nk/\epsilon+n(k/\epsilon)^{\omega-1}+k^{3}/\epsilon)$. Overall,
Algorithm 4.1 requires $\widetilde{O}(nk/\epsilon+k^{\prime 2}/\epsilon^{4})$
queries to $\mathbf{A}$ and runs in time
$\widetilde{O}(n(k/\epsilon)^{\omega-1}+nk^{\prime\omega-1}+(k^{\prime}/\epsilon^{2})^{\omega})$.
∎

In light of Theorem 4.16, to obtain a low rank approximation for $\mathbf{A}$,
it suffices to obtain an SF Projection. In particular, it suffices to obtain a
matrix $\mathbf{Q}\in\mathbb{R}^{n\times k^{\prime}}$, for
$k^{\prime}=\textrm{poly}(k,1/\epsilon)$ such that
$\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$ is an $(\epsilon,k)$-SF projection,
by querying $\widetilde{O}(nk/\epsilon)$ entries in $\mathbf{A}$. One possible
approach to computing such a $\mathbf{Q}$ is to use the following result by
Musco and Woodruff:

######  Theorem 4.17.

(Theorem 25, [MW17].) Given a PSD matrix $\mathbf{A}$, integer $k$,
$\epsilon>0$, there exists an algorithm that reads
$\widetilde{O}(nk/\epsilon^{6}+nk^{2}/\epsilon^{2})$ entries of $\mathbf{A}$
and with probability at least $99/100$, outputs
$\mathbf{M},\mathbf{N}^{\top}\in\mathbb{R}^{n\times k}$ such that

$\|\mathbf{A}-\mathbf{M}\mathbf{N}\|^{2}_{2}\leq(1+\epsilon)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$

Instantiating this theorem with $\epsilon=O(1)$ and $k=k/\epsilon$, we obtain
a matrix $\mathbf{M},\mathbf{N}^{\top}\in\mathbb{R}^{n\times k/\epsilon}$ such
that

$\begin{split}\|\mathbf{A}-\mathbf{M}\mathbf{N}^{\top}\|^{2}_{2}&\leq
O(1)\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{2}+O\left(\frac{\epsilon}{k}\right)\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}\\\
&\leq
O\left(\frac{\epsilon}{k}\right)\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}\end{split}$

where the last inequality follows from observing
$\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{F}\leq\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}$
and

$\|\mathbf{A}-\mathbf{A}_{k}\|^{2}_{F}=\sum^{n}_{j=k+1}\sigma^{2}_{j}(\mathbf{A})\geq\left(\frac{k}{\epsilon}-k\right)\sigma^{2}_{k/\epsilon}\geq\left(\frac{k}{\epsilon}-k\right)\|\mathbf{A}-\mathbf{A}_{k/\epsilon}\|^{2}_{2}$

We can then compute an orthonormal basis for $\mathbf{M}$ and denote it by
$\mathbf{Q}$. Here, we observe $\mathbf{P}=\mathbf{Q}\mathbf{Q}^{\top}$ is an
$(\epsilon,k)$-SF projection matrix. Further, the algorithm of Musco and
Woodruff instantiated with the above parameters queries
$\widetilde{O}(nk^{2}/\epsilon^{2})$ entries in $\mathbf{A}$. As a corollary
of Theorem 4.16, providing the rank-$k/\epsilon$ projection matrix
$\mathbf{Q}$ as input to Algorithm 4.1, implies an algorithm for low rank
approximation which queries $\widetilde{O}(nk^{2}/\epsilon^{2})$ entries in
$\mathbf{A}$. This already improves the $\epsilon$-dependence in the query
complexity of best known algorithm for PSD low-rank approximation, since the
algorithm of Musco and Woodruff requires $\widetilde{O}(nk/\epsilon^{2.5})$
queries [MW17]. Note, this algorithm has worse dependence on $k$. However, our
goal is to obtain linear dependence on both $k$ and $1/\epsilon$. Towards this
end, we focus on obtaining an SF projection with fewer queries to
$\mathbf{A}$.

###  4.2 Spectral Regression

In this subsection, we consider the Spectral Regression problem. This problem
is a natural generalization of least-squares regression, when the response
variable is a matrix. Spectral Regression arises in the context of Regularized
Least Squares Classification, for instance [CLL+10]. Given matrices
$\mathbf{A}\in\mathbb{R}^{n\times d}$, $\mathbf{X}\in\mathbb{R}^{d\times m}$
and $\mathbf{B}\in\mathbb{R}^{n\times m}$, the Spectral Regression problem
considers the following optimization problem:

$\min_{\mathbf{X}}\|\mathbf{A}\mathbf{X}-\mathbf{B}\|_{2}$

We note that this is natural variant of multi-response regression where we
minimize the difference between $\mathbf{A}\mathbf{X}$ and $\mathbf{B}$ in
spectral norm as opposed to the extensively studied and well-understood
Frobenius norm. To the best of our knowledge the only relevant related work on
Spectral Regression is by Clarkson and Woodruff [CW09] and Cohen et. al.
[CNW15]. Both these works provide oblivious sketches to reduce the dimension
of the problem, which unfortunately do not suffice for our application.
Instead of Spectral Regression in its full generality we focus on the
following special case:

Given an $n\times n$ PSD matrix $\mathbf{A}$, a rank parameter $k$, and an
accuracy parameter $\epsilon$, let $\mathbf{C}$ be a
$n\times\sqrt{nk/\epsilon}$ matrix such that it is a column PCP for
$\mathbf{A}$, satisfying the guarantees of Lemma 4.11, instantiated with
$k=k/\epsilon$, and $\epsilon=O(1)$. Let $\mathbf{Z}^{\top}$ be a
$k/\epsilon\times\sqrt{nk/\epsilon}$ matrix with orthonormal rows such that
the corresponding projection matrix $\mathbf{Z}\mathbf{Z}^{\top}$ is an
$(O(1),k/\epsilon)$-SF Projection for $\mathbf{C}$. Then, we consider the
following Spectral Regression problem:

$\min_{\mathbf{W}\in\mathbb{R}^{n\times
k/\epsilon}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|_{2}$ (4.15)

Our main technical contribution here is to obtain a new algorithm to solve
this optimization problem. We subsequently show how understanding this special
case is crucial to obtaining optimal algorithms for low rank approximation of
PSD matrices. The techniques we develop here may be of independent interest
and find applications to other problems. Formally, we prove the following:

######  Theorem 4.18.

(Approximate Spectral Regression.) Let
$\mathbf{C}\in\mathbb{R}^{n\times\sqrt{nk/\epsilon}}$ be a column PCP for
$\mathbf{A}$ satisfying the guarantees of Lemma 4.11 instantiated with
$k=k/\epsilon$ and $\epsilon=O(1)$. Let
$\mathbf{Z}\in\mathbb{R}^{\sqrt{nk/\epsilon}\times k/\epsilon}$ be an
orthonormal matrix such that $\mathbf{Z}\mathbf{Z}^{\top}$ is an
$(O(1),k/\epsilon)$-SF projection for $\mathbf{C}$. Then, Algorithm 4.2
queries $\widetilde{O}(nk/\epsilon)$ entries in $\mathbf{A}$ and with
probability $99/100$ computes $\widehat{\mathbf{W}}$ such that

$\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}\leq\widetilde{O}(1)\left(\min_{\mathbf{W}\in\mathbb{R}^{n\times
k/\epsilon}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{C}-\mathbf{C}_{k/\epsilon}\|^{2}_{F}\right)$

Further, the algorithm runs in time
$\widetilde{O}(nk/\epsilon+(k/\epsilon)^{\omega})$.

Input: A PSD Matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, integer $k$, and
$\epsilon>0$. $\mathbf{C}\in\mathbb{R}^{n\times\sqrt{nk/\epsilon}}$, a column
PCP for $\mathbf{A}$ satisfying the guarantees of Lemma 4.11 instantiated with
$k=k/\epsilon$ and $\epsilon=O(1)$.
$\mathbf{Z}\in\mathbb{R}^{\sqrt{nk/\epsilon}\times k/\epsilon}$ be an
orthonormal matrix such that $\mathbf{Z}\mathbf{Z}^{\top}$ is an
$(O(1),k/\epsilon)$-SF projection for $\mathbf{C}$. 1. Consider the Spectral
Regression problem:
$\min_{\mathbf{W}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}$ Let
$t=\sqrt{nk/\epsilon}$. For all $j\in[t]$, compute
$\tau_{j}(\mathbf{Z}^{\top})=\|\mathbf{Z}_{j,*}\|^{2}_{2}$. Let
$q=\\{q_{1},q_{2},\ldots,q_{t}\\}$ be a distribution of columns of
$\mathbf{C}$ such that for all $j\in[t]$,
$q_{j}=\min(\tau_{j}(\mathbf{Z}^{\top}),1)$. 2. Construct a sampling matrix
$\SS$ such that $\mathbf{C}\SS$ selects each column $\mathbf{C}_{*,j}$
independently with probability $q_{j}$ and scales it by $1/\sqrt{q_{j}}$.
Similarly, construct $\mathbf{Z}^{\top}\SS$. Consider the sketched
optimization problem :
$\min_{\mathbf{W}}\|\mathbf{C}\SS-\mathbf{W}\mathbf{Z}^{\top}\SS\|^{2}_{2}$ 3.
Compute
$(\mathbf{Z}^{\top}\SS)^{\dagger}=\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}$.
Let $\widehat{\mathbf{W}}=\mathbf{C}\SS(\mathbf{Z}^{\top}\SS)^{\dagger}$ be
the solution to the sketched optimization problem. Output:
$\widehat{\mathbf{W}}\in\mathbb{R}^{n\times k/\epsilon}$ such that
$\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}\leq\widetilde{O}(1)\min_{\mathbf{W}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{F}+\widetilde{O}(\epsilon/k)\|\mathbf{C}-\mathbf{C}_{k/\epsilon}\|^{2}_{F}$

We begin by characterizing the optimal solution and optimal cost for the
Spectral Regression problem. We prove a structural result that shows the
optimal solution for Spectral Regression is given by projecting $\mathbf{C}$
away from the span of $\mathbf{Z}^{\top}$. This matches the characterization
of the optimal solution to regression under the Frobenius norm, given by the
well-known normal equations. Recall, by definition of the Moore-Penrose
pseudoinverse, this projection matrix is
$(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}$.

Then, the optimal cost for Equation 4.15 is
$\|\mathbf{C}-\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}\|^{2}_{2}$
and is achieved by $\mathbf{W}^{*}=\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}$.
Intuitively, we show that any feasible $\mathbf{W}$ must incur the above cost
by analyzing
$\|y^{\top}(\mathbf{C}-\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top})\|^{2}_{2}$
for a fixed vector $y$. This enables us to exploit the geometry of Euclidean
space and instantiate $y$ as needed to relate it back to the spectral norm.

######  Lemma 4.19.

(Characterizing Opt for Spectral Regression.) Let $\mathbf{C}$ and
$\mathbf{Z}$ be matrices as defined in Theorem 4.18. Let
$\mathbf{W}^{*}=\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}=\mathbf{C}\mathbf{Z}(\mathbf{Z}^{\top}\mathbf{Z})^{-1}$,
such that $\mathbf{W}^{*}\mathbf{Z}^{\top}$ is the projection of $\mathbf{C}$
on the colspan$(\mathbf{Z})$. Let
$\mathbf{C}^{*}=\mathbf{C}-\mathbf{W}^{*}\mathbf{Z}^{\top}$ be the projection
of $\mathbf{C}$ orthogonal to colspan$(\mathbf{Z})$. Then,

$\|\mathbf{C}^{*}\|^{2}_{2}=\min_{\mathbf{W}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}$

and the corresponding minimizer is $\mathbf{W}^{*}$.

###### Proof.

Note, by definition
$\|\mathbf{C}-\mathbf{W}^{*}\mathbf{Z}^{\top}\|^{2}_{2}=\|\mathbf{C}^{*}\|^{2}_{2}$
and since $\mathbf{W}^{*}$ is feasible,

$\min_{\mathbf{W}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}\leq\|\mathbf{C}^{*}\|^{2}_{2}.$

Therefore, it suffices to show any $\mathbf{W}$ must incur cost at least
$\|\mathbf{C}^{*}\|^{2}_{2}$. By definition, we have

$\begin{split}\mathbf{C}=\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}+\mathbf{C}(\mathbf{I}-(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top})=\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}+\mathbf{C}^{*}\\\
\end{split}$

By definition of spectral norm,
$\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}\geq\|y^{\top}\mathbf{C}-y^{\top}\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}$,
for all $y$ such that $\|y\|_{2}=1$. Next, for any unit vector
$y\in\mathbb{R}^{n}$,

$\begin{split}\|y^{\top}\mathbf{C}-y^{\top}\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}&=\|y^{\top}(\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}+\mathbf{C}^{*})-y^{\top}\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}\\\
&=\|y^{\top}\mathbf{C}^{*}-y^{\top}(\mathbf{W}-\mathbf{W}^{*})\mathbf{Z}^{\top}\|^{2}_{2}\\\
&=\|y^{\top}\mathbf{C}^{*}\|^{2}_{2}+\|y^{\top}(\mathbf{W}-\mathbf{W}^{*})\mathbf{Z}^{\top}\|^{2}_{2}+2\langle
y^{\top}\mathbf{C}^{*},y^{\top}(\mathbf{W}-\mathbf{W}^{*})\mathbf{Z}^{\top}\rangle\end{split}$
(4.16)

We observe that
$\mathbf{W}\mathbf{Z}^{\top}=\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}$
is the projection of $\mathbf{C}$ on the rowspan of $\mathbf{Z}^{\top}$ and
$\mathbf{C}^{*}$ is the projection of $\mathbf{C}$ on the orthogonal
complement of rowspan of $\mathbf{Z}^{\top}$. Therefore,
$\langle\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top},\mathbf{C}^{*}\rangle=0$.
Further, for any $y$, $y^{\top}(\mathbf{W}-\mathbf{W}^{*})\mathbf{Z}^{\top}$
is in the row span of $\mathbf{Z}^{\top}$ and is thus perpendicular to
$y^{\top}\mathbf{C}^{*}$. Plugging this back in to Equation 4.16, we have

$\|y^{\top}\mathbf{C}-y^{\top}\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}=\|y^{\top}\mathbf{C}^{*}\|^{2}_{2}+\|y^{\top}(\mathbf{W}-\mathbf{W}^{*})\mathbf{Z}^{\top}\|^{2}_{2}\geq\|y^{\top}\mathbf{C}^{*}\|^{2}_{2}$
(4.17)

where the inequality follows from non-negativity of norms. Since Equation 4.17
holds for all $y$, we can pick $y$ such that
$\|y^{\top}\mathbf{C}^{*}\|^{2}_{2}=\|\mathbf{C}^{*}\|^{2}_{2}$. Therefore,
$\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}\geq\|y^{\top}\mathbf{C}-y^{\top}\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}\geq\|\mathbf{C}^{*}\|^{2}_{2}$.
This completes the proof. ∎

Next, we sketch the Spectral Regression problem from Equation 4.15 such that
we approximately preserve the spectral norm cost of all
$\mathbf{W}\in\mathbb{R}^{n\times k/\epsilon}$. A natural approach here would
be to follow the Affine Embedding idea for Frobenius norm and hope a similar
guarantee holds for spectral norm as well. However, since $\mathbf{Z}$ could
have rank as large as $k/\epsilon$ and we can no longer obtain a relative-
error $(1+\epsilon)$-approximate Affine Embedding even for Frobenius norm
without incurring a larger dependence on $\epsilon$. Instead, we relax the
notion of approximation for our sketch. We note that it suffices to construct
a sketch $\SS$ such that if

$\widehat{\mathbf{W}}=\arg\min_{\mathbf{W}}\|\mathbf{C}\SS-\mathbf{W}\mathbf{Z}^{\top}\SS\|^{2}_{2}$

then

$\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}\leq\widetilde{O}(1)\left(\min_{\mathbf{W}\in\mathbb{R}^{n\times
k/\epsilon}}\|\mathbf{C}-\mathbf{W}\mathbf{Z}^{\top}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{C}-\mathbf{C}_{k/\epsilon}\|^{2}_{F}\right)$

as stated in Theorem 4.18. Note, here we only need to weakly preserve the cost
of the optimal $\mathbf{W}$ for the sketched problem as opposed to preserving
the cost of all matrices $\mathbf{W}$. At a high level, this comes down to
analyzing the spectrum of $\|\mathbf{C}^{*}\SS\SS^{\top}\|_{2}$ We begin with
the definition of the Approximate Matrix Multiplication (AMM) guarantee and
discuss its application in approximately minimizing Spectral Regression.

######  Definition 4.20.

($(\epsilon,k)$-Spectral AMM.) Given matrices
$\mathbf{A}\in\mathbb{R}^{n\times m}$ and $\mathbf{B}\in\mathbb{R}^{m\times
d}$, a sketch $\mathbf{\Pi}\in\mathbb{R}^{m\times t}$ satisfies
$(\epsilon,k)$-Spectral AMM if with probability at least $1-\delta$,

$\|\mathbf{A}\mathbf{\Pi}\mathbf{\Pi}^{\top}\mathbf{B}-\mathbf{A}\mathbf{B}\|_{2}\leq\epsilon\sqrt{\left(\|\mathbf{A}\|^{2}_{2}+\frac{\|\mathbf{A}\|^{2}_{F}}{k}\right)\cdot\left(\|\mathbf{B}\|^{2}_{2}+\frac{\|\mathbf{B}\|^{2}_{F}}{k}\right)}$

Approximate Matrix Multiplication was introduced by Drineas et al. [DKM06]
with respect to the Frobenius norm, as opposed to the spectral norm above.
Subsequent work by Cohen et al. [CNW15] studied the spectral norm bound and
showed that any sketch $\mathbf{\Pi}$ that is an oblivious subspace embedding
(i.e., satisfies Lemma 4.3 with $\mathbf{\Pi}$ being an oblivious sketch)
implies an AMM guarantee, as long as $\mathbf{\Pi}$ has
$\Theta(k+\log(1/\delta)/\epsilon^{2})$ columns. The Spectral AMM property
combined with an $O(1)$-Subspace Embedding suffice to approximately minimize
the Spectral Regression problem :

######  Theorem 4.21.

(Theorem 3, [CNW15]) Let $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{\Pi}$ be as
defined above. If $\mathbf{\Pi}$ is an
$(\sqrt{\epsilon},\textrm{rank}(\mathbf{A}))$-Spectral AMM for
$\mathbf{U}_{A}$ and $(\mathbf{I}-\mathbf{P}_{\mathbf{A}})\mathbf{B}$, and an
$O(1)$-Subspace Embedding for $\mathbf{A}$, and
$\widehat{\mathbf{X}}=\arg\min_{\mathbf{X}}\|\mathbf{\Pi}\mathbf{A}\mathbf{X}-\mathbf{\Pi}\mathbf{B}\|^{2}_{2}$,
then with probability $99/100$,

$\|\mathbf{A}\widehat{\mathbf{X}}-\mathbf{B}\|^{2}_{2}\leq(1+\epsilon)\|\mathbf{P}_{\mathbf{A}}\mathbf{B}-\mathbf{B}\|^{2}_{2}+\frac{\epsilon}{k}\|\mathbf{P}_{\mathbf{A}}\mathbf{B}-\mathbf{B}\|^{2}_{F}$

where $\mathbf{U}_{\mathbf{A}}$ is an orthonormal basis for $\mathbf{A}$ and
$\mathbf{P}_{\mathbf{A}}$ is the projection onto the span of $\mathbf{A}$.

However, all the constructions presented for the sketch in [CEM+15] are either
oblivious sketches or require sampling proportional to both $\mathbf{A}$ and
$\mathbf{B}$. Applying an oblivious sketch $\SS$ in our problem requires
computing $\mathbf{C}\SS$ which would query
$\Omega(\textsf{nnz}(\mathbf{C}))=\Omega(n^{1.5}\sqrt{k/\epsilon})$ entries in
$\mathbf{A}$. Therefore, the main challenge here is to construct a sampling
matrix $\SS$ while reading $\widetilde{O}(nk/\epsilon)$ entries in
$\mathbf{A}$ such that $\SS$ is an $(\widetilde{O}(1),k/\epsilon)$-Spectral
AMM and an $O(1)$-Subspace Embedding. We construct $\SS$ by sampling
$\widetilde{O}(k/\epsilon)$ columns of $\mathbf{Z}^{\top}$ proportional to the
leverage scores of $\mathbf{Z}^{\top}$. While it is easy to show $\SS$ is a
Subspace Embedding, observe that our sampling probabilities are computed
without reading $\mathbf{C}$.

Proof of Theorem 4.18. As a starting point, we observe that yet again, since
$\mathbf{Z}^{\top}$ has orthonormal rows, the leverage scores are simply the
$\ell^{2}_{2}$ norms of the columns of $\mathbf{Z}^{\top}$. Therefore, one
possible approach is to construct a leverage score sampling sketch $\SS$ for
$\mathbf{C}$, by sampling columns proportional to the leverage scores of
$\mathbf{Z}^{\top}$. We note we can afford to sample at most
$\widetilde{O}(k/\epsilon)$ columns, since our algorithm queries all entries
in the resulting sketched matrix $\mathbf{C}\SS$.

Further, for reasons to be discussed later, it is crucial that we sample
columns of $\mathbf{C}$ independently, as opposed to the standard way of
sampling with replacement we have used thus far. The independent sampling
process can be described as follows: for all $j\in[\sqrt{nk/\epsilon}]$, we
sample $\mathbf{C}^{*}_{*,j}$ with probability
$\min(\|\mathbf{Z}^{\top}_{*,j}\|^{2}_{2},1)$. We use the following lemma from
[CMM17] to show that independently sampling columns satisfies some desirable
properties.

######  Lemma 4.22.

(Lemma 21, [CMM17].) Given a matrix $\mathbf{M}\in\mathbb{R}^{n\times m}$, for
all $j\in[m]$ let
$\bar{\rho}^{k}_{j}(\mathbf{M})=\Theta(\rho^{k}_{j}(\mathbf{M})$ be estimates
of the rank-$k$ column ridge-leverage scores of $\mathbf{M}$ and let
$q_{j}=\min(\bar{\rho}^{k}_{j}(\mathbf{M})\log(k/\delta)/\epsilon^{2},1)$.
Then, construct $\mathbf{M}\SS$ by selecting each column $\mathbf{M}_{*,j}$
with probability $q_{j}$ and scale it by $1/\sqrt{q_{j}}$. Then, with
probability at least $1-\delta$, $\mathbf{M}\SS$ has
$\sum_{j\in[m]}\bar{\rho}^{k}_{j}\cdot\log(k/\delta)/\epsilon^{2}$ columns and

$(1-\epsilon)\mathbf{M}\SS\SS^{\top}\mathbf{M}^{\top}-\frac{\epsilon}{k}\|\mathbf{M}-\mathbf{M}_{k}\|^{2}_{F}\mathbf{I}\preceq\mathbf{M}\mathbf{M}^{\top}\preceq(1+\epsilon)\mathbf{M}\SS\SS^{\top}\mathbf{M}^{\top}+\frac{\epsilon}{k}\|\mathbf{M}-\mathbf{M}_{k}\|^{2}_{F}\mathbf{I}$

The above lemma independently samples columns proportional to the ridge
leverage scores. In our setting, we can set the ridge parameter $\lambda=0$,
and sample according to the exact leverage scores of $\mathbf{Z}^{\top}$.
Formally, let $q=\\{q_{1},q_{2},\ldots,q_{m}\\}$ be the corresponding
distribution over columns of $\mathbf{Z}^{\top}$ such that
$q_{j}=\min(\|\mathbf{Z}^{\top}\|^{2}_{2}\log(k),1)$. Since
$\mathbf{Z}^{\top}$ has $k/\epsilon$ orthonormal rows, the leverage scores sum
up to $\textrm{rank}(\mathbf{Z}^{\top})\leq k/\epsilon$. We then use Lemma
4.22 by setting $\epsilon=1/10,\delta=0.01$ and thus with probability at least
$99/100$, $\mathbf{Z}^{\top}\SS$ has
$\sum_{j\in[\sqrt{nk/\epsilon}]}\tau_{j}(\mathbf{Z}^{\top})\log(n)=\widetilde{O}(k/\epsilon)$
rows and

$\frac{9}{10}\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z}\preceq\mathbf{Z}^{\top}\mathbf{Z}\preceq\frac{11}{10}\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z}$
(4.18)

If the guarantee in Equation 4.18 holds for a sketch $\SS$, we refer to $\SS$
as satisfying an $O(1)$-Subspace Embedding property. Observe, this is
equivalent to $\SS$ preserving all singular values of $\mathbf{Z}^{\top}$ up
to a constant.

We can now obtain a closed form solution for the Spectral Regression problem
in the sketched space. By Lemma 4.19, the optimal solution to the optimization
problem in Step 2 of Algorithm 4.2 is given by
$\widehat{\mathbf{W}}=\mathbf{C}\SS(\mathbf{Z}^{\top}\SS)^{\dagger}$. Since
$\SS$ satisfies the $O(1)$-Subspace Embedding property in Equation 4.18, it
preserves the rank of $\mathbf{Z}^{\top}$. Therefore, $\mathbf{Z}^{\top}\SS$
has full row rank and
$(\mathbf{Z}^{\top}\SS)^{\dagger}=\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}$
and thus
$\widehat{\mathbf{W}}=\mathbf{C}\SS\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}$
is the optimal solution. Next, we bound the cost of $\widehat{\mathbf{W}}$ in
the original problem. Let
$\mathbf{P}_{\mathbf{Z}^{\top}}=\mathbf{Z}^{\dagger}\mathbf{Z}^{\top}$ be the
orthogonal projection matrix onto $\mathbf{Z}^{\top}$. Using the fact that
$\|\mathbf{M}\|^{2}_{2}=\max_{\|y\|_{2}=1}\|y^{\top}\mathbf{M}\|^{2}_{2}$ and
the Pythagorean Theorem for Euclidean space we have

$\begin{split}\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}&=\|\mathbf{C}-\mathbf{C}\SS\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\|^{2}_{2}\\\
&=\max_{\|y\|^{2}_{2}=1}\|y^{\top}\mathbf{C}\mathbf{P}_{\mathbf{Z}^{\top}}-y^{\top}\mathbf{C}\SS\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\mathbf{P}_{\mathbf{Z}^{\top}}\|^{2}_{2}+\\\
&\hskip
41.91672pt\|y^{\top}\mathbf{C}(\mathbf{I}-\mathbf{P}_{\mathbf{Z}^{\top}})-y^{\top}\mathbf{C}\SS\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}(\mathbf{I}-\mathbf{P}_{\mathbf{Z}^{\top}})\|^{2}_{2}\\\
\end{split}$ (4.19)

Here, we observe
$y^{\top}\mathbf{C}\SS\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}$
is a vector in the row space of $\mathbf{Z}^{\top}$
and$(\mathbf{I}-\mathbf{P}_{\mathbf{Z}^{\top}})$ is the projection on the
orthogonal complement of rowspan$(\mathbf{Z}^{\top})$, thus this evaluates to
$0$. Since
$\mathbf{C}(\mathbf{I}-\mathbf{P}_{\mathbf{Z}^{\top}})=\mathbf{C}^{*}$, we can
upper bound
$\|y^{\top}\mathbf{C}(\mathbf{I}-\mathbf{P}_{\mathbf{Z}^{\top}})\|_{2}$ by
$\|\mathbf{C}^{*}\|^{2}_{2}$. Similarly, we can upper bound the first term by
its spectral norm. Therefore, plugging this back into Equation 4.19,

$\begin{split}\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}&\leq\|\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}-\mathbf{C}\SS\SS^{\top}\mathbf{Z}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\mathbf{Z}^{\top}\|^{2}_{2}+\|\mathbf{C}^{*}\|^{2}_{2}\\\
&=\|\left(\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z}-\mathbf{C}\SS\SS^{\top}\mathbf{Z}\right)(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\|^{2}_{2}+\|\mathbf{C}^{*}\|^{2}_{2}\\\
&\leq\|\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z}-\mathbf{C}\SS\SS^{\top}\mathbf{Z}\|^{2}_{2}\|(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\|^{2}_{2}+\|\mathbf{C}^{*}\|^{2}_{2}\\\
\end{split}$ (4.20)

where we use that $\mathbf{Z}^{\top}$ has orthonormal columns and the sub-
multiplicativity of the spectral norm. From Equation 4.18, it follows that for
all $i\in[k/\epsilon]$,
$\sigma_{i}^{2}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})=(1\pm
0.1)^{2}\sigma^{2}_{i}(\mathbf{Z}^{\top}\mathbf{Z})=(1\pm 0.1)^{2}$.
Therefore,
$\|(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})^{-1}\|^{2}_{2}=1/\sigma^{2}_{\min}(\mathbf{Z}^{\top}\SS\SS^{\top}\mathbf{Z})\leq
100/81$. Substituting this back into Equation 4.20, we have

$\begin{split}\|\mathbf{C}-\widehat{\mathbf{W}}\mathbf{Z}^{\top}\|^{2}_{2}&\leq
O(1)\|(\mathbf{C}(\mathbf{Z}^{\top})^{\dagger}\mathbf{Z}^{\top}-\mathbf{C})\SS\SS^{\top}\mathbf{Z}\|^{2}_{2}+\|\mathbf{C}^{*}\|^{2}_{2}\\\
&\leq
O(1)\|\mathbf{C}^{*}\SS\SS^{\top}\mathbf{Z}\|^{2}_{2}+\|\mathbf{C}^{*}\|^{2}_{2}\\\
\end{split}$ (4.21)

where the last inequality follows from the definition of $\mathbf{C}^{*}$. In
order to bound the cost above, we focus on analyzing
$\|\mathbf{C}^{*}\SS\SS^{\top}\mathbf{Z}\|^{2}_{2}$. Since we want to compare
$\|\mathbf{C}^{*}\SS\SS^{\top}\mathbf{Z}\|^{2}_{2}$ to
$\|\mathbf{C}^{*}\mathbf{Z}\|^{2}_{2}$, a natural way to proceed would be to
interpret this term as an instance of Approximate Matrix Product. Therefore,
we next show that the leverage score sampling matrix $\SS$ satisfies the
Spectral AMM property for $\mathbf{C}^{*}$ and $\mathbf{Z}^{\top}$. Here, we
want to analyze how sampling columns of $\mathbf{C}^{*}$ proportional to the
leverage scores of $\mathbf{Z}^{\top}$ affects the spectrum of
$\mathbf{C}^{*}$. An important tool in this analysis is the following result
by Rudelson and Vershynin on how the spectral norm of a matrix degrades when
we sample a uniformly random subset of rows of a matrix:

######  Theorem 4.23.

( Theorem 1.8 in [RV07]) Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$,
let $\mathcal{Q}$ be a uniformly random subset of $[n]$ s.t.
${\bf\mathbb{E}}\left[\mathcal{Q}\right]=q$. Let $\mathbf{A}_{|\mathcal{Q}}$
denote the submatrix restricted to the rows indexed by $\mathcal{Q}$. Then,

${\bf\mathbb{E}}\left[\big{\|}\mathbf{A}_{|\mathcal{Q}}\big{\|}_{2}\right]=O\left(\sqrt{\frac{q}{n}}\|\mathbf{A}\|_{2}+\sqrt{\log(q)}\|\mathbf{A}\|_{(n/q)}\right)$

where $\|\mathbf{A}\|_{(n/q)}$ is the average of the largest $n/q$