Algorithm

# Low-Rank Approximation with $1/\epsilon^{1/3}$ Matrix-Vector Products

Ainesh Bakshi
<EMAIL_ADDRESS>
CMU  Kenneth L. Clarkson
<EMAIL_ADDRESS>
IBM  David P. Woodruff
<EMAIL_ADDRESS>
CMU

###### Abstract

We study iterative methods based on Krylov subspaces for low-rank
approximation under any Schatten-$p$ norm. Here, given access to a matrix
$\mathbf{A}$ through matrix-vector products, an accuracy parameter $\epsilon$,
and a target rank $k$, the goal is to find a rank-$k$ matrix $\mathbf{Z}$ with
orthonormal columns such that
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq(1+\epsilon)\min_{\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{k}}\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top}\right)\right\|_{\mathcal{S}_{p}}$,
where $\|\mathbf{M}\|_{\mathcal{S}_{p}}$ denotes the $\ell_{p}$ norm of the
the singular values of $\mathbf{M}$. For the special cases of $p=2$ (Frobenius
norm) and $p=\infty$ (Spectral norm), Musco and Musco (NeurIPS 2015) obtained
an algorithm based on Krylov methods that uses $\tilde{O}(k/\sqrt{\epsilon})$
matrix-vector products, improving on the naïve $\tilde{O}(k/\epsilon)$
dependence obtainable by the power method, where $\tilde{O}(\cdot)$ suppresses
poly$(\log(dk/\epsilon))$ factors.

Our main result is an algorithm that uses only
$\tilde{O}(kp^{1/6}/\epsilon^{1/3})$ matrix-vector products, and works for
_all_ , not necessarily constant, $p\geq 1$. For $p=2$ our bound improves the
previous $\tilde{O}(k/\epsilon^{1/2})$ bound to $\tilde{O}(k/\epsilon^{1/3})$.
Since the Schatten-$p$ and Schatten-$\infty$ norms of any matrix are the same
up to a $1+\epsilon$ factor when $p\geq(\log d)/\epsilon$, our bound recovers
the result of Musco and Musco for $p=\infty$. Further, we prove a matrix-
vector query lower bound of $\Omega(1/\epsilon^{1/3})$ for _any_ fixed
constant $p\geq 1$, showing that surprisingly
$\tilde{\Theta}(1/\epsilon^{1/3})$ is the optimal complexity for constant $k$.

To obtain our results, we introduce several new techniques, including
optimizing over _multiple Krylov subspaces_ simultaneously, and _pinching
inequalities_ for partitioned operators. Our lower bound for $p\in[1,2]$ uses
the Araki-Lieb-Thirring trace inequality, whereas for $p>2$, we appeal to a
norm-compression inequality for aligned partitioned operators. As our
algorithms only require matrix-vector product access, they can be applied in
settings where alternative techniques such as sketching cannot, e.g., to
covariance matrices, Hessians defined implicitly by a neural network, and
arbitrary polynomials of a matrix.

##  1 Introduction

Iterative methods, and in particular Krylov subspace methods, are ubiquitous
in scientific computing. Algorithms such as power iteration, Golub-Kahan
Bidiagonalization, Arnoldi iteration, and the Lanczos iteration, are used in
basic subroutines for matrix inversion, solving linear systems, linear
programming, low-rank approximation, and numerous other fundamental linear
algebra primitives [Saa81, LS13]. A common technique in the analysis of Krylov
methods is the use of Chebyshev polynomials, which can be applied to the
singular values of a matrix to implement an approximate interval or step
function [MH02, Riv20]. Further, Chebyshev polynomials reduce the degree
required to accurately approximate such functions, leading to significantly
fewer iterations and faster running time. In this paper we investigate the
power of Krylov methods for low-rank approximation in the matrix-vector
product model.

The Matrix-Vector Product Model. In this model, there is an underlying matrix
$\mathbf{A}$, which is often implicit, and for which the only access to
$\mathbf{A}$ is via matrix-vector products. Namely, the algorithm chooses a
query vector $v^{1}$, obtains the product $\mathbf{A}\cdot v^{1}$, chooses the
next query vector $v^{2}$, which is any randomized function of $v^{1}$ and
$\mathbf{A}\cdot v^{1}$, then receives $\mathbf{A}\cdot v^{2}$, and so on. If
$\mathbf{A}$ is a non-symmetric matrix, we assume access to products of the
form $\mathbf{A}^{\top}v$ as well. We refer to the minimal number $q$ of
queries needed by the algorithm to solve a problem with constant probability
as the query complexity. We note that upper bounds on the query complexity
immediately translate to running time bounds for the RAM model, when
$\mathbf{A}$ is explicit, since a matrix-vector product can be implemented in
$\operatorname{\mathtt{nnz}}(\mathbf{A})$ time, i.e., the number of non-zero
entries in the matrix. Since this model captures a large family of iterative
methods, it is natural to ask whether Krylov subspace based methods yield
optimal algorithms, where the complexity measure of interest is the number of
matrix-vector products.

This model and related vector-matrix-vector query models were formalized for a
number of problems in [SWYZ19, RWZ20], though the model is standard for
measuring efficiency in scientific computing and numerical linear algebra,
see, e.g., [BFG96]; in that literature, methods that use only matrix-vector
products are called _matrix-free_. Subsequently, for the problem of estimating
the top eigenvector, nearly tight bounds were obtained in [SAR18, BHSW20].
Also, for the problem of estimating the trace of a positive semidefinite
matrix, tight bounds were obtained in [MMMW21] (see, also [WWZ14], where tight
bounds were shown in the related vector-matrix-vector query model). For
recovering a planted clique from a random graph, upper and lower bounds were
obtained in [RWYZ21]. In the non-adaptive setting, where $v^{1},\ldots,v^{q}$,
are chosen before making any queries to $\mathbf{A}$, this is equivalent to
the sketching model, which is thoroughly studied on its own (see, e.g.,
[Nel11, Woo14]), and in the context of data streams [Mut05, LNW14b].

Why is the matrix $\mathbf{A}$ implicit? A small query complexity $q$ leads to
an algorithm running in time $\mathcal{O}\left(T(\mathbf{A})\cdot
q+P(n,d,q)\right)$, where $T(\mathbf{A})$ is the time to multiply the $n\times
d$ matrix $\mathbf{A}$ by an arbitrary vector, and $P(n,d,q)$ is the time
needed to form the queries and process the query responses, which is typically
small. When the matrix $\mathbf{A}$ is given as a list of
$\operatorname{\mathtt{nnz}}(\mathbf{A})$ non-zero entries, then
$T(\mathbf{A})\leq\operatorname{\mathtt{nnz}}(\mathbf{A})$. However, in many
problems $\mathbf{A}$ is not given explicitly, and it is too expensive to
write $\mathbf{A}$ down. Indeed, one may be given $\mathbf{A}$ but want to
compute a low-rank approximation to the “covariance” (Gram) matrix
$\mathbf{A}^{\top}\mathbf{A}$, and computing $\mathbf{A}^{\top}\mathbf{A}$ is
too slow [MW17a]. More generally, one may be given
$\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ and a function
$f:\mathbb{R}\to\mathbb{R}$, and want to compute matrix-vector products with
the generalized matrix function
$f(\mathbf{A})=\mathbf{U}f(\mathbf{\Sigma})\mathbf{V}^{\top}$, where
$\mathbf{U}$ has orthonormal columns, $\mathbf{V}^{\top}$ has orthonormal
rows, $\mathbf{\Sigma}$ is a diagonal matrix, and $f$ is applied entry-wise to
each entry on the diagonal.

The covariance matrix corresponds to $f(x)=x^{2}$, and other common functions
$f$ include the matrix exponential $f(x)=e^{x}$ and low-degree polynomials.
For instance, when $\mathbf{A}$ is the adjacency matrix of an undirected
graph, $f(x)=x^{3}/6$ is used to count the number of triangles [Tso08, Avr10].
Yet another example is when $\mathbf{A}$ is the Hessian $\mathbf{H}$ of a
neural network with a huge number of parameters, for which it is often
impossible to compute or store the entire Hessian [GKX19]. Typically
$\mathbf{H}\cdot v$, for any chosen vector $v$, is computed using
Pearlmutter’s trick [Pea94]. However, even with Pearlmutter’s trick and
distributed computation on modern GPUs, it takes 20 hours to compute the
eigendensity of a single Hessian $\mathbf{H}$ with respect to the cross-
entropy loss on the CIFAR-10 dataset from a set of fixed weights for ResNet-18
[KH+09], which has approximately 11 million parameters [HZRS16, GKX19]. This
time is directly proportional to the number of matrix-vector products, and
therefore minimizing this quantity is crucial.

Algorithms and Lower Bounds for Low-Rank Approximation. The low-rank
approximation problem is well studied in numerical linear algebra, with
countless applications to clustering, data mining, principal component
analysis, recommendation systems, and many more. (For surveys on low-rank
approximation, see the monographs [KV09, Mah11, Woo14] and references
therein.) In this problem, given an implicit $n\times d$ matrix $\mathbf{A}$,
the goal is to output a matrix $\mathbf{Z}\in\mathbb{R}^{d\times k}$ with
orthonormal columns such that

$\displaystyle\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{X}\leq\left(1+\epsilon\right)\min_{\mathbf{U}:\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{k}}\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top}\right)\right\|_{X},$
(1.1)

where $\|\cdot\|_{X}$ denotes some norm. Note that given $\mathbf{Z}$, one can
compute $\mathbf{A}\mathbf{Z}$ with an additional $k$ queries, which will be
negligible, and then $(\mathbf{A}\mathbf{Z})\cdot\mathbf{Z}^{\top}$ is a
rank-$k$ matrix written in factored form, i.e., as the product of an $n\times
k$ matrix and a $k\times d$ matrix. Among other things, low-rank approximation
provides (1) a compression of $\mathbf{A}$ from $nd$ parameters to $(n+d)k$
parameters, (2) faster matrix-vector products, since
$\mathbf{A}\mathbf{Z}\cdot\mathbf{Z}^{\top}\cdot y$ can be computed in
$O((n+d)k)$ time for an arbitrary vector $y$, as opposed to the $O(nd)$ time
needed to compute $\mathbf{A}\cdot y$, and (3) de-noising, as often matrices
$\mathbf{A}$ are close to low-rank (e.g., they are the product of latent
factors) but only high rank due to noise.

Despite its tremendous importance, the optimal matrix-vector product
complexity of low-rank approximation is unknown for any commonly used norm.
The best known upper bound is due to Musco and Musco [MM15], who achieve
$\tilde{\mathcal{O}}(k/\epsilon^{1/2})$ queries111We let
$\tilde{\mathcal{O}}(f)=f\cdot\textrm{poly}(\log(dk/\epsilon))$. for both the
case when $\|\cdot\|_{X}$ is the commonly studied Frobenius norm
$\|\mathbf{B}\|_{F}=\left(\sum_{i,j}\mathbf{B}_{i,j}^{2}\right)^{1/2}$ as well
as when $\|\cdot\|_{X}$ is the Spectral (operator) norm
$\|\mathbf{B}\|_{2}=\sup_{\|y\|_{2}=1}\|\mathbf{B}y\|_{2}$.

On the lower bound front, there is a trivial lower bound of $k$, since
$\mathbf{A}$ may be full rank and achieving (1.1) requires $k$ matrix-vector
products since one must reconstruct the column span of $\mathbf{A}$ exactly.
However, no lower bounds in terms of the approximation factor $\epsilon$ were
known. We note that Simchowitz, Alaoui and Recht [SAR18] prove lower bounds
for approximating the top $r$ eigenvalues of a symmetric matrix; however these
guarantees are incomparable to those that follow from a low-rank
approximation, even when the norm $\|\cdot\|_{X}$ is the operator norm (see
Appendix A for a brief discussion).

Relationship to the Sketching Literature. Low-rank approximation has been
extensively studied in the sketching literature which, when $\mathbf{A}$ is
given explicitly, can achieve
$\mathcal{O}\left(\operatorname{\mathtt{nnz}}(\mathbf{A})\right)$ time both
for the Frobenius norm [CW13, MM13a, NN13], as well as for Schatten-$p$ norms
[LW20]. However, these works require reading all of the entries in
$\mathbf{A}$, and thus do not apply to any of the settings mentioned above.
Further, the matrix-vector query model is especially important for problems
such as trace estimation, where a low-rank approximation is used to first
reduce the variance [MMMW21]. As trace estimation is often applied to implicit
matrices, e.g., in computing Stochastic Lanczos Quadrature (SLQ) for Hessian
eigendensity estimation [GKX19], in studying the effects of batch
normalization and residual connections in neural networks [YGKM20], and in
computing a disentanglement regularizer for deep generative models [PPZ+20],
sketching algorithms for low-rank approximation often do not apply.

Another important application is low-rank approximation of covariance matrices
[MW17a], for which the covariance matrix is not given explicitly. Here, we
have a data matrix $\mathbf{A}$ and we want a low-rank approximation for
$\mathbf{A}\mathbf{A}^{\top}$. Even when $\SS$ is a sparse sketching matrix,
the matrix $\SS\mathbf{A}$ is no longer sparse, and one needs to multiply
$\SS\mathbf{A}$ by $\mathbf{A}^{\top}$ to obtain a sketch of
$\SS\mathbf{A}\mathbf{A}^{\top}$, which is a dense matrix-matrix
multiplication. Moreover, when viewed in the matrix-vector product model,
sketching algorithms obtain provably worse query complexity than existing
iterative algorithms (see Table 1.1 for a comparison). Further, as modern GPUs
often do not exploit sparsity, even when the matrix $\mathbf{A}$ is given, a
GPU may not be able to take advantage of sparse queries, which means the total
time taken is proportional to the number of matrix-vector products.

Motivating Schatten-$p$ Norms. The Schatten norms for $1\leq p<2$ are more
robust than the Frobenius norm, as they dampen the effect of large singular
values. In particular, the Schatten-$1$ norm, also known as the nuclear norm,
has been widely used for robust PCA [XCS10, CLMW11, YPCC16] as well as a
convex relaxation of matrix rank in matrix completion [CR09, CP10], low-
dimensional Euclidean embeddings [RFP10, TDSL00, RS00], image denoising
[GZZF14, GXM+17] and tensor completion [YZ16]. In contrast, for $p>2$,
Schatten norms are more sensitive to large singular values and provide an
approximation to the operator norm. In particular, for a rank $r$ matrix, it
is easy to see that setting $p=\log(r)/\eta$ yields a $(1+\eta)$-approximation
to the operator norm (i.e., $p=\infty$). While the Block Krylov algorithm of
Musco and Musco [MM15] implies a matrix-vector query upper bound of
$\tilde{\mathcal{O}}\left(k/\epsilon^{1/2}\right)$ for Schatten-$\infty$ low-
rank approximation, the exact complexity of this problem remains an
outstanding open problem. When $p>2$, we can interpolate between Frobenius and
operator norm, and setting $p$ to be a large fixed constant can be a proxy for
Schatten-$\infty$ low-rank approximation, with significantly fewer matrix-
vector products (see Theorem 5.1).

Our Central Question. The main question of our work is:

What is the matrix-vector product complexity of low-rank approximation for the
Frobenius norm, and more generally, for other matrix norms?

###  1.1 Our Results

Problem | Frobenius | Schatten-$p$, $p\in[1,2)$ | Schatten-$p$, $p>2$
---|---|---|---
Sketching [CW09, LW20] | $\Theta(k/\epsilon)$ | $\Omega(k^{2/p}/\epsilon^{4/p+1})$ | $\Omega(\min(n,d)^{1-2/p})$
Block Krylov [MM15] | $\tilde{\mathcal{O}}(k/\epsilon^{1/2})$ | N/A | N/A
Our Upper Bound | $\tilde{\mathcal{O}}(k/\epsilon^{1/3})$ | $\tilde{\mathcal{O}}(k/\epsilon^{1/3})$ | $\tilde{\mathcal{O}}(kp^{1/6}/\epsilon^{1/3})$
Our Lower Bound | $\Omega(1/\epsilon^{1/3})$ | $\Omega(1/\epsilon^{1/3})$ | $\Omega(1/\epsilon^{1/3})$

Figure 1.1: Prior Upper and Lower Bounds on the Matrix Vector Product
Complexity for Frobenius and Schatten-$p$ low-rank Approximation. The
poly$(k/\epsilon)$ factors in prior sketching work for Schatten-$p$ are not
explicit, but we have computed lower bounds on them to illustrate our
improvements. Our bounds are optimal, up to logarithmic factors, for constant
$k$. For $p>\log(d)/\epsilon$, spectral low-rank approximation [MM15] implies
an $\tilde{\mathcal{O}}\left(k/\sqrt{\epsilon}\right)$ upper bound.

We begin by stating our results for Frobenius and more generally, Schatten-$p$
norm low-rank approximation for any $p\geq 1$; see Table 1.1 for a summary.

######  Theorem 1.1 (Query Upper Bound, informal Theorem 5.1).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, a target rank $k\in[d]$,
an accuracy parameter $\epsilon\in(0,1)$ and any (not necessarily constant)
$p\in[1,\mathcal{O}\left(\log(d)/\epsilon\right)]$, there exists an algorithm
that uses $\tilde{\mathcal{O}}\left(kp^{1/6}/\epsilon^{1/3}\right)$ matrix-
vector products and outputs a $d\times k$ matrix $\mathbf{Z}$ with orthonormal
columns such that with probability at least $99/100$,

$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq(1+\epsilon)\min_{\mathbf{U}:\hskip
3.61371pt\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{k}}\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top}\right)\
\right\|_{\mathcal{S}_{p}}.$

When $p\geq\log(d)/\epsilon$, we get
$\tilde{\mathcal{O}}\left(k/\sqrt{\epsilon}\right)$ matrix-vector products.

We note that for Frobenius norm low-rank approximation (Schatten $p$ for
$p=2$), we improve the prior matrix-vector product bound of
$\tilde{\mathcal{O}}(k/\epsilon^{1/2})$ by Musco and Musco [MM15] to
$\tilde{\mathcal{O}}(k/\epsilon^{1/3})$. For Schatten-$p$ low-rank
approximation for $p\in[1,2)$, we improve work of Li and Woodruff [LW20] who
require query complexity at least $\Omega(k^{2/p}/\epsilon^{4/p+1})$, which is
a polynomial factor worse in both $k$ and $1/\epsilon$ than our
$\tilde{\mathcal{O}}(k/\epsilon^{1/3})$ bound.

For $p>2$, [LW20] obtain a query complexity of $\Omega(\min(n,d)^{1-2/p})$. We
drastically improve this to $\tilde{\mathcal{O}}(k/\epsilon^{1/3})$, which
does not depend on $d$ or $n$ at all. Setting $p=\log(d)/\epsilon$ suffices to
obtain a $(1+\epsilon)$-approximation to the spectral norm ($p=\infty)$, and
we obtain an $\tilde{\mathcal{O}}\left(k/\sqrt{\epsilon}\right)$ query
algorithm, matching the best known bounds for spectral low-rank approximation
[MM15]. When $p>\log(d)/\epsilon$, we can simply run Block Krylov for
$p=\infty$.

######  Remark 1.2 (Comments on the RAM Model).

Although our focus is on minimizing the number of matrix-vector products,
which is the key resource in the applications described above, our bounds also
improve the running time of low-rank approximation algorithms when the matrix
$\mathbf{A}$ has a small number of non-zero entries and is explicitly given.
For simplicity, we state our bounds and those of previous work without using
algorithms for fast matrix multiplication; similar improvements hold when
using such algorithms. When $\operatorname{\mathtt{nnz}}(\mathbf{A})=O(n)$,
for Frobenius norm low-rank approximation, work in the sketching literature,
and in particular [ACW17] (building off of [CW13, NN13, Coh16]), achieves
$O(nk^{2}/\epsilon)$ time. In contrast, in this setting our runtime is
$\tilde{\mathcal{O}}(nk^{2}/\epsilon^{2/3})$. Similarly, for Schatten-$p$ low-
rank approximation for $p\in[1,2)$, the previous best [LW20] requires
$\tilde{\Omega}(nk^{4/p}/\epsilon^{(8/p-2)})$ time, while for $p>2$ [LW20]
requires $\tilde{\Omega}(nd^{2(1-2/p)}(k/\epsilon)^{4/p})$ time. In both cases
our runtime is only $\tilde{\mathcal{O}}(nk^{2}p^{1/3}/\epsilon^{2/3})$. We
obtain analogous improvements when the sparsity
$\operatorname{\mathtt{nnz}}(\mathbf{A})$ is allowed to be $n(k/\epsilon)^{C}$
for a small constant $C>0$.

Next, we state our lower bounds on the matrix-vector query complexity of
Schatten-$p$ low-rank approximation.

######  Theorem 1.3 (Query Lower Bound for constant $p$, informal Theorem 6.1
and Theorem 6.4 ).

Given $\varepsilon>0$, and a fixed constant $p\geq 1$, there exists a
distribution $\mathcal{D}$ over $n\times n$ matrices such that for
$\mathbf{A}\sim\mathcal{D}$, any algorithm that with at least constant
probability outputs a unit vector $v$ such that
$\left\|\mathbf{A}\left(\mathbf{I}-vv^{\top}\right)\right\|^{p}_{\mathcal{S}_{p}}\leq(1+\varepsilon)\min_{\left\|u\right\|_{2}=1}\left\|\mathbf{A}\left(\mathbf{I}-uu^{\top}\right)\right\|^{p}_{\mathcal{S}_{p}}$
must perform $\Omega(1/\varepsilon^{1/3})$ matrix-vector queries to
$\mathbf{A}$.

######  Remark 1.4.

We note that this is the first lower bound as a function of $\epsilon$ for
this problem, even for the well-studied case of $p=2$, achieving an
$\Omega(1/\epsilon^{1/3})$ bound, which is tight for any constant $k$,
simultaneously for all constant $p\geq 1$.

######  Remark 1.5.

Braverman, Hazan, Simchowitz and Woodworth [BHSW20] and Simchowitz, Alaoui and
Recht [SAR18] establish eigenvalue estimation lower bounds that we use in our
arguments, but their results do not directly imply low-rank approximation
lower bounds for any matrix norm that we are aware of, including spectral low-
rank approximation, i.e., $p=\infty$ (see Appendix A).

#### Matrix Polynomials and Streaming Algorithms.

Since our algorithms are based on iterative methods, they generalize naturally
to low-rank approximations of matrices of the form
$\mathbf{A}\left(\mathbf{A}^{\top}\mathbf{A}\right)^{\ell}$ and
$\left(\mathbf{A}^{\top}\mathbf{A}\right)^{\ell}$ for any integer $\ell$,
given $\mathbf{A}$ as input. We defer the details to Appendix B.

Since we work in the matrix-vector model, our algorithms naturally extend to
the multi-pass turnstile streaming setting. Notably, for $p>2$, with
$\mathcal{O}\left(\log(d/\epsilon)p^{1/6}/\epsilon^{1/3}\right)$ passes we are
able to improve the
$\tilde{\mathcal{O}}\left(n\left(\frac{kn^{1-2/p}}{\epsilon^{2}}+\frac{k^{2/p}+n^{1-2/p}}{\epsilon^{2+2/p}}\right)\right)$
memory bound of [LW20] to $\tilde{\mathcal{O}}\left(nk/\epsilon^{1/3}\right)$.
We defer the details to Appendix C.

###  1.2 Open Questions

We note that our lower bounds are tight only when the target rank $k$ and
Schatten norm $p$ are fixed constants. In particular, it is open to obtain
matrix-vector lower bounds that grow as a function of $k$, $p$ and
$1/\epsilon$. For the important special case of Spectral low-rank
approximation ($p=\infty$), it is open to obtain any lower bound that grows as
a function of $1/\epsilon$, even when the target rank $k=1$ (see Appendix A
for more details). We also note that improving our upper bound to even
$p^{1/6-o(1)}$ would imply a faster algorithm for Spectral low-rank
approximation, addressing the main open question in [Woo14].

##  2 Technical Overview

For our technical overview, we drop polylogarithmic factors appearing in the
analysis and assume the input $\mathbf{A}$ is a symmetric $n\times n$ matrix
(we handle arbitrary $n\times d$ matrices in Section 5).

###  2.1 Algorithms for Low-Rank Approximation

We first describe our algorithm for the special case of rank-$1$ approximation
in the Frobenius norm, i.e., $p=2$. Our algorithm is inspired by the Block
Krylov algorithm of Musco and Musco [MM15]. Briefly, their algorithm begins
with a random starting vector $g$ (block size is $1$) and computes the Krylov
subspace $\mathbf{K}=[\mathbf{A}g;\mathbf{A}^{2}g;\ldots;\mathbf{A}^{q}g]$,
for $q=\mathcal{O}\left(1/\epsilon^{1/2}\right)$. Next, their algorithm
computes an orthonormal basis for the column span of $\mathbf{K}$, denoted by
a matrix $\mathbf{Q}$, and outputs the top singular vector of
$\mathbf{Q}^{\top}\mathbf{A}^{2}\mathbf{Q}$, denoted by $z$ (see Algorithm 5.6
for a formal description). It follows from Theorem 1, guarantee (1) in [MM15]
that

$\left\|\mathbf{A}\left(\mathbf{I}-zz^{\top}\right)\right\|_{F}^{2}\leq\left(1+\epsilon\right)\min_{\left\|u\right\|_{2}=1}\left\|\mathbf{A}\left(\mathbf{I}-uu^{\top}\right)\right\|_{F}^{2},$
(2.1)

and it is easy to see that this algorithm requires
$\Theta\left(1/\epsilon^{1/2}\right)$ matrix-vector products. A naïve analysis
requires an $\mathcal{O}\left(1/\epsilon\right)$-degree polynomial in the
matrix $\mathbf{A}$ to obtain (2.1), while [MM15] use Chebyshev polynomials to
approximate the threshold function between first and second singular value,
and save a quadratic factor in the degree. The guarantee in (2.1) then follows
from observing that the best vector in the Krylov subspace is at least as good
as the one that exists using Chebyshev polynomial approximation.

######  Algorithm 2.1 (Algorithm Sketch for Frobenius rank-$1$ LRA ).

Input: An $n\times n$ symmetric matrix $\mathbf{A}$, accuracy parameter
$0<\varepsilon<1$. 1. Run Block Krylov for
$\mathcal{O}\left(1/\epsilon^{1/3}\right)$ iterations with a random starting
vector $g$. Let $z_{1}$ be the resulting output. 2. Run Block Krylov for
$\mathcal{O}\left(\log(n/\epsilon)\right)$ iterations, but initialize with an
$n\times b$ random matrix $\mathbf{G}$, where
$b=\mathcal{O}\left(1/\epsilon^{1/3}\right)$. Let $z_{2}$ be the resulting
output. Output:
$z=\arg\max_{z_{1},z_{2}}\left(\left\|\mathbf{A}z_{1}\right\|^{2}_{2},\left\|\mathbf{A}z_{2}\right\|^{2}_{2}\right)$.

Our starting point is the observation that while we require degree
$\Theta\left(1/\epsilon^{1/2}\right)$ to separate the first and second
singular values, if any subsequent singular value is sufficiently separated
from $\sigma_{1}$, a significantly smaller degree polynomial suffices. In the
context of Krylov methods, this translates to the intuition that starting with
a matrix $\mathbf{G}$ with $b$ columns (block size is $b$) should result in
fewer iterations to find some vector in the top $b$ subspace of $\mathbf{A}$.
On the other hand, if no such singular value exists, the norm of the tail must
be large and we can get away with a less accurate solution. We show that we
can indeed exploit this trade-off by running Block Krylov on two different
scales in parallel and then combine the solution. In particular, we use
Algorithm 2.1.

Algorithm 2.1 captures the extreme points of the trade-off between the size of
the starting matrix and the number of iterations, such that the total number
of matrix-vector products is at most
$\tilde{\mathcal{O}}({1/\epsilon^{1/3}})$. Further, we can compute the squared
Euclidean norms of $\mathbf{A}z_{1}$ and $\mathbf{A}z_{2}$ with an additional
matrix-vector product, and it remains to analyze the Frobenius cost of
projecting $\mathbf{A}$ on the subspace $\mathbf{I}-zz^{\top}$, where $z$ is
the unit vector output by Algorithm 2.1.

Using gap-independent guarantees for Block Krylov (see Lemma 5.2 for a formal
statement), it follows that with $\mathcal{O}\left(1/\epsilon^{1/3}\right)$
iterations, we have

$\left\|\mathbf{A}z_{1}\right\|_{2}^{2}\geq\sigma_{1}^{2}(\mathbf{A})-\epsilon^{2/3}\sigma_{2}^{2}(\mathbf{A}).$
(2.2)

In contrast, using gap-dependent guarantees (see Lemma 5.3) for Block Krylov
initialized with block size $b$, it follows that for any $\gamma>0$, running
$q=\log(1/\gamma)\cdot\sqrt{\sigma_{1}(\mathbf{A})/\left(\sigma_{1}(\mathbf{A})-\sigma_{b}(\mathbf{A})\right)}$
iterations results in $z_{2}$ such that

$\left\|\mathbf{A}z_{2}\right\|_{2}^{2}\geq\sigma_{1}^{2}(\mathbf{A})-\gamma\sigma_{2}^{2}(\mathbf{A}).$
(2.3)

If $\sigma_{b}(\mathbf{A})\leq\sigma_{1}(\mathbf{A})/2$, we can set
$\gamma=\epsilon/n$ in Equation (2.3) to obtain a highly accurate solution.
Further, regardless of the input instance, Step 3 in Algorithm 2.1 ensures
that we get the best of both guarantees, (2.2) and (2.3). Then, observing that
$\mathbf{I}-zz^{\top}$ is an orthogonal projection matrix (see Definition 4.1)
and using the Pythagorean Theorem for Euclidean space we have:

$\left\|\mathbf{A}\left(\mathbf{I}-zz^{\top}\right)\right\|_{F}^{2}=\left\|\mathbf{A}\right\|_{F}^{2}-\left\|\mathbf{A}zz^{\top}\right\|_{F}^{2}=\left\|\mathbf{A}\right\|_{F}^{2}-\left\|\mathbf{A}z\right\|_{2}^{2},$
(2.4)

where the second inequality follows from unitary invariance (see Fact 4.8) of
the Frobenius norm and that the squared Frobenius norm of a rank-$1$ matrix
$\mathbf{A}z$ (vector) is equal to its squared Euclidean norm. If it happens
that $\sigma_{2}(\mathbf{A})\leq\sigma_{1}(\mathbf{A})/2$, i.e., a constant
gap exists between the first two singular values, then since guarantee (2.3)
implies that
$\left\|\mathbf{A}z\right\|_{2}^{2}\geq\sigma_{1}^{2}(\mathbf{A})-(\epsilon/n)\sigma_{2}^{2}(\mathbf{A})$,
we can plug this into (2.4) to yield a $(1+\epsilon/n)$-approximate solution.
Hence, we focus on instances where
$\sigma_{2}(\mathbf{A})>\sigma_{1}(\mathbf{A})/2$.

Consider the case where the Frobenius norm of the tail is large, i.e.,
$\|\mathbf{A}-\mathbf{A}_{1}\|_{F}^{2}\geq\sigma_{2}^{2}(\mathbf{A})/\epsilon^{1/3}$,
where $\mathbf{A}_{1}$ is the best rank-$1$ approximation to $\mathbf{A}$.
Then we only require an $\epsilon^{2/3}$-approximate solution (plugging
guarantee (2.2) into (2.4) ) since

$\left\|\mathbf{A}\left(\mathbf{I}-z_{1}z^{\top}_{1}\right)\right\|_{F}^{2}\leq\left\|\mathbf{A}\right\|_{F}^{2}-\sigma_{1}^{2}(\mathbf{A})+\epsilon^{2/3}\sigma_{2}^{2}(\mathbf{A})\leq\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{F}^{2}+\epsilon\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{F}^{2}.$
(2.5)

Otherwise,
$\sum_{i=2}^{n}\sigma_{i}^{2}(\mathbf{A})<\sigma_{2}^{2}(\mathbf{A})/\epsilon^{1/3}$,
which implies that there is a constant gap between the second and $b$-th
singular values, where $b=\mathcal{O}\left(1/\epsilon^{1/3}\right)$. To see
this, observe if $\sigma_{b}(\mathbf{A})>\sigma_{2}(\mathbf{A})/4$, then
$\sum_{i=2}^{n}\sigma_{i}^{2}(\mathbf{A})\geq\sum_{i=2}^{b}\sigma_{i}^{2}(\mathbf{A})\geq
b\sigma_{2}^{2}(\mathbf{A})/4$, which is a contradiction when
$b>10/\epsilon^{1/3}$, and thus
$\sigma_{b}(\mathbf{A})\leq\sigma_{2}(\mathbf{A})/4<\sigma_{1}/2$. Now we can
apply guarantee (2.3) with $q=\mathcal{O}\left(\log(n/\epsilon)\right)$ and
conclude
$\left\|\mathbf{A}z\right\|_{2}^{2}\geq\sigma_{1}^{2}(\mathbf{A})-(\epsilon/n)\sigma_{2}^{2}(\mathbf{A})$,
yielding a highly accurate solution yet again. Overall, this suffices to
obtain a $(1+\epsilon)$-approximate solution with
$\tilde{\mathcal{O}}({1/\epsilon^{1/3}})$ matrix-vector queries.

#### Challenges in generalizing to Schatten $p\neq 2$ and rank $k>1$.

The outline above crucially relies on the norm of interest being Frobenius. In
particular, we use the Pythagorean Theorem to analyze the cost of the
candidate solution in Equation (2.4); however, the Pythagorean Theorem does
not hold for non-Euclidean spaces. Therefore, a priori, it is unclear how to
analyze the Schatten-$p$ norm of a candidate rank-$1$ approximation. A proxy
for the Pythagorean Theorem that holds for Schatten-$p$ norms is Mahler’s
operator inequality (see Fact 4.11), which is in the right direction but holds
only for $p\geq 2$, whereas we would like to handle all $p\geq 1$. Separately,
for $p>2$, the case where the tail is small corresponds to
$\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{\mathcal{S}_{p}}^{p}\leq\sigma_{2}^{p}\left(\mathbf{A}\right)/\epsilon^{1/3}$.
Therefore, naïvely extending the above argument requires picking a block size
that scales proportional to $\mathcal{O}\left(2^{p}/\epsilon^{1/3}\right)$ to
induce a constant gap between $\sigma_{1}$ and $\sigma_{b}$, and the number of
matrix-vector products scales exponentially in $p$.

Finally, in the above outline, we also crucially use that
$\left\|\mathbf{A}zz^{\top}\right\|_{F}^{2}=\left\|\mathbf{A}z\right\|_{2}^{2}$.
Observe that this no longer holds if we replace $z$ with a matrix $\mathbf{Z}$
that has $k$ orthonormal columns. Therefore, it remains unclear how to relate
$\left\|\mathbf{A}\mathbf{Z}\right\|^{p}_{\mathcal{S}_{p}}$ to
$\left\|\mathbf{A}\mathbf{Z}_{*,i}\right\|_{2}^{2}$, yet the vector-by-vector
error guarantee obtained by Block Krylov (see Lemmas 5.2 and 5.3) only bounds
the latter.

#### Handling all Schatten-$p$ Norms and $k>1$.

We modify our algorithm to run Block Krylov on $\mathbf{A}^{\top}$ and obtain
an orthonormal matrix $\mathbf{W}$ such that for all $i\in[k]$,

$\left\|\mathbf{A}^{\top}\mathbf{W}_{*,i}\right\|^{2}\geq\sigma_{i}^{2}(\mathbf{A})-\gamma\sigma^{2}_{k+1}(\mathbf{A}),$
(2.6)

for some $\gamma>0$. We then analyze the cost
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$,
where $\mathbf{Z}$ is a basis for $\mathbf{A}^{\top}\mathbf{W}$. Our key
insight is to interpret the input matrix $\mathbf{A}$ as a partitioned
operator (block matrix) and invoke _pinching inequalities_ for such operators.
Pinching inequalities were originally introduced to understand unitarily
invariant norms over direct sums of Hilbert spaces [VN37, Sch60]. In our
setting, given a block matrix
$\mathbf{M}=\begin{pmatrix}\mathbf{M}^{(1)}&\mathbf{M}^{(2)}\\\ \
\mathbf{M}^{(3)}&\mathbf{M}^{(4)}\end{pmatrix}$, the _pinching inequality_
(see Fact 4.13) implies that for all $p\geq 1$,

$\left\|\mathbf{M}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{M}^{(1)}\right\|_{\mathcal{S}_{p}}^{p}+\left\|\mathbf{M}^{(4)}\right\|_{\mathcal{S}_{p}}^{p}.$
(2.7)

A priori, it is unclear how to use Equation (2.7) to bound
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$.
First, we establish a general inequality for the Schatten norm of a matrix
times an orthogonal projection. Let $\mathbf{P}$ and $\mathbf{Q}$ be any
$n\times n$ orthogonal projection matrices with rank $k$ (see Definition 4.1).
Then, we prove (see Lemma 5.5 for details) that for any matrix $\mathbf{A}$,

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}^{p}+\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}.$
(2.8)

To obtain this inequality, we use a rotation argument along with the fact that
the Schatten-$p$ norms are unitarily invariant to show that
$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}=\left\|\begin{pmatrix}\mathbf{A}^{(1)}&\mathbf{A}^{(2)}\\\
\
\mathbf{A}^{(3)}&\mathbf{A}^{(4)}\end{pmatrix}\right\|_{\mathcal{S}_{p}}^{p}$,
where
$\left\|\mathbf{A}^{(1)}\right\|_{\mathcal{S}_{p}}=\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}$
and
$\left\|\mathbf{A}^{(4)}\right\|_{\mathcal{S}_{p}}=\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}$,
and then we can apply Equation (2.7) to the block matrix above.

Once we have established Equation (2.8), we can set
$\mathbf{P}=\mathbf{W}\mathbf{W}^{\top}$ and set
$\mathbf{Q}=\mathbf{Z}\mathbf{Z}^{\top}$ to be the projection matrix
corresponding to the column span of
$\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{\top}$. Then, we have that
$\mathbf{P}\mathbf{A}\mathbf{Q}=\mathbf{W}\mathbf{W}^{\top}\mathbf{A}$ and
$\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)=\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)$,
and combined with (2.8) this yields

$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{W}\mathbf{W}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}.$
(2.9)

To obtain a bound on
$\left\|\mathbf{W}\mathbf{W}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}$,
we appeal to the per-vector guarantees in Equation (2.6). However, translating
from $\ell_{2}^{2}$ error to
$\sigma_{p}^{p}\left(\mathbf{W}^{\top}\mathbf{A}\right)$ incurs a mixed
guarantee (see Lemma 5.7 for details):

$\left\|\mathbf{W}\mathbf{W}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\mathcal{O}\left(\gamma
p\right)\sum_{i\in[k]}\sigma_{k+1}^{2}\left(\mathbf{A}\right)\sigma_{i}^{p-2}\left(\mathbf{A}\right).$

To use this bound, we require $\sigma_{1}(\mathbf{A})$ to be comparable to
$\sigma_{k+1}(\mathbf{A})$ and thus we require an involved case analysis,
which appears in the proof of Theorem 5.1.

#### Avoiding an exponential dependence on $p$.

Our main insight here is that we do not require a block size that induces a
constant gap between singular values. Instead, we first observe that if the
block size $b$ is large enough such that $\sigma_{b}\leq\sigma_{2}/(1+1/p)$,
then $\mathcal{O}\left(\log(n/\epsilon)\sqrt{p}\right)$ iterations suffice to
obtain a vector $z$ such that
$\left\|\mathbf{A}z\right\|_{2}^{2}\geq\sigma_{1}^{2}\left(\mathbf{A}\right)-\left(\epsilon/n\right)\sigma_{2}^{2}\left(\mathbf{A}\right)$.
Therefore, we can trade-off the threshold for the Schatten norm of the tail
with the number of iterations as follows: if
$\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{\mathcal{S}_{p}}^{p}\leq\frac{1}{p^{1/3}\epsilon^{1/3}}\sigma_{2}^{p}\left(\mathbf{A}\right)$,
then setting $b=(1+1/p)^{p}/(\epsilon p)^{1/3}=\Theta(1/(\epsilon p)^{1/3})$
suffices to induce a gap of $1+1/p$ with block size $b$. The total number of
matrix-vector products is
$\mathcal{O}\left(b\cdot\log(n/\epsilon)\sqrt{p}\right)=\tilde{\mathcal{O}}(p^{1/6}/\epsilon^{1/3})$,
since $p$ can be assumed to be at most $(\log n)/\epsilon$. Otherwise,
$\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{\mathcal{S}_{p}}^{p}>\frac{1}{p^{1/3}\epsilon^{1/3}}\sigma_{2}^{p}\left(\mathbf{A}\right)$,
and we only require a $(1+\epsilon^{2/3}/p^{1/3})$-approximate solution
instead (compare with Equation (2.5)). Using gap-independent bounds (see Lemma
5.2), it suffices to start with block size $1$ and run
$\mathcal{O}\left(\log(n/\epsilon)p^{1/6}/\epsilon^{1/3}\right)$ iterations to
obtain a $(1+\epsilon^{2/3}/p^{1/3})$-approximate solution.

#### Avoiding a Gap-Dependent Bound.

We note that even when there is a constant gap between the first and second
singular values, and the per vector guarantee is highly accurate, i.e., for
all $i\in[k]$,
$\left\|\mathbf{A}\mathbf{Z}_{*,i}\right\|^{2}\geq\sigma^{2}_{i}(\mathbf{A})-\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2}(\mathbf{A})$,
it is not clear how to lower bound
$\left\|\mathbf{A}\mathbf{Z}\right\|_{\mathcal{S}_{p}}^{p}$ in Equation 2.9.
In general, the best bound we can obtain using the above equation is

$\left\|\mathbf{A}\mathbf{Z}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\mathcal{O}\left(\frac{\epsilon}{\textrm{poly}(d)}\right)\sigma_{k+1}^{2}\cdot\sum_{i\in[k]}\sigma_{i}^{p-2},$
(2.10)

which may be vacuous when the top $k$ singular values are significantly larger
than $\sigma_{k+1}$ and $p>2$. One could revert to a gap-dependent bound,
where the error is in terms of the gap between $\sigma_{1}$ and
$\sigma_{k+1}$, which one could account for by running an extra factor of
$\mathcal{O}\left(\log(\sigma_{1}/\sigma_{k+1})\right)$ iterations.

To avoid this gap-dependent bound, we split $\mathbf{A}$ into a head part
$\mathbf{A}_{H}$ and a tail part $\mathbf{A}_{T}$, such that $\mathbf{A}_{H}$
has all singular values that are at least $\left(1+1/d\right)\sigma_{k+1}$ and
$\mathbf{A}_{T}$ has the remaining singular values. We then bound
$\left\|\mathbf{A}_{H}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}$
and
$\left\|\mathbf{A}_{T}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}$
separately. Repeating the above analysis, we can obtain Equation (2.10) for
$\mathbf{A}_{T}$ instead, and since all singular values larger than
$\sigma_{k+1}$ in $\mathbf{A}_{T}$ are bounded, we can obtain
$\left\|\mathbf{A}_{T}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq\mathcal{O}\left(\epsilon
k/\textrm{poly}(d)\right)\sigma_{k+1}^{p}$. To adapt the analysis for
$\mathbf{A}_{T}$ and obtain this bound, we use Cauchy’s interlacing theorem to
relate the $j$-th singular value of
$\mathbf{A}_{T}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)$ to the
$(i^{*}+j)$-th singular value of
$\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)$, where $i^{*}$
is the rank of $\mathbf{A}_{H}$. We lower bound the $(i^{*}+j)$-th singular
value of $\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)$ using
the per vector guarantee of [MM15].

To bound
$\left\|\mathbf{A}_{H}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}$,
we observe it has rank at most $k$ and thus

$\left\|\mathbf{A}_{H}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq\sqrt{k}\cdot\left\|\mathbf{A}_{H}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{F}=\sqrt{k}\cdot\sqrt{\left\|\mathbf{A}_{H}\right\|_{F}^{2}-\left\|\mathbf{A}_{H}\mathbf{Z}\right\|_{F}^{2}},$

and we show how to bound this term in Section 5. Intuitively, while the
$k$-dimensional subspace that we find can “swap out" singular vectors
corresponding to singular values $\sigma_{i}$ for which $\sigma_{i}$ is very
close to $\sigma_{k+1}$, since they serve equally well for a Schatten-$p$ low-
rank approximation, for singular values $\sigma_{i}$ that are a bit larger
than $\sigma_{k+1}$, the $k$-dimensional subspace we find cannot do this. More
precisely, if $y$ is a singular vector of $\mathbf{A}_{H}$ with singular value
$\sigma_{i}$, then the projection of $y$ onto the $k$-dimensional subspace
that our algorithm finds (namely, $\mathbf{Z}$) must be at least
$1-\sigma_{k+1}^{2}/((\sigma_{i}^{2}-\sigma_{k+1}^{2})\textrm{poly}(d))$,
which suffices to bound the above since the additive error is inversely
proportional to $\sigma_{i}^{2}$ when $\sigma_{i}^{2}\gg\sigma_{k+1}^{2}$, and
so the very tiny additive error negates the effect of very large singular
values.

###  2.2  Query Lower Bounds.

Our lower bounds rely on the hardness of estimating the smallest eigenvalue of
a Wishart ensemble (see Definition 4.15), as established in recent work of
Braverman, Hazan, Simchowitz and Woodworth [BHSW20]. In particular, [BHSW20]
show that for a $d\times d$ instance $\mathbf{W}$ of a Wishart ensemble,
estimating $\lambda_{d}(\mathbf{W})$ (minimum eigenvalue) to additive error
$1/d^{2}$ requires $\Omega(d)$ adaptive matrix-vector product queries (see
Theorem 3.1 in [BHSW20]). To obtain hardness for Schatten-$p$ low-rank
approximation, we show that when $d=\Theta\left(1/\epsilon^{1/3}\right)$, any
candidate unit vector $z$ that satisfies
$\left\|\left(\mathbf{I}-\mathbf{W}/5\right)\left(\mathbf{I}-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq\left(1+\epsilon\right)\min_{\left\|u\right\|_{2}=1}\left\|\left(\mathbf{I}-\mathbf{W}/5\right)\left(\mathbf{I}-uu^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p},$
can be used to obtain an estimate
$\hat{\lambda}_{d}=\frac{5}{p}\left(1-\left\|\left(\mathbf{I}-\mathbf{W}/5\right)z\right\|_{2}^{p}\right)$
such that $\hat{\lambda}_{d}=(1\pm
1/d^{2})\lambda_{d}\left(\mathbf{I}-\mathbf{W}/5\right)$. Let
$\mathbf{A}=\left(\mathbf{I}-\mathbf{W}/5\right)$. To show our query lower
bound, in contrast to the analysis of our algorithm, the challenge is now to
lower bound
$\left\|\mathbf{A}\left(\mathbf{I}-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$
in terms of $\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}$ and
$\left\|\mathbf{A}z\right\|_{2}^{p}$ (contrast with Equation (2.9)).

#### Projection Cost via Araki-Lieb-Thirring.

First, we note that the case of $p=2$ is easy given the Pythagorean theorem.
For $p\in[1,2)$, we can establish an inequality fairly straightforwardly:
using the trace inner product definition of Schatten-$p$ (see Definition 4.7 )
norms, we have,

$\left\|\mathbf{A}\left(\mathbf{I}-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}=\mathrm{Tr}\left(\left(\left(\mathbf{I}-zz^{\top}\right)^{2}\mathbf{A}^{2}\left(\mathbf{I}-zz^{\top}\right)^{2}\right)^{p/2}\right),$
(2.11)

Since $p/2\in[1/2,1)$, we can use the reverse _Araki-Lieb-Thirring_ inequality
(see Fact 4.10) to show that

$\begin{split}\mathrm{Tr}\left(\left(\left(\mathbf{I}-zz^{\top}\right)^{2}\mathbf{A}^{2}\left(\mathbf{I}-zz^{\top}\right)^{2}\right)^{p/2}\right)&\geq\mathrm{Tr}\left(\left(\mathbf{I}-zz^{\top}\right)\mathbf{A}^{p}\left(\mathbf{I}-zz^{\top}\right)\right)\\\
&=\mathrm{Tr}\left(\mathbf{A}^{p}\right)-\mathrm{Tr}\left(\left(zz^{\top}\right)^{p/2}\left(\mathbf{A}^{2}\right)^{p/2}\left(zz^{\top}\right)^{p/2}\right)\\\
&\geq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{A}zz^{\top}\right\|_{\mathcal{S}_{p}}^{p}\end{split}$
(2.12)

where we use the cyclicity of the trace and again use reverse _Araki-Lieb-
Thirring_ (Fact 4.10) to show that

$\mathrm{Tr}\left(\left(zz^{\top}\right)^{\frac{p}{2}}\left(\mathbf{A}^{2}\right)^{\frac{p}{2}}\left(zz^{\top}\right)^{\frac{p}{2}}\right)\leq\mathrm{Tr}\left(\left(zz^{\top}\mathbf{A}^{2}zz^{\top}\right)^{p/2}\right)=\left\|\mathbf{A}zz^{\top}\right\|_{\mathcal{S}_{p}}^{p}.$

Since we have
$\left\|\mathbf{A}zz^{\top}\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{A}z\right\|_{2}^{p}$,
we conclude
$\left\|\mathbf{A}\left(\mathbf{I}-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{A}zz^{\top}\right\|_{2}^{p}.$
This approach only works for $p\in[1,2)$; for $p>2$ the application of _Araki-
Lieb-Thirring_ is reversed in Equation 2.12 (since $p/2>1$, see Fact 4.10) and
we no longer get a lower bound on the cost in Equation 2.11. We therefore
require a new approach.

#### Projection Cost via Norm Compression.

Recall, $z$ is the unit vector output by our candidate low-rank approximation
and let $y=\mathbf{A}z/\left\|\mathbf{A}z\right\|_{2}$. We yet again interpret
the input matrix $\mathbf{A}$ as a partitioned operator by considering the
projection of $\mathbf{A}$ onto $zz^{\top}$, $yy^{\top}$ and the projection
away from these rank-$1$ subspaces. In particular, let
$\mathbf{I}-yy^{\top}=\mathbf{Y}\mathbf{Y}^{\top}$, and
$\mathbf{I}-zz^{\top}=\mathbf{Z}\mathbf{Z}^{\top}$, where $\mathbf{Y}$ and
$\mathbf{Z}$ have orthonormal columns. Then, using a rotation argument, we
show that

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}=\left\|\begin{pmatrix}y^{\top}\mathbf{A}z&y^{\top}\mathbf{A}\mathbf{Z}\\\
\
\mathbf{Y}^{\top}\mathbf{A}z&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\end{pmatrix}\right\|_{\mathcal{S}_{p}}.$

We define the $p$-compression of $\mathbf{A}$, $\mathbf{C}_{\mathbf{A},p}$:

$\mathbf{C}_{\mathbf{A},p}=\begin{pmatrix}\left\|y^{\top}\mathbf{A}z\right\|_{\mathcal{S}_{p}}&\left\|y^{\top}\mathbf{A}\mathbf{Z}\right\|_{\mathcal{S}_{p}}\\\
\
\left\|\mathbf{Y}^{\top}\mathbf{A}z\right\|_{\mathcal{S}_{p}}&\left\|\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\right\|_{\mathcal{S}_{p}}\end{pmatrix}.$

To relate the norms of $\mathbf{A}$ and $\mathbf{C}_{\mathbf{A},p}$, we
consider Audenaert’s Norm Compression Conjecture [Aud08], a question in
functional analysis concerning operator inequalities (see also [AK12]):

######  Conjecture 2.2 (Schatten-$p$ Norm Compression).

Let $\mathbf{M}$ be a partitioned operator (block matrix) such that
$\mathbf{M}=\begin{pmatrix}\mathbf{M}_{1}&\mathbf{M}_{2}\\\
\mathbf{M}_{3}&\mathbf{M}_{4}\end{pmatrix}$. Let
$\mathbf{C}_{\mathbf{M},p}=\begin{pmatrix}\left\|\mathbf{M}_{1}\right\|_{\mathcal{S}_{p}}&\left\|\mathbf{M}_{2}\right\|_{\mathcal{S}_{p}}\\\
\left\|\mathbf{M}_{3}\right\|_{\mathcal{S}_{p}}&\left\|\mathbf{M}_{4}\right\|_{\mathcal{S}_{p}}\end{pmatrix}$
be a $2\times 2$ matrix that denotes the Schatten-$p$ compression of
$\mathbf{M}$ for any $p\geq 1$. Then,
$\left\|\mathbf{M}\right\|_{\mathcal{S}_{p}}\geq\left\|\mathbf{C}_{\mathbf{M},p}\right\|_{\mathcal{S}_{p}}$
if $1\leq p\leq 2$, and
$\left\|\mathbf{M}\right\|_{\mathcal{S}_{p}}\leq\left\|\mathbf{C}_{\mathbf{M},p}\right\|_{\mathcal{S}_{p}}$
if $2\leq p<\infty.$

We could simply appeal to this conjecture to obtain that for all $p>2$,

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}\leq\left\|\mathbf{C}_{\mathbf{A},p}\right\|_{\mathcal{S}_{p}}=\left\|\begin{pmatrix}\left\|yy^{\top}\mathbf{A}zz^{\top}\right\|_{\mathcal{S}_{p}}&\left\|yy^{\top}\mathbf{A}\left(I-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}\\\
\
\left\|\left(\mathbf{I}-yy^{\top}\right)\mathbf{A}zz^{\top}\right\|_{\mathcal{S}_{p}}&\left\|\left(\mathbf{I}-yy^{\top}\right)\mathbf{A}\left(\mathbf{I}-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}\end{pmatrix}\right\|_{\mathcal{S}_{p}}.$
(2.13)

However, for our choice of $y$,
$\left\|yy^{\top}\mathbf{A}\left(I-zz^{\top}\right)\right\|_{\mathcal{S}_{p}}=0$.
With padding and rotation arguments, we can then reduce our problem to a block
matrix where the blocks in each row are aligned, i.e., each row is a scalar
multiple of a fixed matrix (see Lemma 6.6). Then, we can use one of the few
special cases of Conjecture 2.2 for aligned operators which has actually been
proved, and appears in Fact 4.14. We can thus unconditionally obtain the
inequality in Equation (2.13).

Now that we have reduced to the case where we have a $2\times 2$ matrix with
$3$ non-zero entries, we would like to bound its Schatten-$p$ norm. We
explicitly compute the singular values of $\mathbf{C}_{\mathbf{A},p}$ (see
Fact 6.7 ), and then use the structure of the instance to directly lower bound
$\left\|\mathbf{A}z\right\|_{2}^{p}$ as follows:

$\left\|\mathbf{A}z\right\|_{2}^{p}+\left(1+\mathcal{O}\left(\epsilon^{2p/3}\right)\right)\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{C}_{\mathbf{A},p}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p},$
(2.14)

where the last inequality follows from Equation (2.13). Since we understand
the spectrum of the matrix $\mathbf{A}$, we can explicitly compute all the
terms in (2.14) above and show that we can obtain an accurate estimate of the
minimum singular value of $\mathbf{A}$ from
$\left\|\mathbf{A}z\right\|_{2}^{p}$. See details in Section 6.2.

##  3 Additional Related Work

Existing approaches to solve low-rank approximation problems under several
norms fall into two broad categories: iterative methods and linear sketching.
Iterative methods, such as Krylov subspace based methods, are captured by the
matrix-vector product framework, whereas linear sketching allows for the
choice of a matrix $\SS\in\mathbb{R}^{t\times n}$, where $t$ is the number of
“queries”, and then observes the product $\SS\cdot\mathbf{A}$ and so on (see
[Woo14] and references therein). The model has important applications to
streaming and distributed algorithms and several recent works have focused on
estimating spectral norms and the top singular values [AN13, LNW14a, LW16b,
BBK+21], estimating Schatten and Ky-Fan norms [LW16b, LW17, LW16a, BKKS19] and
low-rank approximation [CW13, MM13b, NN13, BDN15, Coh16].

In addition to studying unitarily invariant norms, such as the Schatten norm,
there also has been significant amount of work on studying low-rank
approximation under matrix $\ell_{p}$ norms [SWZ17, BBB+19, SWZ20, MW21] and
weighted low-rank approximation [SJ03, RSW16, BWZ19], settings in which the
problem is known to be NP-Hard. Finally, there has been a recent flurry of
work on sublinear time algorithms for low-rank approximation under various
structural assumptions on the input [MW17b, BW18, IVWW19, SW19, BCW20] and in
quantum-inspired models [KP16, CLW18, Tan19, RSML18, GLT18, GSLW19, CCHW20].

##  4 Preliminaries

Given an $n\times d$ matrix $\mathbf{A}$ with rank $r$, and $n\geq d$, we can
compute its singular value decomposition, denoted by
${SVD}(\mathbf{A})=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, such that
$\mathbf{U}$ is an $n\times r$ matrix with orthonormal columns,
$\mathbf{V}^{\top}$ is an $r\times d$ 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 any unitarily invariant norm, but in particular for any
Schatten-$p$ norm (defined below), we have
$\mathbf{A}_{k}=\min_{\textrm{rank}(\mathbf{X})=k}\|\mathbf{A}-\mathbf{X}\|_{\mathcal{S}_{p}}$.
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.

We use the notation $\mathbf{I}_{k}$ to denote a _truncated identity matrix_ ,
that is, a square matrix with its top $k$ diagonal entries equal to one, and
all other entries zero. The dimension of $\mathbf{I}_{k}$ will be determined
by context.

######  Definition 4.1 (Orthogonal Projection Matrices).

Given a $d\times d$ symmetric matrix $\mathbf{P}$ and $k\in[d]$, $\mathbf{P}$
is a rank-$k$ orthogonal projection matrix if $\textrm{rank}(\mathbf{P})=k$
and $\mathbf{P}^{2}=\mathbf{P}$.

It follows from the above definition that $\mathbf{P}$ has eigenvalues that
are either $0$ or $1$ and admits a singular value decomposition of the form
$\mathbf{U}\mathbf{U}^{\top}$ where $\mathbf{U}$ has $k$ orthonormal columns.

######  Definition 4.2 (Unitary Matrices).

Given a symmetric matrix $\mathbf{U}\in\mathbb{R}^{d\times d}$ we say
$\mathbf{U}$ is a unitary matrix if
$\mathbf{U}^{\top}\mathbf{U}=\mathbf{U}\mathbf{U}^{\top}=\mathbf{I}$.

######  Definition 4.3 (Rotation Matrices).

Given a symmetric matrix $\mathbf{R}\in\mathbb{R}^{d\times d}$ we say
$\mathbf{R}$ is a rotation matrix if $\mathbf{R}$ is unitary and
$\textrm{det}\left(\mathbf{R}\right)=1$.

######  Fact 4.4 (Courant-Fischer for Singular Values).

Given an $n\times d$ matrix $\mathbf{A}$ with singular values
$\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{d}$, the following holds: for
all $i\in[d]$,

$\sigma_{i}=\max_{S:\hskip 3.61371pt\textrm{dim}(S)=i}\hskip
7.22743pt\min_{x\in S:\hskip 3.61371pt\left\|x\right\|_{2}=1}\hskip
7.22743pt\left\|x^{\top}\mathbf{A}\right\|_{2}.$

######  Fact 4.5 (Weyl’s Inequality for Singular Values (see Exercise 22
[Tao20])).

Given $n\times d$ matrices $\mathbf{X},\mathbf{Y}$, for any $i,(j-1)\in[d]$
such that $i+j\leq d$ ,

$\sigma_{i+j}\left(\mathbf{X}+\mathbf{Y}\right)\leq\sigma_{i}(\mathbf{X})+\sigma_{j+1}(\mathbf{Y}).$

######  Fact 4.6 (Bernoulli’s Inequality).

For any $x,p\in\mathbb{R}$ such that $x\geq-1$ and $p\geq 1$,
$\left(1+x\right)^{p}\geq 1+px$.

#### Schatten Norms and Trace Inequalities.

We recall some basic facts for Schatten-$p$ norms. We also require the
following trace and operator inequalities.

######  Definition 4.7 (Schatten-$p$ Norm).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, let
$\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{d}$ be the singular values of
$\mathbf{A}$. Then, for any $p\in[0,\infty)$, the Schatten-$p$ norm of
$\mathbf{A}$ is defined as

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}=\mathrm{Tr}\left(\left(\mathbf{A}^{\top}\mathbf{A}\right)^{p/2}\right)^{1/p}=\left(\sum_{i\in[d]}\sigma_{i}^{p}(\mathbf{A})\right)^{1/p}.$

######  Fact 4.8 (Schatten-$p$ norms are Unitarily Invariant).

Given an $n\times d$ matrix $\mathbf{M}$, for any $m\times n$ matrix
$\mathbf{U}$ with orthonormal columns, a norm $\|\cdot\|_{X}$ is defined to be
unitarily invariant if $\|\mathbf{U}\mathbf{M}\|_{X}=\|\mathbf{M}\|_{X}$. The
Schatten-$p$ norm is unitarily invariant for all $p\geq 1$.

There exists a closed-form expression for the low-rank approximation problem
under Schatten-$p$ norms:

######  Fact 4.9 (Schatten-$p$ Low-Rank Approximation).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and an integer
$k\in\mathbb{N}$,

$\mathbf{A}_{k}=\arg\min_{\textrm{rank}(\mathbf{X})\leq
k}\left\|\mathbf{A}-\mathbf{X}\right\|_{\mathcal{S}_{p}},$

where $\mathbf{A}_{k}$ is the truncated SVD of $\mathbf{A}$.

######  Fact 4.10 (Araki–Lieb–Thirring Inequality [Ara90]).

Given PSD matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{d\times d}$, for any
$r\geq 1$, the following inequality holds:

$\mathrm{Tr}\left(\left(\mathbf{B}\mathbf{A}\mathbf{B}\right)^{r}\right)\leq\mathrm{Tr}\left(\mathbf{B}^{r}\mathbf{A}^{r}\mathbf{B}^{r}\right).$

Further, for $0<r<1$, the reverse holds

$\mathrm{Tr}\left(\left(\mathbf{B}\mathbf{A}\mathbf{B}\right)^{r}\right)\geq\mathrm{Tr}\left(\mathbf{B}^{r}\mathbf{A}^{r}\mathbf{B}^{r}\right).$

######  Fact 4.11 (Mahler’s Orthogonal Operator Inequality, Theorem 1.7 in
[Mah90]).

Given $p\geq 2$, and matrices $\mathbf{P}$ and $\mathbf{Q}$ such that the row
(column) span of $\mathbf{P}$ is orthogonal to the row (column) span of
$\mathbf{Q}$, the following inequality holds:

$\left\|\mathbf{P}\right\|^{p}_{\mathcal{S}_{p}}+\left\|\mathbf{Q}\right\|^{p}_{\mathcal{S}_{p}}\leq\left\|\mathbf{P}+\mathbf{Q}\right\|^{p}_{\mathcal{S}_{p}}.$

######  Fact 4.12 (Hölder’s Inequality for Schatten-$p$ Norms, Corollary 4.2.6
[Bha13]).

Given matrices $\mathbf{A},\mathbf{B}^{\top}\in\mathbb{R}^{n\times d}$ and
$p\in[1,\infty)$, the following holds

$\left\|\mathbf{A}\mathbf{B}\right\|_{\mathcal{S}_{p}}\leq\left\|\mathbf{A}\right\|_{\mathcal{S}_{q}}\cdot\left\|\mathbf{B}\right\|_{\mathcal{S}_{r}},$

for any $q,r$ such that $\frac{1}{p}=\frac{1}{q}+\frac{1}{r}$.

We also require _pinching inequalities_ that were originally introduced to
relate norms for partitioned operators over direct sums of Hilbert spaces. In
our context, these inequalities simplify to norm inequalities for block
matrices:

######  Fact 4.13 (Pinching Inequalities for Schatten-$p$ Norms, [BKL02]).

Let $\mathbf{M}\in\mathbb{R}^{td\times td}$ be the following block matrix

$\mathbf{M}=\left[\begin{array}[]{cccc}\mathbf{M}_{(1,1)}&\mathbf{M}_{(1,2)}&\cdots&\mathbf{M}_{(1,t)}\\\
\mathbf{M}_{(2,1)}&\mathbf{M}_{(2,2)}&\cdots&\mathbf{M}_{(1,t)}\\\
\vdots&&\ddots&\vdots\\\
\mathbf{M}_{(t,1)}&\mathbf{M}_{(t,2)}&\cdots&\mathbf{M}_{(t,t)}\\\
\end{array}\right],$

where for all $i,j\in[t]$, $\mathbf{M}_{(i,j)}\in\mathbb{R}^{d\times d}$. For
all $p\geq 1$, the following inequality holds:

$\left(\sum_{i\in[t]}\left\|\mathbf{M}_{(i,i)}\right\|_{\mathcal{S}_{p}}^{p}\right)^{1/p}\leq\left\|\mathbf{M}\right\|_{\mathcal{S}_{p}}.$

We also require a norm compression inequality that is a special case of
Conjecture 2.2 (and known to be true), when each block is aligned in the
following sense:

######  Fact 4.14 (Aligned Norm Compression Inequality, Section 4.3 in
[Aud08]).

Let $\mathbf{M}=\begin{pmatrix}\mathbf{M}_{1}&\mathbf{M}_{2}\\\
\mathbf{M}_{3}&\mathbf{M}_{4}\end{pmatrix}$ such that there exist scalars
$\alpha_{1},\alpha_{2},\beta_{1},\beta_{2}$ such that
$\mathbf{M}_{1}=\alpha_{1}\mathbf{X}$, $\mathbf{M}_{2}=\alpha_{2}\mathbf{X}$,
$\mathbf{M}_{3}=\beta_{1}\mathbf{Y}$and $\mathbf{M}_{4}=\beta_{2}\mathbf{Y}$.
Then, for any $p\geq 2$,

$\left\|\mathbf{M}\right\|_{\mathcal{S}_{p}}\leq\left\|\begin{pmatrix}\left\|\mathbf{M}_{1}\right\|_{\mathcal{S}_{p}}&\left\|\mathbf{M}_{2}\right\|_{\mathcal{S}_{p}}\\\
\left\|\mathbf{M}_{3}\right\|_{\mathcal{S}_{p}}&\left\|\mathbf{M}_{4}\right\|_{\mathcal{S}_{p}}\end{pmatrix}\right\|_{\mathcal{S}_{p}}.$

#### Random Matrix Theory.

Next, we recall some basic facts for Wishart ensembles from random matrix
theory (we refer the reader to [Tao12] for a comprehensive overview).

######  Definition 4.15 (Wishart Ensemble).

An $n\times n$ matrix $\mathbf{W}$ is sampled from a Wishart Ensemble,
$\textsf{Wishart}(n)$, if $\mathbf{W}=\mathbf{X}\mathbf{X}^{\top}$ such that
for all $i,j\in[n]$
$\mathbf{X}_{i,j}\sim\mathcal{N}\left(0,\frac{1}{n}\mathbf{I}\right)$.

######  Fact 4.16 (Norms of a Wishart Ensemble).

Let $\mathbf{W}\sim\textsf{Wishart}(n)$ such that
$n=\Omega(1/\varepsilon^{3})$. Then, with probability $99/100$,
$\left\|\mathbf{W}\right\|_{\textrm{op}}\leq 5$ and for any fixed constant
$p$,
$\left\|\mathbf{I}-\frac{1}{5}\mathbf{W}\right\|^{p}_{\mathcal{S}_{p}}=\Theta\left(\frac{1}{\varepsilon^{1/3}}\right)$.

##  5 Algorithms for Schatten-$p$ LRA

In this section, we focus on obtaining algorithms for low-rank approximation
in Schatten-$p$ norm, simultaneously for all real, not necessarily constant,
$p\in[1,\mathcal{O}\left(\log(d)/\epsilon\right)]$. For the special case of
$p\in\\{2,\infty\\}$, Musco and Musco [MM15] showed an algorithm with matrix-
vector query complexity $\tilde{O}(k/\epsilon^{1/2})$, given below as
Algorithm 5.6. We show that the number of matrix-vector products we require
scales proportional to $\tilde{O}\left(kp^{1/6}/\epsilon^{1/3}\right)$
instead. Finally, recall when $p>\log(d)/\epsilon$, it suffices to run Block
Krylov for $p=\infty$, which requires
$\mathcal{O}\left(\log(d/\epsilon)k/\sqrt{\epsilon}\right)$ matrix-vector
products.

######  Theorem 5.1 (Optimal Schatten-$p$ Low-Rank Approximation).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, a target rank $k\in[d]$,
an accuracy parameter $\epsilon\in(0,1)$ and any
$p\in[1,\mathcal{O}\left(\log(d)/\epsilon\right)]$, Algorithm 5.4 performs
$\mathcal{O}\left(\frac{kp^{1/6}\log(d/\epsilon)}{\epsilon^{1/3}}+\log(d/\epsilon)k\sqrt{p}\right)$
matrix-vector products and outputs a $d\times k$ matrix $\mathbf{Z}$ with
orthonormal columns such that with probability at least $9/10$,

$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq\left(1+\epsilon\right)\min_{\mathbf{U}:\hskip
3.61371pt\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{k}}\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top}\right)\
\right\|_{\mathcal{S}_{p}}.$

Further, in the RAM model, the algorithm runs in time
$\mathcal{O}\left(\frac{\operatorname{\mathtt{nnz}}(\mathbf{A})p^{1/6}k\log^{2}(d/\epsilon)}{\epsilon^{1/3}}+\frac{np^{(\omega-1)/6}k^{\omega-1}}{\epsilon^{(\omega-1)/3}}\right)$.

We first introduce the following lemmas from Musco and Musco [MM15] that
provide convergence bounds for the performance of Block Krylov Iteration
(Algorithm 5.6) :

######  Lemma 5.2 (Gap Independent Block Krylov with Arbitrary Accuracy).

Let $\mathbf{A}$ be an $n\times d$ matrix, $k$ be the target rank and
$\gamma>0$ be an accuracy parameter. Then, initializing Algorithm 5.6 with
block size $k$ and running for
$q=\Omega\left(\log(d/\gamma)/\sqrt{\gamma}\right)$ iterations outputs a
$d\times k$ matrix $\mathbf{Z}$ such that with probability $99/100$, for all
$i\in[k]$,

$\left\|\mathbf{A}\mathbf{Z}_{*,i}\right\|^{2}_{2}=\sigma_{i}^{2}\pm\gamma\sigma^{2}_{k+1}.$

Further, the total number of matrix-vector products is
$\mathcal{O}\left(kq\right)$ and the running time in the RAM model is
$\mathcal{O}\left(\operatorname{\mathtt{nnz}}(\mathbf{A})kq+n\left(kq\right)^{2}+\left(kq\right)^{\omega}\right)$.

The aforementioned lemma follows directly from Theorem 1 in [MM15], using the
per-vector error guarantee (3).

######  Lemma 5.3 (Gap Dependent Block Krylov, Theorem 13 [MM15]).

Let $\mathbf{A}$ be an $n\times d$ matrix and $\gamma>0$, be an accuracy
parameter and $p,k\in\mathcal{N}$ be such that $b\geq k$. Let
$\sigma_{1},\sigma_{2}\ldots\sigma_{d}$ be the singular values of
$\mathbf{A}$. Then, initializing Algorithm 5.6 with block size $b$ and running
for
$q=\Omega\left(\log(n/\gamma)\sqrt{\sigma_{k}}/\sqrt{\sigma_{k}-\sigma_{b}}\right)$
iterations outputs a $d\times k$ matrix $\mathbf{Z}$ such that with
probability $99/100$, for all $i\in[k]$

$\left\|\mathbf{A}\mathbf{Z}_{*,i}\right\|^{2}_{2}=\sigma^{2}_{i}\pm\gamma\sigma^{2}_{k+1}.$

Further, the total number of matrix-vector products is
$\mathcal{O}\left(pq\right)$ and the running time in the RAM model is
$\mathcal{O}\left(\operatorname{\mathtt{nnz}}(\mathbf{A})bq+n\left(bq\right)^{2}+\left(bq\right)^{\omega}\right)$.

######  Algorithm 5.4 (Optimal Schatten-$p$ Low-rank Approximation).

Input: An $n\times d$ matrix $\mathbf{A}$, target rank $k\leq d$, accuracy
parameter $0<\varepsilon<1$, and $p\geq 1$. 1. Let
$\gamma_{1}=\varepsilon^{2/3}/p^{1/3}$. Run Block Krylov Iteration (Algorithm
5.6) on $\mathbf{A}^{\top}$ with block size $k$, and number of iterations
$q=\mathcal{O}\left(\log(d/\gamma_{1})/\sqrt{\gamma_{1}}+\log(d/\epsilon)\sqrt{p}\right)$.
Let $\mathbf{W}_{1}\in\mathbb{R}^{n\times k}$ be the corresponding output with
orthonormal columns. 2. Let $\gamma_{2}=\varepsilon$ and let
$s=\mathcal{O}\left(p^{-1/3}k/\varepsilon^{1/3}\right)$. Run Block Krylov
Iteration (Algorithm 5.6) on $\mathbf{A}^{\top}$ with block size $s$, and
number of iterations $q=\mathcal{O}\left(\log(d/\gamma_{2})\sqrt{p}\right)$.
Let $\mathbf{W}_{2}\in\mathbb{R}^{n\times k}$ be the corresponding output with
orthonormal columns. 3. Run Block Krylov on $\mathbf{A}$ with target rank
$k+1$ and number of iterations
$q=\mathcal{O}\left((\log(dp)+\log(d/\epsilon))\sqrt{p}\right)$, and let
$\hat{\mathbf{Z}}_{1}$ be the resulting $d\times(k+1)$ output matrix. Compute
$\hat{\sigma}^{2}_{1}=\left\|\mathbf{A}(\hat{\mathbf{Z}}_{1})_{*,1}\right\|_{2}^{2}$
and
$\hat{\sigma}^{2}_{k+1}=\left\|\mathbf{A}(\hat{\mathbf{Z}}_{1})_{*,k+1}\right\|_{2}^{2}$,
rough estimates of the $1$-st and $(k+1)$-st singular values of $\mathbf{A}$.
Run Block Krylov on $\mathbf{A}$ with target rank $s$, where
$s=\mathcal{O}\left(p^{-1/3}k/\varepsilon^{1/3}\right)$ and iterations
$q=\mathcal{O}\left(\log(d/\epsilon)\sqrt{p}\right)$, and let
$\hat{\mathbf{Z}}_{2}$ be the resulting $d\times s$ output matrix. Compute
$\hat{\sigma}^{2}_{s}=\left\|\mathbf{A}(\hat{\mathbf{Z}}_{2})_{*,s}\right\|_{2}^{2}$,
an estimate to the $s$-th singular value of $\mathbf{A}$. 4. If
$\hat{\sigma}_{1}^{2}\geq(1+0.5/p)\hat{\sigma}_{k+1}^{2}$, set
$\mathbf{Z}=\mathbf{Z}_{1}$. Else, if
$\hat{\sigma}^{2}_{s}\leq\hat{\sigma}^{2}_{k+1}/\left(1+0.5/p\right)$, set
$\mathbf{Z}$ to be an orthonormal basis for
$\mathbf{A}^{\top}\mathbf{W}_{2}\mathbf{W}_{2}^{\top}$ and otherwise set
$\mathbf{Z}$ to be an orthonormal basis for
$\mathbf{A}^{\top}\mathbf{W}_{1}\mathbf{W}_{1}^{\top}$. Output: A matrix
$\mathbf{Z}\in\mathbb{R}^{d\times k}$ with orthonormal columns such that
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq\left(1+\epsilon\right)\min_{\mathbf{U}:\hskip
3.61371pt\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{k}}\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}.$

Next, we prove the following key lemma relating the Schatten-$p$ norm of row
and column projections applied to a matrix $\mathbf{A}$ to the Schatten-$p$
norm of the matrix itself. We can interpret this lemma as an extension of the
Pythagorean Theorem to Schatten-$p$ spaces and believe this lemma is of
independent interest. We note that we appeal to _pinching inequality_ for
partitioned operators to obtain this lemma.

######  Lemma 5.5 (Schatten-$p$ Norms for Orthogonal Projections).

Let $\mathbf{A}$ be an $n\times d$ matrix, let $\mathbf{P}$ be an $n\times n$
matrix, and let $\mathbf{Q}$ be a $d\times d$ matrix such that both
$\mathbf{P}$ and $\mathbf{Q}$ are orthogonal projection matrices of rank $k$
(see Definition 4.1). Then, the following inequality holds for all $p\geq 1$:

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}^{p}+\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}.$

###### Proof.

Let $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ be the SVD of
$\mathbf{A}$, where $\mathbf{U}\in\mathbb{R}^{n\times d}$ and
$\mathbf{V}^{\top}\in\mathbb{R}^{d\times d}$ have orthonormal columns and rows
respectively. We construct unitary matrices $\mathbf{R}$ and $\SS$, such that
$\mathbf{R}\in\mathbb{R}^{n\times n}$ and $\SS\in\mathbb{R}^{d\times d}$ that
satisfy the following constraints:

  1. 1.

$\mathbf{R}^{\top}\mathbf{I}_{k}\mathbf{R}\mathbf{A}\SS^{\top}\mathbf{I}_{k}\SS=\mathbf{P}\mathbf{A}\mathbf{Q}$,
and

  2. 2.

$\mathbf{R}^{\top}\left(\mathbf{I}-\mathbf{I}_{k}\right)\mathbf{R}\mathbf{A}\SS^{\top}\left(\mathbf{I}-\mathbf{I}_{k}\right)\SS=\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)$,

where the trunctated Identity matrix, $\mathbf{I}_{k}$, left multiplying
$\mathbf{A}$ is $n\times n$ and right multiplying $\mathbf{A}$ is $d\times d$.

Recall, since $\mathbf{P}$ is a rank-$k$ projection matrix, it admits a
decomposition $\mathbf{P}=\mathbf{X}\mathbf{X}^{\top}$ such that $\mathbf{X}$
has $k$ orthonormal columns and similarly
$\mathbf{I}-\mathbf{P}=\mathbf{Y}\mathbf{Y}^{\top}$, where $\mathbf{Y}$ has
$n-k$ orthonormal columns. Further, since $\mathbf{X}$ and $\mathbf{Y}$ span
disjoint subspaces, and the union of their span is $\mathbb{R}^{n}$, the
matrix $\left(\mathbf{X}\mid\mathbf{Y}\right)$, obtained by concatenating
their columns, is unitary. Then, it suffices to set
$\mathbf{R}=\left(\mathbf{X}\mid\mathbf{Y}\right)^{\top}$. To see this,
observe,

$\mathbf{R}^{\top}\mathbf{I}_{k}\mathbf{R}=\left(\mathbf{X}\mid
0\right)\cdot\begin{pmatrix}\mathbf{X}^{\top}\\\
0\end{pmatrix}=\mathbf{X}\mathbf{X}^{\top}=\mathbf{P},$

and similarly,

$\mathbf{R}^{\top}\left(\mathbf{I}-\mathbf{I}_{k}\right)\mathbf{R}=\mathbf{Y}\mathbf{Y}^{\top}=\mathbf{I}-\mathbf{P}.$

We repeat the above argument for the projection matrix $\mathbf{Q}$. Let
$\mathbf{Q}=\mathbf{W}\mathbf{W}^{\top}$, where $\mathbf{W}$ is $d\times k$
and has orthonormal columns, and
$\mathbf{I}-\mathbf{Q}=\mathbf{Z}\mathbf{Z}^{\top}$, where $\mathbf{Z}$ is
$d\times(d-k)$ and has orthonormal columns. Observe, it suffices to set
$\SS=\left(\mathbf{W}\mid\mathbf{Z}\right)^{\top}$, since $\SS$ is unitary and
$\SS^{\top}\mathbf{I}_{k}\SS=\mathbf{Q}$ and
$\SS^{\top}\left(\mathbf{I}-\mathbf{I}_{k}\right)\SS=\mathbf{I}-\mathbf{Q}$.
Note, by construction, we satisfy the two aforementioned constraints.

Let $\hat{\mathbf{A}}=\mathbf{R}\mathbf{A}\SS^{\top}$. Since $\mathbf{R}$ and
$\SS$ are unitary, it follows from unitary invariance of the Schatten-$p$ norm
that

$\left\|\hat{\mathbf{A}}\right\|_{\mathcal{S}_{p}}=\left\|\mathbf{R}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}\SS^{\top}\right\|_{\mathcal{S}_{p}}=\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}$
(5.1)

Further, observe for any $n\times d$ matrix $\mathbf{M}$, we have have the
following block decomposition

$\begin{split}\mathbf{M}&=\mathbf{I}_{k}\mathbf{M}\mathbf{I}_{k}+\mathbf{I}_{k}\mathbf{M}\left(\mathbf{I}-\mathbf{I}_{k}\right)+\left(\mathbf{I}-\mathbf{I}_{k}\right)\mathbf{M}\mathbf{I}_{k}+\left(\mathbf{I}-\mathbf{I}_{k}\right)\mathbf{M}\left(\mathbf{I}-\mathbf{I}_{k}\right)\\\
&=\begin{pmatrix}\mathbf{M}_{1:k,1:k}&\mathbf{M}_{1:k,k+1:d}\\\
\mathbf{M}_{k+1:n,1:k}&\mathbf{M}_{k+1:n,k+1:d}\end{pmatrix},\end{split}$

where the notation $\mathbf{M}_{i:i^{\prime},j:j^{\prime}}$ picks the
$(i^{\prime}-i+1)\times(j^{\prime}-j+1)$ sized sub-matrix corresponding to the
rows indices $[i,i^{\prime}]$ and column indices $[j,j^{\prime}]$. Since
appending rows and columns of $0$’s does not change the singular values, we
have
$\left\|\mathbf{I}_{k}\mathbf{M}\mathbf{I}_{k}\right\|_{\mathcal{S}_{p}}=\left\|\mathbf{M}_{1:k,1:k}\right\|_{\mathcal{S}_{p}}$
and
$\left\|\left(\mathbf{I}-\mathbf{I}_{k}\right)\mathbf{M}\left(\mathbf{I}-\mathbf{I}_{k}\right)\right\|_{\mathcal{S}_{p}}=\left\|\mathbf{M}_{k+1:n,k+1:d}\right\|_{\mathcal{S}_{p}}$.
Setting $\mathbf{M}=\hat{\mathbf{A}}$, we have

$\begin{split}\left\|\hat{\mathbf{A}}\right\|_{\mathcal{S}_{p}}^{p}&=\left\|\begin{pmatrix}\hat{\mathbf{A}}_{1:k,1:k}&\hat{\mathbf{A}}_{1:k,k+1:d}\\\
\hat{\mathbf{A}}_{k+1:n,1:k}&\hat{\mathbf{A}}_{k+1:n,k+1:d}\end{pmatrix}\right\|_{\mathcal{S}_{p}}^{p}\\\
&\geq\left\|\hat{\mathbf{A}}_{1:k,1:k}\right\|_{\mathcal{S}_{p}}^{p}+\left\|\hat{\mathbf{A}}_{k+1:n,k+1:d}\right\|_{\mathcal{S}_{p}}^{p}\\\
&=\left\|\mathbf{I}_{k}\hat{\mathbf{A}}\mathbf{I}_{k}\right\|_{\mathcal{S}_{p}}^{p}+\left\|\left(\mathbf{I}-\mathbf{I}_{k}\right)\hat{\mathbf{A}}\left(\mathbf{I}-\mathbf{I}_{k}\right)\right\|_{\mathcal{S}_{p}}^{p},\end{split}$
(5.2)

where the inequality follows from using the _pinching inequality_ on the block
matrix (see Fact 4.13). By the unitary invariance of the Schatten-$p$ norm, we
have

$\left\|\mathbf{I}_{k}\hat{\mathbf{A}}\mathbf{I}_{k}\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{R}^{\top}\mathbf{I}_{k}\hat{\mathbf{A}}\mathbf{I}_{k}\SS\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}^{p},$

and similarly,

$\left\|\left(\mathbf{I}-\mathbf{I}_{k}\right)\hat{\mathbf{A}}\left(\mathbf{I}-\mathbf{I}_{k}\right)\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{R}^{\top}\left(\mathbf{I}-\mathbf{I}_{k}\right)\hat{\mathbf{A}}\left(\mathbf{I}-\mathbf{I}_{k}\right)\SS\right\|_{\mathcal{S}_{p}}^{p}=\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}.$

Plugging these two bounds back into Equation (5.2), along with Equation (5.1),
we can conclude,

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}^{p}+\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}.$

∎

######  Algorithm 5.6 (Block Krylov Iteration, [MM15]).

Input: An $n\times d$ matrix $\mathbf{A}$, target rank $k$, iteration count
$q$ and a block size parameter $s$ such that $k\leq s\leq d$. 1. Let
$\mathbf{U}$ be a $n\times s$ matrix such that each entry is drawn i.i.d. from
$\mathcal{N}(0,1)$. Let
$\mathbf{K}=\left[\mathbf{A}^{\top}\mathbf{U};(\mathbf{A}^{\top}\mathbf{A})\mathbf{A}^{\top}\mathbf{U};(\mathbf{A}^{\top}\mathbf{A})^{2}\mathbf{A}^{\top}\mathbf{U};\ldots;(\mathbf{A}^{\top}\mathbf{A})^{q}\mathbf{A}^{\top}\mathbf{U}\right]$
be the $d\times s(q+1)$ Krylov matrix obtained by concatenating the matrices
$\mathbf{A}^{\top}\mathbf{U},\ldots,\left(\mathbf{A}^{\top}\mathbf{A}\right)^{q}\mathbf{A}^{\top}\mathbf{U}$.
2. Compute an orthonomal basis $\mathbf{Q}$ for the column span of
$\mathbf{K}$. Let
$\mathbf{M}=\mathbf{Q}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{Q}$. 3.
Compute the top $k$ left singular vectors of $\mathbf{M}$, and denote them by
$\mathbf{Y}_{k}$. Output: $\mathbf{Z}=\mathbf{Q}\mathbf{Y}_{k}$

Note, despite establishing Lemma 5.5, it is not immediately apparent how to
lower bound
$\left\|\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{\mathcal{S}_{p}}^{p}$,
where $\mathbf{Z}$ is a candidate solution. Next, we show how to translate a
guarantee on the Euclidean norm of $\mathbf{A}$ times a column of $\mathbf{Z}$
to a lower bound on
$\left\|\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{\mathcal{S}_{p}}^{p}$.

######  Lemma 5.7 (Per-Vector Guarantees to Schatten Norms).

Let $\mathbf{A}$ be an $n\times d$ matrix with singular values denoted by
$\\{\sigma_{i}\left(\mathbf{A}\right)\\}_{i\in[d]}$. Let $\mathbf{Z}$ be a
$d\times k$ matrix with orthonormal columns that is output by Algorithm 5.6,
such that for all $i\in[k]$, with probability at least 99/100,
$\left\|\mathbf{A}\mathbf{Z}_{*,i}\right\|_{2}^{2}\geq\sigma_{i}^{2}\left(\mathbf{A}\right)-\gamma_{i}\sigma_{k+1}^{2}\left(\mathbf{A}\right)$,
for some $\gamma\in(0,1)$. Then, for any $p\geq 1$, we have

$\left\|\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\sum_{i\in[k]}\mathcal{O}\left(\gamma_{i}p\right)\sigma_{k+1}^{2}\left(\mathbf{A}\right)\sigma_{i}^{p-2}\left(\mathbf{A}\right).$

###### Proof.

First, we observe that it suffices to show that
$\sigma_{i}(\mathbf{A}\mathbf{Z})^{2}\geq\left\|\mathbf{A}z_{i}\right\|_{2}^{2}$,
where $z_{i}$ is shorthand for $\mathbf{Z}_{*,i}$, the $i$-th columm of
$\mathbf{Z}$. Assuming this inequality holds, we can complete the proof as
follows: we know that for all $i\in[k]$,

$\begin{split}\sigma_{i}^{2}(\mathbf{A}\mathbf{Z})\geq\left\|\mathbf{A}z_{i}\right\|_{2}^{2}&\geq\sigma_{i}^{2}(\mathbf{A})-\gamma\sigma_{k+1}^{2}(\mathbf{A})\\\
&=\sigma_{i}^{2}(\mathbf{A})\left(1-\gamma\frac{\sigma_{k+1}^{2}(\mathbf{A})}{\sigma_{i}^{2}(\mathbf{A})}\right)\end{split}$
(5.3)

Then, taking $p/2$-th powers in (5.3),

$\begin{split}\sigma_{i}^{p}(\mathbf{A}\mathbf{Z})&\geq\sigma_{i}^{p}(\mathbf{A})\left(1-\gamma\frac{\sigma_{k+1}^{2}(\mathbf{A})}{\sigma_{i}^{2}(\mathbf{A})}\right)^{p/2}\\\
&\geq\sigma_{i}^{p}(\mathbf{A})\left(1-\mathcal{O}\left(\frac{\gamma
p\sigma_{k+1}^{2}(\mathbf{A})}{\sigma_{i}^{2}(\mathbf{A})}\right)\right)\\\
&=\sigma_{i}^{p}(\mathbf{A})-\mathcal{O}\left(\gamma
p\right)\sigma_{k+1}^{2}\left(\mathbf{A}\right)\sigma_{i}^{p-2}\left(\mathbf{A}\right)\end{split}$
(5.4)

where the second inequality follows from the generalized Bernoulli inequality
(see Fact 4.6). Summing over all $i\in[k]$, we can conclude

$\left\|\mathbf{A}\mathbf{Z}\right\|_{\mathcal{S}_{p}}^{p}\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\sum_{i\in[k]}\mathcal{O}\left(\gamma
p\right)\sigma_{k+1}^{2}\left(\mathbf{A}\right)\sigma_{i}^{p-2}\left(\mathbf{A}\right).$

Therefore, it remains to show that
$\sigma_{i}(\mathbf{A}\mathbf{Z})^{2}\geq\left\|\mathbf{A}z_{i}\right\|_{2}^{2}$.
First, we recall that Algorithm 5.6 outputs $\\{z_{i}\\}_{i\in[k]}$ such that
$z_{i}=\mathbf{Q}\tilde{z}_{i}$, where $\mathbf{Q}$ is an orthonormal basis
for the Krylov space $\mathbf{K}$ (an $d\times s(q+1)$ matrix) and
$\tilde{z}_{i}$ is the $i$-th singular vector of
$\mathbf{Q}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{Q}$. Note that the
$\tilde{z}_{i}$’s are $s(q+1)$-dimensional vectors. Let
$\mathbf{W}\mathbf{\Omega}\mathbf{W}^{\top}$ be the SVD of
$\mathbf{Q}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{Q}$. Then,
$\mathbf{Q}\mathbf{W}\mathbf{\Omega}\mathbf{W}^{\top}\mathbf{Q}^{\top}$ is the
SVD of
$\mathbf{Q}\mathbf{Q}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{Q}\mathbf{Q}^{\top}$.
To see this, let the $i$-th column of $\mathbf{Q}\mathbf{W}$ be denoted by
$\mathbf{Q}\mathbf{W}_{*,i}$. Then,

$\left\langle\mathbf{Q}\mathbf{W}_{*,i},\mathbf{Q}\mathbf{W}_{*,i}\right\rangle=\mathbf{W}_{*,i}^{\top}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{W}_{*,i}=1$

and similarly for any $j\neq i$,

$\left\langle\mathbf{Q}\mathbf{W}_{*,i},\mathbf{Q}\mathbf{W}_{*,j}\right\rangle=\mathbf{W}_{*,i}^{\top}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{W}_{*,j}=0$

where we use that $\mathbf{Q}^{\top}\mathbf{Q}=\mathbf{I}$ and the columns of
$\mathbf{W}$ are orthonormal, which holds by definition. Therefore,
$z_{i}=\mathbf{Q}\tilde{z}_{i}$ is the $i$-th singular vector of
$\mathbf{Q}\mathbf{Q}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{Q}\mathbf{Q}^{\top}$.
Let $\tilde{\mathbf{Z}}$ be the matrix obtained by stacking the vectors
$\tilde{z}_{i}$ together. Then, we have

$\begin{split}\sigma_{i}(\mathbf{A}\mathbf{Z})^{2}=\sigma_{i}^{2}(\mathbf{A}\mathbf{Q}\tilde{\mathbf{Z}})&=\sigma_{i}^{2}(\mathbf{A}\mathbf{Q})\\\
&=\sigma_{i}^{2}(\mathbf{A}\mathbf{Q}\mathbf{Q}^{\top})\\\
&=z_{i}^{T}\mathbf{Q}\mathbf{Q}^{\top}\mathbf{A}^{\top}\mathbf{A}\mathbf{Q}\mathbf{Q}^{\top}z_{i}\\\
&=z_{i}^{\top}\mathbf{A}^{\top}\mathbf{A}z_{i}\end{split}$ (5.5)

where the first equality follows from the definition of $\tilde{\mathbf{Z}}$,
the second follows from observing that $\tilde{\mathbf{Z}}$ are the singular
vectors of $\mathbf{A}\mathbf{Q}$ as shown above, the third follows from
$\mathbf{Q}^{\top}$ having orthonormal rows, the fourth from $z_{i}$ being the
$i$-th singular vector of $\mathbf{A}\mathbf{Q}\mathbf{Q}^{\top}$ and the last
from observing that $z_{i}$ is in the column span of $\mathbf{Q}$ and thus
$\mathbf{Q}\mathbf{Q}^{\top}z_{i}=z_{i}$. This concludes the proof. ∎

Next, we show a lemma relating a high-accuracy per vector guarantee to cost on
the residual subspace.

######  Lemma 5.8 (High-Accuracy Per-Vector Guarantee to Residual Cost).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, integer $p\geq 1$,
$k\in[d]$, $\ell\in[k]$ and orthonormal vectors $\\{w_{i}\\}_{i\in[\ell]}$
such that
$\left\|\mathbf{A}^{\top}w_{i}\right\|_{2}^{2}\geq\sigma_{i}^{2}-\textrm{poly}\left(\epsilon/d\right)\sigma_{k+1}^{2}$
and $(\sigma_{\ell}-\sigma_{\ell+1})/\sigma_{\ell}\geq\epsilon/d$. Let
$\mathbf{W}$ be the matrix formed by stacking together the $w_{i}$’s as
columns. Then,

$\left\|\mathbf{A}_{\ell}^{\top}\left(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top}\right)\right\|^{2}_{F}\leq\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2},$

where $\mathbf{A}_{\ell}$ is the matrix obtained by truncating all but the top
$\ell$ singular values of $\mathbf{A}$.

###### Proof.

By Pythagorean Theorem,

$\begin{split}\left\|\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{\top}\right\|_{F}^{2}&=\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\mathbf{W}^{\top}\right\|_{F}^{2}+\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)^{\top}\mathbf{W}\mathbf{W}^{\top}\right\|_{F}^{2}\\\
&\leq\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\mathbf{W}^{\top}\right\|_{F}^{2}+\sigma_{\ell+1}^{2}\left(\left\|\mathbf{W}\right\|^{2}_{F}-\left\|\mathbf{U}_{\ell}^{\top}\mathbf{W}\right\|^{2}_{F}\right),\end{split}$
(5.6)

where $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$. Further,

$\begin{split}\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\mathbf{W}^{\top}\right\|_{F}^{2}=\left\|\mathbf{\Sigma}_{\ell}\mathbf{U}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}\leq\sum_{i\in[\ell]}\sigma_{i}^{2}-\sigma_{\ell}^{2}\left(\ell-\left\|\mathbf{U}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}\right),\end{split}$
(5.7)

where the last inequality is obtained by making the Euclidean norm of all of
the $\mathbf{U}_{\ell}\mathbf{W}$’s in $[\ell-1]$to be $1$, and the $\ell$-th
row to be
$\left\|(\mathbf{U}_{\ell}^{\top}\mathbf{W})_{\ell,*}\right\|^{2}_{2}=\left\|\mathbf{U}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}-(\ell-1)$.
Rearranging Equation (5.7), we have

$\ell-\left\|\mathbf{U}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}\leq\frac{\left\|\mathbf{A}_{\ell}^{\top}\right\|_{F}^{2}-\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}}{\sigma_{\ell}^{2}}.$
(5.8)

Now, observe $\left\|\mathbf{W}\right\|_{F}^{2}=\ell$, and substituting (5.8)
back into (5.7),

$\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}\leq\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}+\frac{\sigma_{\ell+1}^{2}}{\sigma_{\ell}^{2}}\left(\left\|\mathbf{A}_{\ell}^{\top}\right\|_{F}^{2}-\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}\right).$
(5.9)

Next, we can use the guarantee’s on the $w_{i}$ to obtain a lower bound on
$\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}$ as follows:

$\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}=\sum_{i\in[\ell]}\left\|\mathbf{A}^{\top}w_{i}\right\|_{2}^{2}\geq\left\|\mathbf{A}_{\ell}^{\top}\right\|_{F}^{2}-\ell\cdot\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2},$
(5.10)

Combining equations (5.9) and (5.10), we have

$\begin{split}\left\|\mathbf{A}_{\ell}^{\top}-\mathbf{A}_{\ell}^{\top}\mathbf{W}\mathbf{W}^{\top}\right\|_{F}^{2}&=\left(\left\|\mathbf{A}_{\ell}^{\top}\right\|_{F}^{2}-\left\|\mathbf{A}_{\ell}^{\top}\mathbf{W}\right\|_{F}^{2}\right)\\\
&\leq\frac{\ell\textrm{poly}(\epsilon/d)\sigma_{k+1}^{2}}{1-(\sigma_{\ell+1}/\sigma_{\ell})^{2}}\\\
&\leq\textrm{poly}(\epsilon/d)\sigma_{k+1}^{2},\end{split}$ (5.11)

which concludes the proof. ∎

Next, we need a lemma relating the Schatten-$p$ norm of $\mathbf{A}\mathbf{Z}$
to that of $\mathbf{W}^{\top}\mathbf{A}$, where $\mathbf{Z}$ is an arbitrary
orthonormal basis and $\mathbf{W}$ is an orthonormal basis for
$\mathbf{A}\mathbf{Z}$.

######  Lemma 5.9.

Given a full-rank $n\times d$ matrix $\mathbf{A}$, let $\mathbf{W}$ be a
$n\times k$ matrix with orthonormal columns. Further, let $\mathbf{Z}$ be an
$d\times k$ matrix with orthonormal columns such that $\mathbf{Z}$ is a basis
for $\mathbf{A}^{\top}\mathbf{W}$. Then, for all $i\in[k]$,

$\sigma_{i}\left(\mathbf{A}\mathbf{Z}\right)^{p}\geq\sigma_{i}\left(\mathbf{A}^{\top}\mathbf{W}\right)^{p}$

###### Proof.

We use the following fact that for two matrices $\mathbf{A}$ and $\mathbf{B}$,
we have that for all $i$,
$\sigma_{i}(\mathbf{A}\cdot\mathbf{B})\leq\sigma_{i}(\mathbf{A})\cdot\sigma_{1}(\mathbf{B})$;
see, e.g., (2) in [LC15] and references [33-36] therein.

Using this fact, we have

$\sigma_{i}(\mathbf{A}^{\top}\mathbf{W})=\sigma_{i}(\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{T})=\sigma_{i}(\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{T})\leq\sigma_{i}(\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}^{\top})\cdot\sigma_{1}(\mathbf{W}\mathbf{W}^{T})=\sigma_{i}(\mathbf{Z}\mathbf{Z}^{T}\mathbf{A}^{\top})=\sigma_{i}(\mathbf{Z}^{T}\mathbf{A}),$

where we have used that $\sigma_{1}(\mathbf{W}\mathbf{W}^{T})=1$ since
$\mathbf{W}\mathbf{W}^{T}$ is a projection matrix, and the fact that
$\mathbf{Z}\mathbf{Z}^{T}$ is a basis for the column span of
$\mathbf{A}^{\top}\mathbf{W}$. Raising both sides to the $p$-th power
establishes the lemma.

∎

Finally, we can combine the two aforementioned lemmas to obtain the following
corollary:

######  Corollary 5.10 (Changing Basis for high-accuracy vectors).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, integer $p\geq 1$,
$k\in[d]$, $\ell\in[k]$ and orthonormal vectors $\\{w_{i}\\}_{i\in[\ell]}$
such that
$\left\|\mathbf{A}^{\top}w_{i}\right\|_{2}^{2}\geq\sigma_{i}^{2}-\textrm{poly}\left(\epsilon/d\right)\sigma_{k+1}^{2}$
and $(\sigma_{\ell}-\sigma_{\ell+1})/\sigma_{\ell}\geq\epsilon/d$. Let
$\mathbf{W}$ be the matrix formed by stacking together the $w_{i}$’s as
columns and let $\mathbf{Z}$ be an orthonormal basis for
$\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{\top}$. Then,

$\left\|\mathbf{A}_{\ell}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|^{2}_{F}\leq\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2}.$

###### Proof.

We closely follow the proof in Lemma 5.8. By Pythagorean Theorem,

$\begin{split}\left\|\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{F}^{2}&=\left\|\mathbf{A}_{\ell}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{F}^{2}+\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\mathbf{Z}\mathbf{Z}^{\top}\right\|_{F}^{2}\\\
&\leq\left\|\mathbf{A}_{\ell}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{F}^{2}+\sigma_{\ell+1}^{2}\left(\left\|\mathbf{Z}\right\|^{2}_{F}-\left\|\mathbf{V}_{\ell}\mathbf{Z}\right\|^{2}_{F}\right),\end{split}$
(5.12)

where $\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$. Further,

$\begin{split}\left\|\mathbf{A}_{\ell}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{F}^{2}=\left\|\mathbf{\Sigma}_{\ell}\mathbf{V}_{\ell}^{\top}\mathbf{Z}\right\|_{F}^{2}\leq\sum_{i\in[\ell]}\sigma_{i}^{2}-\sigma_{\ell}^{2}\left(\ell-\left\|\mathbf{V}_{\ell}^{\top}\mathbf{Z}\right\|_{F}^{2}\right),\end{split}$
(5.13)

where the last inequality is obtained by making the Euclidean norm of all of
the $\mathbf{V}_{\ell}^{\top}\mathbf{Z}$’s in $[\ell-1]$to be $1$, and the
$\ell$-th row to be
$\left\|(\mathbf{V}_{\ell}^{\top}\mathbf{Z})_{\ell,*}\right\|^{2}_{2}=\left\|\mathbf{V}_{\ell}^{\top}\mathbf{Z}\right\|_{F}^{2}-(\ell-1)$.
Rearranging Equation (5.13), we have

$\ell-\left\|\mathbf{V}_{\ell}^{\top}\mathbf{Z}\right\|_{F}^{2}\leq\frac{\left\|\mathbf{A}_{\ell}\right\|_{F}^{2}-\left\|\mathbf{A}_{\ell}\mathbf{Z}\right\|_{F}^{2}}{\sigma_{\ell}^{2}}.$
(5.14)

Now, observe $\left\|\mathbf{Z}\right\|_{F}^{2}=\ell$, and substituting (5.14)
back into (5.13),

$\left\|\mathbf{A}\mathbf{Z}\right\|_{F}^{2}\leq\left\|\mathbf{A}_{\ell}\mathbf{Z}\right\|_{F}^{2}+\frac{\sigma_{\ell+1}^{2}}{\sigma_{\ell}^{2}}\left(\left\|\mathbf{A}_{\ell}\right\|_{F}^{2}-\left\|\mathbf{A}_{\ell}\mathbf{Z}\right\|_{F}^{2}\right).$
(5.15)

Next, we can use the guarantee’s on the $w_{i}$ to obtain a lower bound on
$\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}$ as follows:

$\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}=\sum_{i\in[\ell]}\left\|\mathbf{A}^{\top}w_{i}\right\|_{2}^{2}\geq\left\|\mathbf{A}_{\ell}^{\top}\right\|_{F}^{2}-\ell\cdot\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2},$
(5.16)

Observe
$\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}=\sum_{i\in[\ell]}\sigma_{i}^{2}(\mathbf{A}^{\top}\mathbf{W})$.
By Lemma 5.9, we know that for all $i$,
$\sigma_{i}^{2}(\mathbf{A}^{\top}\mathbf{W})\leq\sigma_{i}^{2}(\mathbf{A}\mathbf{Z})$.
Therefore, we can restate the above equation as follows:

$\left\|\mathbf{A}\mathbf{Z}\right\|_{F}^{2}\geq\left\|\mathbf{A}^{\top}\mathbf{W}\right\|_{F}^{2}\geq\left\|\mathbf{A}_{\ell}\right\|_{F}^{2}-\ell\cdot\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2},$
(5.17)

Combining equations (5.15) and (5.17), we have

$\begin{split}\left\|\mathbf{A}_{\ell}-\mathbf{A}_{\ell}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{F}^{2}&=\left(\left\|\mathbf{A}_{\ell}\right\|_{F}^{2}-\left\|\mathbf{A}_{\ell}\mathbf{Z}\right\|_{F}^{2}\right)\\\
&\leq\frac{\ell\textrm{poly}(\epsilon/d)\sigma_{k+1}^{2}}{1-(\sigma_{\ell+1}/\sigma_{\ell})^{2}}\\\
&\leq\textrm{poly}(\epsilon/d)\sigma_{k+1}^{2},\end{split}$ (5.18)

which concludes the proof. ∎

We now have all the ingredients we need to complete the proof of Theorem 5.1.

###### Proof of Theorem 5.1.

Observe, using Lemma 5.2 with probability at least $97/100$, Step 1 of
Algorithm 5.4 outputs $\hat{\sigma}_{i}$’s such that for all $i\in[k+1]$,
$\hat{\sigma}_{i}^{2}=\left(1\pm 0.1/p\right)\sigma_{i}^{2}$ and
$\hat{\sigma}_{s}^{2}=\left(1\pm 0.1/p\right)\sigma_{s}^{2}$, for
$s=\mathcal{O}\left(kp^{-1/3}/\epsilon^{1/3}\right)$. Condition on this event.

At a high level, we proceed via a case analysis: either the Schatten-$p$ norm
of the tail is large compared to the $(k+1)$-st singular value, and we don’t
require a highly accurate solution, or the Schatten-$p$ norm of the tail is
small, and increasing the block size induces a gap. We formalize this
intuition into a proof.

#### No large singular values.

Let us first consider the case where
$\sigma_{1}<\left(1+1/p\right)\sigma_{k+1}$. We yet again split into cases,
and consider the case where the Schatten-$p$ norm of the tail is small, i.e.
$\left\|\mathbf{A}-\mathbf{A}_{k}\right\|^{p}_{\mathcal{S}_{p}}\leq\frac{k}{p^{1/3}\epsilon^{1/3}}\cdot\sigma^{p}_{k+1}$.
Observe, for any $t\in[1,d-k-1]$,

$\frac{k}{p^{1/3}\epsilon^{1/3}}\cdot\sigma^{p}_{k+1}\geq\left\|\mathbf{A}-\mathbf{A}_{k}\right\|^{p}_{\mathcal{S}_{p}}\geq\sum_{i=k+1}^{k+1+t}\sigma_{i}^{p}\geq
t\sigma_{k+1+t}^{p}.$ (5.19)

Then, setting
$t=\frac{\left(1+1/p\right)^{p}k}{\epsilon^{1/3}p^{1/3}}=\Theta\left(\frac{k}{\epsilon^{1/3}p^{1/3}}\right)$,
we have $\sigma_{k+1+t}\leq\sigma_{k+1}/\left(1+1/p\right)$. It suffices to
show that we can detect this gap for some $s\geq k+1+t$. Recall, we know that
$\hat{\sigma}_{k+1}=(1\pm 0.1/p)\sigma_{k+1}$ and $\hat{\sigma}_{s}=(1\pm
0.1/p)\sigma_{s}$. Then, we have

$\hat{\sigma}_{s}\leq\left(1+\frac{0.1}{p}\right)\sigma_{s}\leq\left(1+\frac{0.1}{p}\right)\sigma_{k+1+t}\leq\left(1+\frac{0.1}{p}\right)\cdot\left(\frac{1}{1+1/p}\right)\sigma_{k+1}\leq\frac{1}{\left(1+\frac{0.5}{p}\right)}\hat{\sigma}_{k+1}.$
(5.20)

Therefore, Algorithm 5.4 outputs $\mathbf{Z}$, an orthonormal basis for
$\mathbf{A}^{\top}\mathbf{W}_{2}$, where $\mathbf{W}_{2}$ is obtained by
running Algorithm 5.6 on $\mathbf{A}^{\top}$, initialized with a block size of
$\Theta\left(\frac{k}{\epsilon^{1/3}p^{1/3}}\right)$ and run for
$\mathcal{O}\left(\log(d/\epsilon)\sqrt{p}\right)$ iterations. Observe, since
$\sigma_{k+1+t}\leq\sigma_{k+1}/\left(1+1/p\right)$, this suffices to
demonstrate a gap that depends on $p$ as follows:
$\frac{\sigma_{k}}{\sigma_{k}-\sigma_{k+t+1}}\leq p$. Recall, we account for
this gap by running $\mathcal{O}\left(\log(d)\sqrt{p}\right)$ iterations.
Using the gap dependent analysis (Lemma 5.3), we can conclude that with
probability at least $99/100$, for all $i\in[k]$,

$\left\|\mathbf{A}^{\top}(\mathbf{W}_{2})_{*,i}\right\|_{2}^{2}\geq\sigma_{i}^{2}-\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma^{2}_{k+1}.$
(5.21)

Then, applying Lemma 5.7 with $\mathbf{W}_{2}\mathbf{W}_{2}^{\top}$ satisfying
the guarantee in (5.21), we have

$\begin{split}\left\|\mathbf{A}^{\top}\mathbf{W}_{2}\mathbf{W}_{2}^{\top}\right\|_{\mathcal{S}_{p}}^{p}&\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\textrm{poly}\left(\frac{\epsilon}{d}\right)\sum_{i\in[k]}\sigma_{k+1}^{2}\sigma_{i}^{p-2}\\\
&\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{p}.\end{split}$
(5.22)

where the last inequality uses that $\sigma_{1}<(1+1/p)\sigma_{k+1}$ and
$(1+1/p)^{p-2}=\mathcal{O}\left(1\right)$. Next, we use Lemma 5.5 to relate
$\left\|\mathbf{A}^{\top}\mathbf{W}_{2}\mathbf{W}_{2}^{\top}\right\|_{\mathcal{S}_{p}}^{p}$
to
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$,
where $\mathbf{Z}$ is an orthonormal basis for
$\mathbf{A}^{\top}\mathbf{W}_{2}\mathbf{W}_{2}^{\top}$ as output by the
algorithm. Setting $\mathbf{Q}=\mathbf{Z}\mathbf{Z}^{\top}$ and
$\mathbf{P}=\mathbf{W}_{2}\mathbf{W}_{2}^{\top}$, we observe that
$\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{A}^{\top}\mathbf{W}_{2}\mathbf{W}_{2}^{\top}\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{W}_{2}\mathbf{W}_{2}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}$
and
$\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$.
Then, invoking Lemma 5.5 and plugging in Equation (5.22), we have

$\begin{split}\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}&\leq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{A}^{\top}\mathbf{W}_{2}\mathbf{W}_{2}^{\top}\right\|_{\mathcal{S}_{p}}^{p}\\\
&\leq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}+\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{p}\\\
&\leq\left(1+\textrm{poly}\left(\frac{\epsilon}{d}\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p},\end{split}$
(5.23)

which concludes the analysis in this case.

As shown in Equation 5.20, we can detect a gap between $\sigma_{k+1+t}$ and
$\sigma_{k+1}$ by comparing $\hat{\sigma}_{s}$ and $\hat{\sigma}_{k+1}$. When
5.20 does not hold, we know that
$\hat{\sigma}_{s}\geq\left(1+0.5/p\right)\hat{\sigma}_{k+1}$ and Algorithm 5.4
outputs $\mathbf{Z}$, an orthonormal basis for
$\mathbf{A}^{\top}\mathbf{W}_{1}\mathbf{W}_{1}^{\top}$. Since we have $(1\pm
0.1/p)$-approximate estimates to these quantities, we can conclude that
$\sigma_{s}\geq\left(1+0.1/p\right)\sigma_{k+1}$. Then, we have

$\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}\geq
s\cdot\sigma_{s}^{p}=\Omega\left(\frac{k}{\epsilon^{1/3}p^{1/3}}\right)\sigma_{k+1}^{p}.$
(5.24)

It therefore remains to consider the case where
$\left\|\mathbf{A}-\mathbf{A}_{k}\right\|^{p}_{\mathcal{S}_{p}}>\frac{ck}{p^{1/3}\epsilon^{1/3}}\cdot\sigma^{p}_{k+1}$,
for a fixed universal constant $c$. Here, we note that the tail is large
enough that an additive error of
$\mathcal{O}\left(\epsilon^{2/3}p^{1/3}\right)\sigma_{k+1}^{2}$ on each of the
top-$k$ singular values suffices. Formally, it follows from Lemma 5.2 (setting
$\gamma=\epsilon^{2/3}p^{-1/3}$, and invoking it for $\mathbf{A}^{\top}$) that
initializing Algorithm 5.6 with block size $k$ and running for
$\mathcal{O}\left(\log(d/\epsilon)p^{1/6}/\epsilon^{1/3}\right)$ iterations
suffices to output a $n\times k$ matrix $\mathbf{W}_{1}$ such that with
probability at least $99/100$, for all $i\in[k]$,

$\left\|\mathbf{A}^{\top}\left(\mathbf{W}_{1}\right)_{*,i}\right\|_{2}^{2}\geq\sigma_{i}^{2}-\epsilon^{2/3}p^{-1/3}\sigma_{k+1}^{2}.$
(5.25)

Then, invoking Lemma 5.7 with $\mathbf{A}^{\top}$ and $\mathbf{W}_{1}$ as
defined above, we have

$\begin{split}\left\|\mathbf{A}^{\top}\mathbf{W}_{1}\mathbf{W}_{1}^{\top}\right\|_{\mathcal{S}_{p}}^{p}&=\left\|\mathbf{W}_{1}\mathbf{W}_{1}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}\\\
&\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\sum_{i\in[k]}\mathcal{O}\left(\epsilon^{2/3}p^{-1/3}p\right)\sigma_{k+1}^{2}\sigma_{i}^{p-2}\\\
&\geq\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}-\mathcal{O}\left(k\epsilon^{2/3}p^{2/3}\right)\sigma_{k+1}^{p}\end{split}$
(5.26)

where the last inequality uses that $\sigma_{1}<(1+1/p)\sigma_{k+1}$ and
$(1+1/p)^{p}=\mathcal{O}\left(1\right)$. Recall, in this case, Algorithm 5.4
outputs $\mathbf{Z}\mathbf{Z}^{\top}$ where $\mathbf{Z}$ is an orthonormal
basis for $\mathbf{A}^{\top}\mathbf{W}_{1}\mathbf{W}^{\top}_{1}$. Next, we
invoke Lemma 5.5 to relate
$\left\|\mathbf{A}^{\top}\mathbf{W}_{1}\mathbf{W}_{1}^{\top}\right\|_{\mathcal{S}_{p}}^{p}$
to
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$.
Setting $\mathbf{Q}=\mathbf{Z}\mathbf{Z}^{\top}$ and
$\mathbf{P}=\mathbf{W}_{1}\mathbf{W}^{\top}_{1}$, we observe that
$\left\|\mathbf{P}\mathbf{A}\mathbf{Q}\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{W}_{1}\mathbf{W}_{1}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}$
and
$\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}$.
Then, invoking Lemma 5.5 and plugging in Equation (5.26), we have

$\begin{split}\left\|\left(\mathbf{I}-\mathbf{P}\right)\mathbf{A}\left(\mathbf{I}-\mathbf{Q}\right)\right\|_{\mathcal{S}_{p}}^{p}=\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}&\leq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{W}_{1}\mathbf{W}_{1}^{\top}\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}\\\
&\leq\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}^{p}-\left\|\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}+\mathcal{O}\left(k\epsilon^{2/3}p^{2/3}\right)\sigma_{k+1}^{p}\\\
&\leq\left(1+\mathcal{O}\left(p\epsilon\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p},\end{split}$
(5.27)

where the last inequality follows from our assumption on the Schatten-$p$ norm
of the tail, given the case we are in. Taking the $(1/p)$-th root, and
recalling that $\epsilon<1/2$, we obtain

$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq\left(1+\mathcal{O}\left(\epsilon\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{p},$
(5.28)

which concludes the case where $\ell=0$.

#### Large Singular Values.

Next, we consider the case where $\sigma_{1}>\left(1+1/p\right)\sigma_{k+1}$.
Then, let $\ell\in[k]$ be the largest integer such that
$\sigma_{\ell}>\left(1+0.5/p\right)\sigma_{k+1}$ and
$\sigma_{\ell+1}<(1-\epsilon/d)\sigma_{\ell}$. Observe, such an $\ell$ is
guaranteed to exist.

We then note that in all settings Algorithm 5.4 runs
$\Omega(\log(d/\epsilon)\sqrt{p}$ iterations on $\mathbf{A}^{\top}$ and since
exists a gap of size $p$ between $\sigma_{\ell}$ and $\sigma_{k+1}$ it follows
from Lemma 5.3 that running Block Krylov Iteration that Algorithm 5.4 always
outputs an orthonormal matrix $\mathbf{W}$ s.t. for all $i\in[\ell]$,

$\left\|\mathbf{A}^{\top}\mathbf{W}_{*,i}\right\|^{2}\geq\sigma_{i}^{2}-\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2}.$
(5.29)

Further, for all $i\in[\ell+1,k]$, we have

$\left\|\mathbf{A}^{\top}\mathbf{W}_{*,i}\right\|^{2}\geq\sigma_{i}^{2}-\gamma_{i}\sigma_{k+1}^{2},$
(5.30)

where $\gamma_{i}$ is determined by whether $\mathbf{W}=\mathbf{W}_{1}$ or
$\mathbf{W}=\mathbf{W}_{2}$, as we discuss later.

We note that we cannot simply take $p/2$-th powers here (for large $p$) as
this would introduce cross terms that scale proportional to
$\sigma_{i}(\mathbf{A})$, which can be significantly larger than
$\sigma_{k+1}(\mathbf{A})$. Instead, we require a finer analysis by splitting
$\mathbf{A}$ into a head and tail term. Further, we let $\mathbf{Z}$ be an
orthonormal basis for $\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{\top}$.

We are now ready to bound
$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}$.
By the triangle inequality,

$\begin{split}\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}&\leq\left\|\mathbf{A}_{\ell}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}+\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\\\
\end{split}$ (5.31)

By Corollary 5.10, we know that
$\left\|\mathbf{A}_{\ell}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{F}^{2}\leq\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}^{2}$,
and since all Schatten norms are within a $\sqrt{d}$ factor of the Frobenius
norm, we have

$\left\|\mathbf{A}_{\ell}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}.$

Substituting this back into Equation (5.31), we have

$\left\|\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}\leq\mathcal{O}\left(\frac{\epsilon}{d}\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}+\underbrace{\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}}_{\ref{eqn:bounding-
all-but-one-cost}.1}.$ (5.32)

It remains to bound term 5.32.1 above. By triangle inequality, we have

$\begin{split}&\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}^{p}\\\
&\leq\left(\left\|(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top})\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}+\left\|\mathbf{W}\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}\right)^{p}\end{split}$
(5.33)

We bound the two terms on the RHS independently. To upper bound
$\left\|\mathbf{W}\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}^{p}$
we use the relation between Frobenius and Schatten norms, and recall that by
definition,
$\mathbf{W}\mathbf{W}^{\top}\mathbf{A}\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)=\mathbf{W}\mathbf{W}^{\top}\mathbf{A}-\mathbf{W}\mathbf{W}^{\top}\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}=0$,
and thus

$\begin{split}\left\|\mathbf{W}\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}&\leq\sqrt{k}\left\|\mathbf{W}\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{F}\\\
&=\sqrt{k}\left\|\mathbf{W}\mathbf{W}^{\top}\mathbf{A}_{\ell}(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{F}\\\
&\leq\sqrt{k}\left\|\mathbf{A}_{\ell}(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{F}\\\
&\leq\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1},\end{split}$
(5.34)

where the last inequality follows from Corollary 5.10. Therefore, combining
the above, we have

$\begin{split}\left\|\mathbf{W}\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}\leq\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1},\end{split}$
(5.35)

It remains to upper bound
$\left\|(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top})\left(\mathbf{A}-\mathbf{A}_{\ell}\right)(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top})\right\|_{\mathcal{S}_{p}}^{p}$.
Recall, $\mathbf{W}$ is the orthonormal basis output by Block Krylov run on
$\mathbf{A}^{\top}$, and in Algorithm 5.4 $\mathbf{W}$ is either
$\mathbf{W}_{1}$ or $\mathbf{W}_{2}$. Let $\mathbf{Z}$ is a basis for
$\mathbf{A}^{\top}\mathbf{W}\mathbf{W}^{\top}$. Then, applying Lemma 5.5 with
$\mathbf{Q}=\mathbf{Z}\mathbf{Z}^{\top}$ and
$\mathbf{P}=\mathbf{W}\mathbf{W}^{\top}$, we have

$\begin{split}\left\|\left(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top}\right)\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|^{p}_{\mathcal{S}_{p}}&\leq\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\right\|^{p}_{\mathcal{S}_{p}}-\left\|\mathbf{W}\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|^{p}_{\mathcal{S}_{p}}\\\
&=\sum_{j\in[\ell+1,d]}\sigma_{j}^{p}-\sum_{j\in[k]}\sigma_{j}^{p}\left(\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\mathbf{Z}\right)\end{split}$
(5.36)

Next, we show that for all $j\in[k]$,
$\sigma_{j}\left(\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\mathbf{Z}\right)\geq\sigma_{j+\ell}\left(\mathbf{W}^{\top}\mathbf{A}\right)$.
Here, we invoke Fact 4.5 for
$\mathbf{X}=\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\mathbf{Z}$
and $\mathbf{Y}=\mathbf{W}^{\top}\mathbf{A}_{\ell}\mathbf{Z}$, with $i=j$ and
$j=\ell$. Note, the precondition on the indices $i,j$ in Fact 4.5 is satisfied
since $\mathbf{X},\mathbf{Y}$ are $n\times k$ matrices, and $j\in[k]$ and
$\ell<k$. Then, we have

$\begin{split}\sigma_{j+\ell}\left(\mathbf{W}^{\top}\mathbf{A}\mathbf{Z}\right)&=\sigma_{j+\ell}\left(\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\mathbf{Z}\right)+\mathbf{W}^{\top}\mathbf{A}_{\ell}\mathbf{Z}\right)\\\
&\leq\sigma_{j}\left(\mathbf{W}^{\top}\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\mathbf{Z}\right)+\sigma_{\ell+1}\left(\mathbf{W}^{\top}\mathbf{A}_{\ell}\mathbf{Z}\right),\end{split}$

but $\mathbf{W}^{\top}\mathbf{A}_{\ell}\mathbf{Z}$ is a rank $\leq\ell$
matrix, and thus
$\sigma_{\ell+1}\left(\mathbf{W}^{\top}\mathbf{A}_{\ell}\mathbf{Z}\right)=0$.
Therefore, we can conclude,

$\left\|\left(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top}\right)\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq\sum_{j\in[\ell+1,d]}\sigma^{p}_{j}-\sum_{j\in[\ell+1,k+\ell]}\sigma_{j}^{p}\left(\mathbf{W}^{\top}\mathbf{A}\mathbf{Z}\right).$
(5.37)

However, now we observe that
$\mathbf{W}\mathbf{W}^{\top}\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}=\mathbf{W}\mathbf{W}^{\top}\mathbf{A}$,
and thus
$\sigma_{j}\left(\mathbf{W}^{\top}\mathbf{A}\mathbf{Z}\right)=\sigma_{j}\left(\mathbf{W}\mathbf{A}\right)$.
Recall, for all $j\in[k]$, it follows from Equation (5.4) in the proof of
Lemma 5.7 that
$\sigma_{j}^{p}(\mathbf{W}^{\top}\mathbf{A})=\sigma_{j}^{p}(\mathbf{A}^{\top}\mathbf{W})\geq\sigma_{j}^{p}(\mathbf{A})-\mathcal{O}\left(\gamma_{j}p\right)\sigma_{k+1}^{2}\sigma_{j}^{p-2}$.
Further, by definition, for $j\in[\ell+1,k]$,
$\sigma_{j}\leq\left(1+1/p\right)\sigma_{k+1}$ and thus, for all
$j\in[\ell+1,k]$,

$\begin{split}\sigma_{j}^{p}(\mathbf{W}^{\top}\mathbf{A})&\geq\sigma_{j}^{p}-\mathcal{O}\left(\gamma_{j}p\left(1+1/p\right)^{p-2}\right)\sigma_{k+1}^{p}\\\
&\geq\sigma_{j}^{p}-\mathcal{O}\left(\gamma_{j}p\right)\sigma_{k+1}^{p}.\end{split}$
(5.38)

Recall, we can correctly determine whether
$\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}\leq\frac{k}{p^{1/3}\epsilon^{1/3}}\sigma_{k+1}^{p}$
or not, up to error in estimating the singular values as shown earlier, in the
analysis for the case where there are no large singular values. In particular,
let us first consider the case where this is true. Repeating the argument in
equations (5.19), (5.20), we can conclude that Algorithm 5.4 outputs
$\mathbf{W}=\mathbf{W}_{2}$ and thus for all $j\in[\ell+1,k]$,
$\gamma_{j}=\textrm{poly}(\epsilon/d)$. Therefore, substituting this back into
Equation (5.37), we have

$\begin{split}\left\|\left(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top}\right)\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}&\leq\sum_{j\in[\ell+1,d]}\sigma^{p}_{j}-\left(\sum_{j\in[\ell+1,k]}\sigma^{p}_{j}-\mathcal{O}\left(\gamma_{j}p\right)\sigma_{k+1}^{p}\right)\\\
&\leq\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}+\sum_{j\in[\ell+1,k]}\mathcal{O}\left(\gamma_{j}p\right)\sigma_{k+1}^{p}\\\
&\leq\left(1+\textrm{poly}\left(\frac{\epsilon}{d}\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p},\end{split}$
(5.39)

concluding the analysis in this case.

Next, consider the case where
$\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}>\frac{k}{p^{1/3}\epsilon^{1/3}}\sigma_{k+1}^{p}$.
Then, repeating the analysis in Equation (5.24), we know that Algorithm 5.4
outputs $\mathbf{W}=\mathbf{W}_{2}$, and thus for all $j\in[\ell+1,k]$,
$\gamma_{j}=\epsilon^{2/3}p^{-1/3}$ as shown in Equation (5.25). Again,
substituting this back into Equation (5.37), we have

$\begin{split}\left\|\left(\mathbf{I}-\mathbf{W}\mathbf{W}^{\top}\right)\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}&\leq\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}+\sum_{j\in[\ell+1,k]}\mathcal{O}\left(\gamma_{j}p\right)\sigma_{k+1}^{p}\\\
&\leq\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p}+\mathcal{O}\left(kp\cdot\epsilon^{2/3}p^{-1/3}\right)\sigma_{k+1}^{p}\\\
&\leq\left(1+\mathcal{O}\left(\epsilon
p\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}^{p},\end{split}$
(5.40)

where the last inequality follows from our assumption. Taking the $(1/p)$-th
root and substituting equations (5.35) and (5.40) back into (5.33), we can
conclude

$\left\|\left(\mathbf{A}-\mathbf{A}_{\ell}\right)\left(\mathbf{I}-\mathbf{Z}\mathbf{Z}^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq\left(1+\mathcal{O}\left(\epsilon\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}}+\textrm{poly}\left(\frac{\epsilon}{d}\right)\sigma_{k+1}\leq\left(1+\mathcal{O}\left(\epsilon\right)\right)\left\|\mathbf{A}-\mathbf{A}_{k}\right\|_{\mathcal{S}_{p}},$

which when substituted into Equation (5.31) concludes the analysis.

Next, we analyze the running time and matrix-vector products. Running
Algorithm 5.6 with block size $k$ for
$q=\mathcal{O}\left(\log(d)p^{1/6}/\epsilon^{1/3}\right)$ iterations requires
$\mathcal{O}\left(\frac{\operatorname{\mathtt{nnz}}(\mathbf{A})kp^{1/6}\log(d)}{\epsilon^{1/3}}\right)$
time and $\mathcal{O}\left(\frac{kp^{1/6}\log(d)}{\epsilon^{1/3}}\right)$
matrix-vector products. Similarly, running with block size
$\mathcal{O}\left(k/\left(\epsilon p\right)^{1/3}\right)$ for
$q=\mathcal{O}\left(\log(d/\epsilon)\sqrt{p}\right)$ iterations requires
$\mathcal{O}\left(\frac{\operatorname{\mathtt{nnz}}(\mathbf{A})kp^{1/6}\log(d/\epsilon)}{\epsilon^{1/3}}\right)$
time and $\mathcal{O}\left(\frac{kp^{1/6}\log(d)}{\epsilon^{1/3}}\right)$
matrix-vector products. Finally, we observe that to obtain a
$\left(1+1/p\right)$-approximation to $\sigma_{1}$ and $\sigma_{k+1}$, we need
$\mathcal{O}\left(\log(d)\sqrt{p}\right)$ iterations with blocksize $k+1$ and
this requires $\mathcal{O}\left(\log(d)\sqrt{p}k\right)$ matrix-vector
products. Note, our setting of the exponent of $p$ and $\epsilon$ was chosen
to balance the two cases, and this concludes the proof.

∎

##  6 Query Lower Bounds

Next, we show that the $\epsilon$-dependence obtained by our algorithms for
Schatten-$p$ low-rank approximation is optimal in the restricted computation
model of matrix-vector products. The matrix-vector product model is defined as
follows: given a matrix $\mathbf{A}$, our algorithm is allowed to make
adaptive matrix-vector queries to $\mathbf{A}$, where one matrix-vector query
is of the form $\mathbf{A}v$, for any $v\in\mathbb{R}^{d}$. Our lower bounds
are information-theoretic and rely on the hardness of estimating the smallest
eigenvalue of a Wishart ensemble, as established in recent work of Braverman,
Hazan, Simchowitz and Woodworth [BHSW20].

We split the lower bounds into the case of $p\in[1,2]$ and $p>2$. For
$p\in[1,2]$, we have a simple argument based on the Araki-Lieb-Thirring
inequality (Fact 4.10), whereas for $p>2$, our lower bounds require an
involved argument using a norm compression inequality for partitioned
operators (Fact 4.14).

###  6.1 Lower Bounds for $p\in[1,2]$

The main lower bound we prove in this sub-section is as follows:

######  Theorem 6.1 (Query Lower Bound for $p\in\textrm{[}1,2\textrm{]}$).

Given $\varepsilon>0$, and $p\in[1,2]$, there exists a distribution
$\mathcal{D}$ over $n\times n$ matrices such that for
$\mathbf{A}\sim\mathcal{D}$, any randomized algorithm that with probability at
least $9/10$ outputs a rank-$1$ matrix $\mathbf{B}$ such that
$\left\|\mathbf{A}-\mathbf{B}\right\|^{p}_{\mathcal{S}_{p}}\leq(1+\varepsilon)\left\|\mathbf{A}-\mathbf{A}_{1}\right\|^{p}_{\mathcal{S}_{p}}$
must make $\Omega(1/\varepsilon^{1/3})$ matrix-vector queries to $\mathbf{A}$.

We require the following theorem on the hardness of computing the minimum
eigenvalue of a Wishart Matrix, introduced recently by Braverman, Hazan,
Simchowitz and Woodworth [BHSW20]:

######  Theorem 6.2 (Computing Min Eigenvalue of Wishart, Theorem 3.1
[BHSW20]).

Given $\epsilon\in(0,1)$, there exists a function
$\mathbf{d}:(0,1)\to\mathbb{N}$ such that for all $d\geq\mathbf{d}(\epsilon)$,
the following holds. Let $\mathbf{W}\sim\textrm{Wishart}(d)$ be a Wishart
matrix and $\\{\lambda_{i}\\}_{i\in[d]}$ be the eigenvalues of $\mathbf{W}$,
in descending order. Then, there exists a universal constant $c^{*}$ such
that:

  1. 1.

Let $\zeta_{1}$ be the event that $\lambda_{d}(\mathbf{W})\leq c_{1}/d^{2}$,
$\zeta_{2}$ be the event that
$\lambda_{d-1}(\mathbf{W})-\lambda_{d}(\mathbf{W})\geq c_{2}/d^{2}$ and
$\zeta_{3}$ be the event that $\left\|\mathbf{W}\right\|_{op}\leq 5$, where
$c_{1}$ and $c_{2}$ are constants that depend only on $\epsilon$. Then,
$\Pr_{\mathbf{W}}\left[\zeta_{1}\cap\zeta_{2}\cap\zeta_{3}\right]\geq
1-\frac{c^{*}\sqrt{\epsilon}}{2}$.

  2. 2.

Any randomized algorithm that makes at most $(1-\epsilon)d$ adaptive matrix-
vector queries and outputs an estimate $\hat{\lambda}_{d}$ must satisfy

$\Pr_{\mathbf{W}}\left[\left|\hat{\lambda}_{d}-\lambda_{d}\right|\geq\frac{1}{4d^{2}}\right]\geq
c^{*}\sqrt{\epsilon}.$

We also use the following lemma from [BHSW20] bounding the minimum eigenvalue
of a Wishart ensemble:

######  Lemma 6.3 (Non-Asymptotic Spectra of Wishart Ensembles, Corollary 3.3
[BHSW20]).

Let $\mathbf{W}\sim\textsf{Wishart}(n)$ be such that
$n=\Omega(1/\varepsilon^{3})$. Then, there exists a universal constant
$c_{2}>0$ such that

$\Pr\left[\lambda_{n}\left(\mathbf{W}\right)\geq\frac{1}{n^{2}}\right]\geq
c_{2},\hskip 14.45377pt\textrm{and}\hskip
14.45377pt\Pr\left[\lambda_{n}\left(\mathbf{W}\right)<\frac{1}{2n^{2}}\right]\geq\frac{c_{2}}{2}.$

We are now ready to prove Theorem 6.1. Our high level approach is to show that
we can take any solution that is a $(1+\varepsilon)$-relative-error
Schatten-$p$ low-rank approximation to the hard instance
$\mathbf{I}-\frac{1}{5}\mathbf{W}$, where $\mathbf{W}$ is a Wishart ensemble,
and extract from it an accurate estimate of the minimum eigenvalue of
$\mathbf{W}$, thus appealing to the hardness stated in (2) of Theorem 6.2
above.

###### Proof of Theorem 6.1.

Let $n=\Theta\left(1/\epsilon^{1/3}\right)$ and let
$\mathbf{A}=\mathbf{I}-\frac{1}{5}\mathbf{W}$ be an $n\times n$ instance where
$\mathbf{W}\sim\textsf{Wishart}(n)$. Let $\zeta_{1}$ be the event that
$\left\|\mathbf{W}\right\|_{\textrm{op}}\leq 5$. It follows from Fact 4.16
that $\zeta_{1}$ holds with probability at least $99/100$, and we condition on
this event. Let $\zeta_{2}$ be the event that
$\lambda_{n}\left(\mathbf{W}\right)\geq\frac{1}{n^{2}}=\frac{\varepsilon^{2/3}}{c^{*}}$
and $\zeta_{3}$ be the event that
$\lambda_{n}\left(\mathbf{W}\right)<\frac{1}{2n^{2}}=\frac{\varepsilon^{2/3}}{2c^{*}}$.

Then, conditioning on $\zeta_{2}$, we have that

$1-\frac{1}{5}\lambda_{n}(\mathbf{W})\leq 1-\frac{\varepsilon^{2/3}}{5c^{*}}.$
(6.1)

Similarly, conditioning on $\zeta_{3}$, we have that

$1-\frac{1}{5}\lambda_{n}(\mathbf{W})\geq
1-\frac{\varepsilon^{2/3}}{10c^{*}}.$ (6.2)

We observe that for $p\in[1,2]$, using Bernoulli’s inequality (Fact 4.6) we
have

$\left(1-\frac{1}{5}\lambda_{n}(\mathbf{W})\right)^{p}\geq
1-\frac{p}{5}\lambda_{n}(\mathbf{W})$

and since $(1-x)^{p}\leq(1-x)$ for any $x\in(0,1)$, we also have that,

$\left(1-\frac{1}{5}\lambda_{n}(\mathbf{W})\right)^{p}\leq
1-\frac{1}{5}\lambda_{n}(\mathbf{W})$

Therefore, we can conclude,
$\left(1-\frac{1}{5}\lambda_{n}(\mathbf{W})\right)^{p}=1-\Theta\left(\lambda_{n}(\mathbf{W})\right)$.
Further, it follows from part (1) of Fact 4.16 that
$0\preceq\mathbf{I}-\frac{1}{5}\mathbf{W}\preceq\mathbf{I}$, and thus

$\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}=\sum_{i\in[n]}\lambda_{i}^{p}\left(\mathbf{I}-\frac{1}{5}\mathbf{W}\right)\leq\sum_{i\in[n]}\lambda_{i}\left(\mathbf{I}-\frac{1}{5}\mathbf{W}\right)\leq\mathcal{O}\left(\frac{1}{\epsilon^{1/3}}\right)$
(6.3)

where the last inequality follows from the fact that
$n=\sqrt{c^{*}}/\epsilon^{1/3}$. Let $\mathbf{A}_{1}$ denote the best rank-$1$
approximation to $\mathbf{A}$. Then, it follows from Equation (6.3) that

$\epsilon\left\|\mathbf{A}-\mathbf{A}_{1}\right\|^{p}_{\mathcal{S}_{p}}\leq\epsilon\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}\leq\mathcal{O}\left(\epsilon^{2/3}\right)$
(6.4)

Observe, any $(1+\epsilon)$-approximate relative-error Schatten-$p$ low-rank
approximation algorithm for $k=1$ outputs a matrix $vv^{\top}$ such that

$\begin{split}\left\|\mathbf{A}\left(\mathbf{I}-vv^{\top}\right)\right\|^{p}_{\mathcal{S}_{p}}&\leq(1+\epsilon)\left\|\mathbf{A}-\mathbf{A}_{1}\right\|^{p}_{\mathcal{S}_{p}}\\\
&\leq\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}-\left\|\mathbf{A}\right\|^{p}_{\textrm{op}}+\Theta(\epsilon^{2/3})\end{split}$
(6.5)

By definition of the Schatten-$p$ norm we have:

$\begin{split}\left\|\mathbf{A}\left(\mathbf{I}-vv^{\top}\right)\right\|^{p}_{\mathcal{S}_{p}}&=\mathrm{Tr}\left(\left(\left(\mathbf{I}-vv^{\top}\right)^{2}\mathbf{A}^{2}\left(\mathbf{I}-vv^{\top}\right)^{2}\right)^{p/2}\right)\\\
&\geq\mathrm{Tr}\left(\left(\mathbf{I}-vv^{\top}\right)^{p}\mathbf{A}^{p}\left(\mathbf{I}-vv^{\top}\right)^{p}\right)\\\
&=\mathrm{Tr}\left(\mathbf{A}^{p}-\mathbf{A}^{p}vv^{\top}\right)\\\
&=\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}-\mathrm{Tr}\left(\left(vv^{\top}\right)^{p/2}\left(\mathbf{A}^{2}\right)^{p/2}\left(vv^{\top}\right)^{p/2}\right)\\\
&\geq\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}-\mathrm{Tr}\left(\left(vv^{\top}\mathbf{A}^{2}vv^{\top}\right)^{p/2}\right)\\\
&=\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}-\left\|\mathbf{A}vv^{\top}\right\|^{p}_{\mathcal{S}_{p}}\\\
&=\left\|\mathbf{A}\right\|^{p}_{\mathcal{S}_{p}}-\left\|\mathbf{A}v\right\|^{p}_{2}\end{split}$
(6.6)

where the first and last inequality follows from the reverse Araki-Lieb-
Thirring inequality (Fact 4.10). Combining equations (6.5) and (6.6), we have
that

$\left\|\mathbf{A}\right\|^{p}_{\textrm{op}}\geq\left\|\mathbf{A}v\right\|^{p}_{2}\geq\left\|\mathbf{A}\right\|^{p}_{\textrm{op}}-\Theta(\epsilon^{2/3})$
(6.7)

Next, we observe that $\mathbf{A}v=\left(\mathbf{I}-1/5\mathbf{W}\right)v$ can
be computed with one additional matrix-vector product and

$\left\|\mathbf{A}\right\|_{\textrm{op}}^{p}=\left(1-\frac{1}{5}\lambda_{n}(\mathbf{W})\right)^{p}=1-\frac{p}{5}\lambda_{n}(\mathbf{W})+\mathcal{O}\left(\lambda^{2}_{n}(\mathbf{W})\right)$
(6.8)

Consider the estimator
$\hat{\lambda}(\mathbf{W})=\frac{5}{p}\left(1-\left\|\left(\mathbf{I}-\frac{1}{5}\mathbf{W}\right)v\right\|_{2}^{p}\right)$.
Combining equations (6.7) and (6.8), we can conclude

$\hat{\lambda}(\mathbf{W})=\lambda_{\min}(\mathbf{W})\pm\Theta(\epsilon^{2/3}).$

obtaining an additive error estimate to the minimum eigenvalue of $\mathbf{W}$
by computing an additional matrix-vector product. It follows that we satisfy
conditions (1) and (2) in Theorem 6.2 and thus any algorithm for computing a
rank-$1$ approximation to the matrix
$\mathbf{A}=\mathbf{I}-\frac{1}{5}\mathbf{W}$ in Schatten $p$ norm must make
at least $\frac{1}{\epsilon^{1/3}}$ queries to the aforementioned matrix,
completing the proof. The claim follows from Theorem 6.2. ∎

###  6.2  Lower Bound for $p>2$

We now consider the case when $p>2$. We note that the previous approach no
longer works since we cannot lower bound the cost of
$\|\left(\mathbf{I}-\mathbf{W}/5\right)\left(\mathbf{I}-vv^{\top}\right)\|_{\mathcal{S}_{p}}$,
as the Araki-Lieb-Thirring inequality reverses (see application in Equation
6.6). Therefore, we require a new approach, and appeal to a special case of
Conjecture 2.2 that is known to be true, i.e. the Aligned Norm Compression
inequality (see Fact 4.14). The main theorem we prove in this sub-section is
as follows:

######  Theorem 6.4 (Query Lower Bound for $p>2$).

Given $\varepsilon>0$, and $p\geq 2$ such that $p=\mathcal{O}\left(1\right)$,
there exists a distribution $\mathcal{D}$ over $n\times n$ matrices such that
for $\mathbf{A}\sim\mathcal{D}$, any randomized algorithm that with
probability at least $99/100$ outputs a unit vector $u$ such that
$\left\|\mathbf{A}-\mathbf{A}uu^{\top}\right\|^{p}_{\mathcal{S}_{p}}\leq(1+\varepsilon)\left\|\mathbf{A}-\mathbf{A}_{1}\right\|^{p}_{\mathcal{S}_{p}}$
must make $\Omega\left(1/\varepsilon^{1/3}\right)$ matrix-vector queries to
$\mathbf{A}$.

We first introduce a sequence of key lemmas required for our proof.

######  Corollary 6.5 (Special Case of Lemma 5.2).

Given $\gamma\in[0,1]$, a vector $v\in\mathbb{R}^{d}$ and an $n\times d$
matrix $\mathbf{A}$, let $t=\log(n/\gamma)/(c\sqrt{\gamma})$, for a fixed
universal constant $c$. Then, there exists an algorithm that computes $t$
matrix-vector products with $\mathbf{A}$ and outputs a unit vector $u$ such
that with probability at least $99/100$,

$\left\|\mathbf{A}\right\|_{\textrm{op}}^{2}-\left\|\mathbf{A}u\right\|_{2}^{2}\leq
O\left(\gamma\sigma_{2}^{2}\right).$

where $\sigma_{2}$ is the second largest singular value of $\mathbf{A}$.

Next, we prove a key lemma relating the norm of a matrix to norms of
orthogonal projections applied to the matrix. We note that this lemma is
straight forward and holds for arbitrary vectors unit $u,v$ if Conjecture 2.2
holds. However, we show that we can transform our matrix to have structure
such that we can apply Fact 4.14 instead.

######  Lemma 6.6 (Orthogonal Projectors to Block Matrices ).

Given an $n\times d$ matrix $\mathbf{A}$, $p>2$ and unit vectors
$u\in\mathbb{R}^{d},v\in\mathbb{R}^{n}$, such that
$\left(\mathbf{I}-vv^{\top}\right)\mathbf{A}uu^{\top}=0$. Then, we have

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}\leq\left\|\begin{pmatrix}\left\|vv^{\top}\mathbf{A}uu^{\top}\right\|_{\mathcal{S}_{p}}&\left\|vv^{\top}\mathbf{A}\left(\mathbf{I}-uu^{\top}\right)\right\|_{\mathcal{S}_{p}}\\\
\
0&\left\|\left(\mathbf{I}-vv^{\top}\right)\mathbf{A}\left(\mathbf{I}-uu^{\top}\right)\right\|_{\mathcal{S}_{p}}\end{pmatrix}\right\|_{\mathcal{S}_{p}}.$

###### Proof.

Let $\mathbf{I}-vv^{\top}=\mathbf{Y}\mathbf{Y}^{\top}$, where $\mathbf{Y}$ has
$n-1$ orthonormal columns. Further, since $v$ and $\mathbf{Y}$ span disjoint
subspaces, and the union of their span is $\mathbb{R}^{n}$, the matrix
$\left(v\mid\mathbf{Y}\right)$, obtained by concatenating their columns is
unitary. Then, let $\mathbf{R}=\left(v\mid\mathbf{Y}\right)^{\top}$ and
observe, $\mathbf{R}$ has orthonormal rows and columns (since $\mathbf{R}$ is
unitary). Next, let $\mathbf{I}-uu^{\top}=\mathbf{Z}\mathbf{Z}^{\top}$, where
$\mathbf{Z}$ is $d\times(d-1)$ and has orthonormal columns. Let
$\SS=\left(u\mid\mathbf{Z}\right)^{\top}$, and observe $\SS$ has orthonormal
rows and columns.

Let $\hat{\mathbf{A}}=\mathbf{R}\mathbf{A}\SS^{\top}$, which admits the
following block-matrix form:

$\begin{split}\hat{\mathbf{A}}=\begin{pmatrix}v^{\top}\\\
\mathbf{Y}^{\top}\end{pmatrix}\cdot\mathbf{A}\cdot\left(u\mid\mathbf{Z}\right)=\begin{pmatrix}v^{\top}\\\
\mathbf{Y}^{\top}\end{pmatrix}\left(\mathbf{A}u\mid\mathbf{A}\mathbf{Z}\right)=\begin{pmatrix}v^{\top}\mathbf{A}u&v^{\top}\mathbf{A}\mathbf{Z}\\\
\mathbf{Y}^{\top}\mathbf{A}u&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\end{pmatrix}\end{split}$

Since $\mathbf{R}$ and $\SS$ are unitary, it follows from unitary invariance
of the Schatten-$p$ norm that

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}=\left\|\hat{\mathbf{A}}\right\|_{\mathcal{S}_{p}}=\left\|\begin{pmatrix}v^{\top}\mathbf{A}u&v^{\top}\mathbf{A}\mathbf{Z}\\\
\mathbf{Y}^{\top}\mathbf{A}u&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\end{pmatrix}\right\|_{\mathcal{S}_{p}}=\left\|\begin{pmatrix}v^{\top}\mathbf{A}u&v^{\top}\mathbf{A}\mathbf{Z}\\\
0&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\end{pmatrix}\right\|_{\mathcal{S}_{p}},$
(6.9)

where the last equality follows from observing that
$\left\|\mathbf{Y}^{\top}\mathbf{A}u\right\|_{F}=\left\|\mathbf{Y}\mathbf{Y}^{\top}\mathbf{A}uu^{\top}\right\|_{F}=\left\|\left(\mathbf{I}-vv^{\top}\right)\mathbf{A}uu^{\top}\right\|_{F}=0$
and therefore $\mathbf{Y}^{\top}\mathbf{A}u$ is a matrix of all $0$s. Next, we
append a set of $d-2$ columns of $0$’s to make the top left and top right
block the same size. Since this does not change the singular values, we have

$\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}=\left\|\begin{pmatrix}v^{\top}\mathbf{A}u&0&v^{\top}\mathbf{A}\mathbf{Z}\\\
0&0&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\end{pmatrix}\right\|_{\mathcal{S}_{p}}$
(6.10)

Next, we construct a rotation matrix $\mathbf{R}$ such that on right
multiplying a row vector by $\mathbf{R}$, the first $d-1$ coordinates remain
the same and on the remaining coordinates, the vector
$v^{\top}\mathbf{A}\mathbf{Z}$ gets mapped to $ce_{1}^{\top}$ for some scalar
$c$. Let $\SS$ be the $d-1\times d-1$ rotation matrix such that
$v^{\top}\mathbf{A}\mathbf{Z}\SS=ce_{1}^{\top}$. Then,
$\mathbf{R}=\begin{pmatrix}\mathbf{I}&0\\\ 0&\SS\end{pmatrix}$ and it is easy
to verify that $\mathbf{R}$ is unitary. Therefore,

$\begin{pmatrix}v^{\top}\mathbf{A}u&0&v^{\top}\mathbf{A}\mathbf{Z}\\\
0&0&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\end{pmatrix}\cdot\mathbf{R}=\begin{pmatrix}v^{\top}\mathbf{A}u&0&ce_{1}^{\top}\\\
0&0&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\SS\end{pmatrix}$

Now, we observe the final matrix above has a block matrix form we can apply
the Aligned Norm Compression inequality from Fact 4.14, with
$\alpha_{1}=v^{\top}\mathbf{A}u$, $\alpha_{2}=c$, $\beta_{1}=0$ and
$\beta_{2}=0$, and therefore

$\begin{split}\left\|\mathbf{A}\right\|_{\mathcal{S}_{p}}=\left\|\begin{pmatrix}v^{\top}\mathbf{A}u&0&ce_{1}^{\top}\\\
0&0&\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\SS\end{pmatrix}\right\|_{\mathcal{S}_{p}}&\leq\left\|\begin{pmatrix}\left\|v^{\top}\mathbf{A}u\right\|_{\mathcal{S}_{p}}&0&\left\|ce_{1}^{\top}\right\|_{\mathcal{S}_{p}}\\\
0&0&\left\|\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\SS\right\|_{\mathcal{S}_{p}}\end{pmatrix}\right\|_{\mathcal{S}_{p}}\\\
&=\left\|\begin{pmatrix}\left\|vv^{\top}\mathbf{A}uu^{\top}\right\|_{\mathcal{S}_{p}}&\left\|vv^{\top}\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{\mathcal{S}_{p}}\\\
0&\left\|\mathbf{Y}\mathbf{Y}^{\top}\mathbf{A}\mathbf{Z}\mathbf{Z}^{\top}\right\|_{\mathcal{S}_{p}}\end{pmatrix}\right\|_{\mathcal{S}_{p}}\end{split}$
(6.11)

where the last equality follows from unitary invariance and substituting the
definition of $\mathbf{Y}\mathbf{Y}^{\top}$ and $\mathbf{Z}\mathbf{Z}^{\top}$
completes the proof.

∎

######  Fact 6.7 (SVD of a $2\times 2$ Matrix).

Given a $2\times 2$ matrix $\mathbf{M}=\begin{pmatrix}a&b\\\ c&d\end{pmatrix}$
let $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$ be the SVD of $\mathbf{M}$.
Then,

$\mathbf{\Sigma}_{1,1}=\sqrt{\frac{a^{2}+b^{2}+c^{2}+d^{2}+\sqrt{\left(a^{2}+b^{2}-c^{2}-d^{2}\right)^{2}+4\left(ac+bd\right)^{2}}}{2}},$

and

$\mathbf{\Sigma}_{2,2}=\sqrt{\frac{a^{2}+b^{2}+c^{2}+d^{2}-\sqrt{\left(a^{2}+b^{2}-c^{2}-d^{2}\right)^{2}+4\left(ac+bd\right)^{2}}}{2}}.$

Now, we are ready to prove Theorem 6.4.

###### Proof of Theorem 6.4.

Let $\mathbf{A}=\mathbf{I}-\frac{1}{5}\mathbf{W}$ where $\mathbf{W}$ is an
$n\times n$ Wishart matrix as in the proof of Theorem 6.1 and we have by
hypothesis that there is an algorithm that with probability at least $99/100$,
outputs a unit vector $u$ such that
$\left\|\mathbf{A}\left(\mathbf{I}-uu^{\top}\right)\right\|_{\mathcal{S}_{p}}^{p}\leq(1+\varepsilon)\left\|\mathbf{A}-\mathbf{A}_{1}\right\|_{\mathcal{S}_{p}}^{p}$.