Sublinear Time Low-Rank Approximation of Positive Semidefinite Matrices

TL;DR

This paper presents a novel sublinear time algorithm for computing relative-error low-rank approximations of positive semidefinite matrices, significantly improving efficiency over previous methods and nearly matching theoretical lower bounds.

Contribution

The authors introduce a sublinear time algorithm for PSD matrix low-rank approximation that surpasses prior approaches and approaches optimality based on proven lower bounds.

Findings

Achieves relative-error low-rank approximation in sublinear time.

Outperforms previous $nnz(A)$ time algorithms based on oblivious subspace embeddings.

Extends techniques to spectral norm approximation and ridge regression for PSD matrices.

Abstract

We show how to compute a relative-error low-rank approximation to any positive semidefinite (PSD) matrix in sublinear time, i.e., for any $n \times n$ PSD matrix $A$, in $\tilde O(n \cdot poly(k/\epsilon))$ time we output a rank-$k$ matrix $B$, in factored form, for which $\|A-B\|_F^2 \leq (1+\epsilon)\|A-A_k\|_F^2$, where $A_k$ is the best rank-$k$ approximation to $A$. When $k$ and $1/\epsilon$ are not too large compared to the sparsity of $A$, our algorithm does not need to read all entries of the matrix. Hence, we significantly improve upon previous $nnz(A)$ time algorithms based on oblivious subspace embeddings, and bypass an $nnz(A)$ time lower bound for general matrices (where $nnz(A)$ denotes the number of non-zero entries in the matrix). We prove time lower bounds for low-rank approximation of PSD matrices, showing that our algorithm is close to optimal. Finally, we extend our techniques to give sublinear time algorithms for low-rank approximation of $A$ in the (often stronger) spectral norm metric $\|A-B\|_2^2$ and for ridge regression on PSD matrices.

Tables

Table 1. Table 1: Comparison of our results to prior work on low-rank approximation of PSD matrices. ‖ M ‖ ∗ = ∑ i = 1 n σ i ​ ( M ) subscript norm 𝑀 superscript subscript 𝑖 1 𝑛 subscript 𝜎 𝑖 𝑀 \|M\|_{*}=\sum_{i=1}^{n}\sigma_{i}(M) denotes the nuclear norm of matrix M 𝑀 M . The cited results all output B 𝐵 B (in factored form), which is itself PSD. In Theorem 12 we show how to modify our algorithm to satisfy this condition, and run in n ​ k ω ⋅ missing p ​ o ​ l ​ y ​ ( log ⁡ n / ϵ ) ⋅ 𝑛 superscript 𝑘 𝜔 missing 𝑝 𝑜 𝑙 𝑦 𝑛 italic-ϵ nk^{\omega}\cdot\mathop{\mathrm{missing}}{poly}(\log n/\epsilon) time. The table shows results that do not require incoherence assumptions on A 𝐴 A . For general PSD matrices, all known incoherence based results (see e.g., [ Git11 , GM13 ] ) degrade to Θ ​ ( n ω ) Θ superscript 𝑛 𝜔 \Theta(n^{\omega}) runtime. Additionally, as discussed, any general low-rank approximation algorithm can be applied to PSD matrices, with state-of-the-art approaches running in input-sparsity time [ CW13 , MM13 , NN13 ] . [ CW17 ] extends these results to the case where the output B 𝐵 B is restricted to be PSD. [ TYUC16 ] does the same but outputs B 𝐵 B with rank 2 ​ k 2 𝑘 2k . For a more in depth discussion of bounds obtained in prior work, see Section 1.4 .

Source	Runtime	Approximation Bound
[DM05]	$n{k}^{\omega -1}\cdot missingpoly\left(1/ϵ\right)$	${\Vert A-B\Vert }_F\le {\Vert A-A_k\Vert }_F+ϵ{\Vert A\Vert }_{\ast }$222Note that this bound is stated incorrectly as ${\Vert A-B\Vert }_F\le {\Vert A-A_k\Vert }_F+ϵ{\sum }_{i=1}^n{\left(A_{ii}\right)}^2$ in [DM05].
[KMT09]	$n{k}^{\omega -1}\cdot missingpoly\left(1/ϵ\right)$	${\Vert A-B\Vert }_F\le {\Vert A-A_k\Vert }_F+ϵn\cdot {\max }_i{A}_{ii}$
[GM13]	$\tilde{O}\left(n^2\right)+n{k}^{\omega -1}missingpoly\left(\log n/ϵ\right)$	${\Vert A-B\Vert }_F\le {\Vert A-A_k\Vert }_F+ϵ{\Vert A-A_k\Vert }_{\ast }$
[AGR16] + [BW09]	$n\log n\cdot missingpoly\left(k\right)$	${\Vert A-B\Vert }_F\le \left(k+1\right){\Vert A-A_k\Vert }_{\ast }$
[MM16]	$n{k}^{w-1}\cdot missingpoly\left(\log k/ϵ\right)$	${\Vert A^{1/2}-B\Vert }_F^2\le \left(1+ϵ\right){\Vert A^{1/2}-A_k^{1/2}\Vert }_F^2$
[CW17]	$O\left(\mathrm{nnz}\left(A\right)\right)+nmissingpoly\left(k/ϵ\right)$	${\Vert A-B\Vert }_F^2\le \left(1+ϵ\right){\Vert A-A_k\Vert }_F^2$
Our Results
Theorem 9	$n{k}^{\omega -1}\cdot missingpoly\left(\log n/ϵ\right)$	${\Vert A-B\Vert }_F^2\le \left(1+ϵ\right){\Vert A-A_k\Vert }_F^2$
Theorem 25	$n{k}^{\omega }\cdot missingpoly\left(\log n/ϵ\right)$	${\Vert A-B\Vert }_2^2\le \left(1+ϵ\right){\Vert A-A_k\Vert }_2^2+\frac{ϵ}{k}{\Vert A-A_k\Vert }_F^2$

Equations

$$\Vert A-B{\Vert }_F^2\le (1+ϵ)\Vert A-A_k{\Vert }_F^2,$$

$$(1-ϵ)\Vert A-PA{\Vert }_F^2\le \Vert C-PC{\Vert }_F^2\le (1+ϵ)\Vert A-PA{\Vert }_F^2.$$

$$(1-ϵ)\Vert C-CP{\Vert }_F^2\le \Vert R-RP{\Vert }_F^2+\Delta \le (1+ϵ)\Vert C-CP{\Vert }_F^2$$

$${\tau }_i^k(A)=a_i^T{\left(A{A}^T+\frac{\Vert A-A_k{\Vert }_F^2}{k}I\right)}^{+}{a}_i.$$

$$(1-ϵ)\Vert A-PA{\Vert }_F^2\le \Vert C-PC{\Vert }_F^2\le (1+ϵ)\Vert A-PA{\Vert }_F^2.$$

$${\tau }_i^k(A^{1/2})\le {\tilde{\tau }}_i^k(A^{1/2})\le 3{\tau }_i^k(A^{1/2}).$$

$$\frac{1}{2}\left(A+\lambda I\right)⪯\left(A^{1/2}S{S}^T{A}^{1/2}+\lambda I\right)⪯\frac{3}{2}\left(A+\lambda I\right)$$

$$x_i^T{\left(A^{1/2}S{S}^T{A}^{1/2}+\lambda I\right)}^{-1}{x}_i=\frac{1}{\lambda }{\left(A-AS(S^{T}AS+\lambda I{)}^{-1}{S}^{T}A\right)}_{i,i}.$$

$${\tau }_i^k(A)\le 2\sqrt{\frac{n}{k}}\cdot {\tau }_i^k(A^{1/2}).$$

$${\tau }_i^k(A^{1/2})$$

$${\tau }_i^k(A)$$

$$=x_i^T{A}^{1/2}{\left(A^2+\frac{\Vert A-A_k{\Vert }_F^2}{k}I\right)}^{-1}{A}^{1/2}{x}_i$$

$$=x_i^{T}U\hat{\Lambda }{U}^T{x}_i$$

$${\lambda }_i^2+\frac{1}{k}{\lambda }_i\sum _{j=k+1}^n{\lambda }_j\le 2\sqrt{\frac{n}{k}}\left({\lambda }_i^2+\frac{1}{k}\sum _{j=k+1}^n{\lambda }_j^2\right).$$

$${\lambda }_i^2+\frac{1}{k}{\lambda }_i\sum _{j=k+1}^n{\lambda }_j$$

$$\le {\lambda }_i^2+\frac{1}{\sqrt{n}\cdot k^{3/2}}\cdot n\sum _{j=k+1}^n{\lambda }_j^2$$

$$\le \sqrt{\frac{n}{k}}\left({\lambda }_i^2+\frac{1}{k}\sum _{j=k+1}^n{\lambda }_j^2\right)$$

$$A{S}_1{\left(S_1^T{A}^{T}A{S}_1+\frac{\Vert A{S}_1-(A{S}_1{)}_k{\Vert }_F^2}{k}I\right)}^{+}{S}_1^T{A}^T=U\bar{\Sigma }{U}^T$$

$$(1-ϵ)\Vert A{S}_1(I-P){\Vert }_F^2\le \Vert \tilde{A}(I-P){\Vert }_F^2+\Delta \le (1+ϵ)\Vert A{S}_1(I-P){\Vert }_F^2$$

$$(1-ϵ)\Vert A{S}_1(I-P){\Vert }_2^2-\frac{ϵ}{k}\Vert A-A_k{\Vert }_F^2\le \Vert \tilde{A}(I-P){\Vert }_2^2\le (1+ϵ)\Vert A{S}_1(I-P){\Vert }_2^2+\frac{ϵ}{k}\Vert A-A_k{\Vert }_F^2.$$

$$W=\underset{W|\mathrm{rank}(W)=k}{\mathrm{argmin}}\Vert S_4^{T}A{S}_1{S}_{3}W{Z}^T-S_4^{T}A{S}_1{\Vert }_F^2.$$

$$N=\underset{N\in {\mathbb{R}}^{n\times k}}{\mathrm{argmin}}\Vert S_5^{T}Q{N}^T-S_5^{T}A{\Vert }_F^2.$$

$$\Vert A-M{N}^T{\Vert }_F^2\le (1+ϵ)\Vert A-A_k{\Vert }_F^2.$$

$$\Vert A{S}_1-A{S}_{1}Z{Z}^T{\Vert }_2^2\le \frac{ϵ}{k}\Vert A-A_k{\Vert }_F^2.$$

$$\Vert \tilde{A}-{\tilde{A}}_{k_1}{\Vert }_F^2$$

$$\le \frac{3}{2}\Vert A{S}_1-(A{S}_1{)}_{k_1}{\Vert }_F^2-\Delta \le 603\Vert A-A_k{\Vert }_F^2$$

$$\frac{1}{2}\Vert A{S}_1(I-P){\Vert }_2^2-\frac{1}{2k_1}\Vert A-A_{k_1}{\Vert }_F^2\le \Vert \tilde{A}(I-P){\Vert }_2^2\le \frac{3}{2}\Vert A{S}_1(I-P){\Vert }_2^2+\frac{1}{2k_1}\Vert A-A_{k_1}{\Vert }_F^2.$$

$$\Vert \tilde{A}-{\tilde{A}}_{k_1}{\Vert }_2^2\le \Vert \tilde{A}-\tilde{A}{U}_{k_1}{U}_{k_1}^T{\Vert }_2^2\le \frac{3}{2}\Vert A{S}_1-(A{S}_1{)}_{k_1}{\Vert }_2^2+\frac{1}{2k_1}\Vert A-A_{k_1}\Vert \le \frac{4ϵ}{k}\Vert A-A_k{\Vert }_F^2$$

$$\Vert A{S}_1-A{S}_{1}Z{Z}^T{\Vert }_2^2$$

$$\le 4\Vert \tilde{A}-{\tilde{A}}_{k_1}{\Vert }_2^2+\frac{4}{k_1}\Vert \tilde{A}-{\tilde{A}}_{k_1}{\Vert }_F^2+\frac{1}{k_1}\Vert A-A_{k_1}{\Vert }_F^2$$

$$\le \frac{(4+4\cdot 603+1)ϵ}{k}\cdot \Vert A-A_k{\Vert }_F^2$$

$$\Vert A{S}_1-Q{Q}^{T}A{S}_1{\Vert }_F^2\le (1+ϵ)\Vert A-A_k{\Vert }_F^2.$$

$$M^{\ast }=\underset{M|\mathrm{rank}(M)=k}{\mathrm{argmin}}\Vert M{Z}^T-A{S}_1{\Vert }_F^2.$$

Full text

Sublinear Time Low-Rank Approximation of

Positive Semidefinite Matrices

Cameron Musco

MIT

cnmusco@mit.edu

  

David P. Woodruff

Carnegie Mellon University

dwoodruf@cs.cmu.edu

We show how to compute a relative-error low-rank approximation to any positive semidefinite (PSD) matrix in sublinear time, i.e., for any $n\times n$ PSD matrix $A$, in $\tilde{O}(n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon))$ time we output a rank-$k$ matrix $B$, in factored form, for which $\|A-B\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$, where $A_{k}$ is the best rank-$k$ approximation to $A$. When $k$ and $1/\epsilon$ are not too large compared to the sparsity of $A$, our algorithm does not need to read all entries of the matrix. Hence, we significantly improve upon previous $\operatorname{nnz}(A)$ time algorithms based on oblivious subspace embeddings, and bypass an $\operatorname{nnz}(A)$ time lower bound for general matrices (where $\operatorname{nnz}(A)$ denotes the number of non-zero entries in the matrix). We prove time lower bounds for low-rank approximation of PSD matrices, showing that our algorithm is close to optimal. Finally, we extend our techniques to give sublinear time algorithms for low-rank approximation of $A$ in the (often stronger) spectral norm metric $\|A-B\|_{2}^{2}$ and for ridge regression on PSD matrices.

1 Introduction

A fundamental task in numerical linear algebra is to compute a low-rank approximation of a matrix. Such an approximation can reveal underlying low-dimensional structure, can provide a compact way of storing a matrix in factored form, and can be quickly applied to a vector. Countless applications include clustering [DFK*+*04, FSS13, LBKW14, CEM*+*15], datamining [AFK*+*01], information retrieval [PRTV00], learning mixtures of distributions [AM05, KSV08], recommendation systems [DKR02], topic modeling [Hof03], and web search [AFKM01, Kle99].

One of the most well-studied versions of the problem is to compute a near optimal low-rank approximation with respect to the Frobenius norm. That is, given an $n\times n$ input matrix $A$ and an accuracy parameter $\epsilon>0$, output a rank-$k$ matrix $B$ for which:

[TABLE]

where for a matrix $C$, $\|C\|_{F}^{2}=\sum_{i,j}C_{i,j}^{2}$ is its squared Frobenius norm, and $A_{k}=\textrm{argmin}_{\textrm{rank-}k\textrm{ }B}\|A-B\|_{F}$. $A_{k}$ can be computed exactly using the singular value decomposition, but takes $O(n^{3})$ time in practice and $n^{\omega}$ time in theory, where $\omega\approx 2.373$ is the exponent of matrix multiplication.

In seminal work, Frieze, Kannan, and Vempala [FKV04] and Achlioptas and McSherry [AM07] show that using randomization and approximation, much faster runtimes are possible. Specifically, [FKV04] gives an algorithm that, assuming access to the row norms of $A$, outputs rank-$k$ $B$, in factored form, such that with good probability, $\|A-B\|_{F}^{2}\leq\|A-A_{k}\|_{F}^{2}+\epsilon\|A\|_{F}^{2}$. The algorithm runs in just $n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. However $\operatorname{nnz}(A)$ additional time is required to compute the row norms, where $\operatorname{nnz}(A)$ denotes the number of non-zero entries of $A$. Further, the guarantee achieved can be significantly weaker than (1), since the error is of the form $\epsilon\|A\|_{F}^{2}$ rather than $\epsilon\|A-A_{k}\|_{F}^{2}$. Note that $\|A-A_{k}\|_{F}^{2}\ll\|A\|_{F}^{2}$ precisely when $A$ is well-approximated by a rank-$k$ matrix. Related additive error algorithms with additional assumptions were given for tensors in [SWZ16].

Sarlós [Sar06] showed how to achieve (1) with constant probability in $\tilde{O}(\operatorname{nnz}(A)\cdot k/\epsilon)+n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. This was improved by Clarkson and Woodruff [CW13] who achieved $O(\operatorname{nnz}(A))+n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. See also work by Bourgain, Dirksen, and Nelson [BDN15], Cohen [Coh16], Meng and Mahoney [MM13], and Nelson and Nguyen [NN13] which further improved the degree in the $\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ term. For a survey, see [Woo14].

In the special case that $A$ is rank-$k$ and so $\|A-A_{k}\|_{F}^{2}=0$, (1) is equivalent to the well studied low-rank matrix completion problem [CR09]. Much attention has focused on completing incoherent low-rank matrices, whose singular directions are represented uniformly throughout the rows and columns and hence can be identified via uniform sampling and without fully accessing the matrix. Under incoherence (and often condition number) assumptions, a number of methods are able to complete a rank-$k$ matrix in $\tilde{O}(n\cdot\mathop{\mathrm{missing}}{poly}(k))$ time [JNS13, Har14].

For general matrices, without incoherence, it is not hard to see that $\Omega(\operatorname{nnz}(A))$ is a time lower bound: if one does not read a constant fraction of entries of $A$, with constant probability one can miss an entry much larger than all others, which needs to be included in the low-rank approximation.

1.1 Low-rank Approximation of Positive Semidefinite Matrices

An important class of matrices for which low-rank approximation is often applied is the set of positive semidefinite (PSD) matrices. These are real symmetric matrices with all non-negative eigenvalues. They arise for example as covariance matrices, graph Laplacians, Gram matrices (in particular, kernel matrices), and random dot product models [YS07]. In multidimensional scaling, low-rank approximation of PSD matrices in the Frobenius norm error metric (1) corresponds to the standard ‘strain minimization’ problem [CC00]. Completion of low-rank, or nearly low-rank (i.e., when $\|A-A_{k}\|_{F}^{2}\approx 0$), PSD matrices from few entries is important in applications such as quantum state tomography [GLF*+*10] and global positioning using local distances [SY05, YH38].

Due to its importance, a vast literature studies low-rank approximation of PSD matrices [DM05, ZTK08, KMT09, BW09, LKL10, GM13, WZ13, DLWZ14, WLZ16, TYUC16, LJS16, MM16, CW17]. However, known algorithms either run in at least $\operatorname{nnz}(A)$ time, do not achieve the relative-error guarantee of (1), or require strong incoherence assumptions111 Many of these algorithms satisfy the additional constraint that the low-rank approximation $B$ is PSD. This is also now known to be possible in $O(\operatorname{nnz}(A))$ time using sketching-based algorithms for general matrices [CW17]. (see Table 1).

At the same time, the simple $\Omega(\operatorname{nnz}(A))$ time lower bound for general matrices does not hold in the PSD case. Positive semidefiniteness ensures that for all $i,j$, $|A_{i,j}|\leq\max(A_{i,i},A_{j,j})$. So ‘hiding’ a large entry in $A$ requires creating a corresponding large diagonal entry. By reading the $n$ diagonal elements, an algorithm can avoid being tricked by this approach. While far from an algorithm, this argument raises the possibility that improved runtimes could be possible for PSD matrices.

1.2 Our Results

We give the first sublinear time relative-error low-rank approximation algorithm for PSD matrices. Our algorithm reads just $nk\cdot\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)$ entries of $A$ and runs in $nk^{\omega-1}\cdot\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)$ time (Theorem 9). With probability $99/100$ it outputs a matrix $B$ in factored form which satisfies (1). We critically exploit the intuition that large entries cannot ‘hide’ in PSD matrices, but surprisingly require no additional assumptions on $A$, such as incoherence or bounded condition number.

We complement our algorithm with an $\Omega(nk/\epsilon)$ time lower bound. The lower bound is information-theoretic, showing that any algorithm which reads fewer than this number of entries in the input PSD matrix cannot achieve the relative-error guarantee of (1) with constant probability. As our algorithm only reads $nk\cdot\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)$ entries of $A$, this is nearly optimal for constant $\epsilon$. We note that the actual time complexity of our algorithm is slower by a factor of $k^{\omega-2}$.

Finally, we show that our techniques can be extended to compute $B$ satisfying the spectral norm guarantee: $\|A-B\|_{2}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{2}^{2}+\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}$ using just $nk^{2}\cdot\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)$ accesses to $A$ and $nk^{\omega}\cdot\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)$ time (Theorem 25). This guarantee is often stronger than (1) when $\|A-A_{k}\|_{F}^{2}$ is large, and is important in many applications. For example, we use this result to solve the ridge regression problem $\min_{x\in\mathbb{R}^{n}}\|Ax-y\|_{2}^{2}+\lambda\|x\|_{2}^{2}$ up to $(1+\epsilon)$ relative error in $\tilde{O}\left(\frac{ns_{\lambda}^{\omega}}{\epsilon^{2\omega}}\right)$ time, where $s_{\lambda}=\operatorname{tr}((A^{2}+\lambda I)^{-1}A^{2})$ is the statistical dimension of the problem (see Theorem 27). Typically $s_{\lambda}\ll n$, so our runtime is sublinear and improves significantly on existing input-sparsity time results [ACW16]. For a summary of our results and comparison to prior work see, Table 1.

1.3 Algorithm Overview

The starting point for our approach is the fundamental fact that any matrix $A$ contains a subset of $O(k/\epsilon)$ columns, call them $C$, that span a relative-error rank-$k$ approximation to $A$ [DRVW06, DV06, DMM06]. Computing the best low-rank approximation to $A$ using an SVD requires access to all $\Theta(n^{2})$ dot products between the columns of the matrix. However, given $C$, just $n\cdot O(k/\epsilon)$ dot products are needed – to project the remaining columns of the matrix to the span of the subset.

Additionally, a subset of size $\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ can be identified using an intuitive approach known as adaptive sampling [DV06]: columns are iteratively added to the subset, with each new column being sampled with probability proportional to its norm outside the column span of the current subset. Formally, column $a_{i}$ is selected with probability $\frac{\|a_{i}-P_{C}a_{i}\|_{2}^{2}}{\|A-P_{C}A\|_{F}^{2}}$ where $P_{C}$ is the projection onto the current subset $C$. Computing these sampling probabilities requires knowing the norm of each $a_{i}$ along with its dot product with each column currently in $C$. So, overall this approach gives a relative-error low-rank approximation using just $n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ dot products between columns of $A$.

The above observation is surprising – not only does every matrix contain a small column subset witnessing a near optimal low-rank approximation, but also, such a witness can be found using significantly less information about the column span of the matrix than is required by a full SVD.

This fact is not immediately algorithmically useful, as computing the required dot products takes $\operatorname{nnz}(A)\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. However, given PSD $A$, we can write the eigendecomposition $A=U\Lambda U^{T}$ where $\Lambda$ is a non-negative diagonal matrix of eigenvalues, and let $A^{1/2}=U\Lambda^{1/2}U^{T}$ be the matrix square root of $A$. Since $A^{1/2}A^{1/2}=A$, the entry $A_{i,j}$ is just the dot product between the $i^{th}$ and $j^{th}$ columns of $A^{1/2}$. So with $A$ in hand, the dot products have been ‘precomputed’ and the above approach yields a low-rank approximation algorithm for $A^{1/2}$ running in just $n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. Note that, aligning with our initial intuition that reading the diagonal entries of $A$ is necessary to avoid the $\operatorname{nnz}(A)$ time lower bound for general matrices, the diagonal entries of $A$ are the column norms of $A^{1/2}$, and hence their values are critical to computing the adaptive sampling probabilities.

By the above argument, given PSD $A$, we can compute in $n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time a rank-$k$ orthogonal projection matrix $P\in\mathbb{R}^{n\times n}$ (in factored form) for which $\|A^{1/2}-A^{1/2}P\|_{F}^{2}\leq(1+\epsilon)\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}$. This approach can be implemented using adaptive sampling [DV06], sublinear time volume sampling [AGR16], or as shown in [MM16], recursive ridge leverage score sampling. The ridge leverage scores are a natural interpolation between adaptive sampling and the widely studied leverage scores, which, as we will see, have a number of additional algorithmically useful properties. As discussed in [MM16], the guarantee for $A^{1/2}$ is useful for a number of kernel learning methods such as kernel ridge regression. However, it is very different from our final goal. In fact, one can show that projecting to $P$ can yield an arbitrarily bad low-rank approximation to $A$ itself (see Appendix A).

We note that, since $P$ is constructed via column selection methods, it is possible to efficiently compute a factorization of $A^{1/2}PA^{1/2}$ (see Appendix A). Further, this matrix gives a near optimal low-rank approximation of $A$ if we use error parameter $\epsilon^{\prime}=\epsilon/\sqrt{n}$. This approach gives a first sublinear time algorithm, but it is significantly suboptimal. Namely, it requires reading $\tilde{O}(nk/\epsilon^{\prime})=\tilde{O}(n^{3/2}k/\epsilon)$ entries of $A$ and takes $n^{1.69}\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time using fast matrix multiplication.

To improve the dependence on $n$, we need a better understanding of how to perform ridge leverage score sampling on $A$ itself. We start by showing that the ridge leverage scores of $A^{1/2}$ are within a factor of $O(\sqrt{n/k})$ of the ridge leverage scores of $A$. By this bound, if we over-sample columns of $A$ by a factor of $O(\sqrt{n/k})$ using the ridge leverage scores of $A^{1/2}$ (computable via [MM16]), obtaining a sample of $\tilde{O}(\sqrt{n/k}\cdot k/\epsilon^{2})$ columns, the sample will be a so-called projection-cost preserving sketch (PCP) of $A$. The notion of a PCP was introduced in [CEM*+*15], where a matrix $C$ is defined to be an $(\epsilon,k)$-column PCP of $A$ if for all rank-$k$ projection matrices $P$:

[TABLE]

One important property of a PCP is that good low-rank approximations to $C$ translate to good low-rank approximations of $A$. More precisely, if $U$ is an $n\times k$ matrix with orthonormal columns for which $\|C-UU^{T}C\|_{F}^{2}\leq(1+\epsilon)\|C-C_{k}\|_{F}^{2}$, then $\|A-UU^{T}A\|_{F}^{2}\leq\frac{(1+\epsilon)^{2}}{(1-\epsilon)}\|A-A_{k}\|_{F}^{2}$.

Letting $C$ be the $n\times\tilde{O}(\sqrt{nk}/\epsilon^{2})$ submatrix which we sample via ridge leverage scores, we can apply an $\operatorname{nnz}(C)$ time algorithm to compute a subspace $U\in\mathbb{R}^{n\times k}$ whose columns span a near-optimal low-rank approximation of $C$, and hence of $A$ by the PCP property. Using standard sampling techniques, we can approximately project the columns of $A$ to $U$, producing our final solution. This gives time complexity $n^{3/2}\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$, improving slightly upon our first approach.

To reduce the time to linear in $n$, we must further reduce the size of $C$ by sampling a small subset of its rows, which themselves form a PCP. To find these rows, we cannot afford to use oblivious sketching techniques, which would take at least $\operatorname{nnz}(C)$ time, nor can we use our previous method for providing $O(\sqrt{n/k})$ overestimates to the ridge leverage scores, since $C$ is no longer PSD. In fact, the row ridge leverage scores of $C$ can be arbitrarily large compared to those of $A^{1/2}$.

The key idea to getting around this issue is that, since $C$ is a column PCP of $A$, projecting its columns onto $A$’s top eigenvectors gives a near optimal low-rank approximation. Further, we can show that the ridge leverage scores of $A^{1/2}$ (appropriately scaled) upper bound the standard leverage scores of this low-rank approximation. Sampling by these leverage scores is not enough to give a guarantee like (2) – they ignore the entire component of $C$ not falling in the span of $A$’s top eigenvectors and so may significantly distort projection costs over the matrix. Further, at this point, we have no idea how to estimate the row norms of $C$, or even its Frobenius norm, with $n\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ samples, which are necessary to implement any kind of adaptive sampling approach.

Fortunately, using that sampling at least preserves the matrix in expectation, along with a few other properties of the ridge leverage scores of $A^{1/2}$, we show that, with good probability, sampling $\tilde{O}(\sqrt{nk}/\mathop{\mathrm{missing}}{poly}(\epsilon))$ rows of $C$ by these scores yields $R$ satisfying for all rank-$k$ projection matrices $P$:

[TABLE]

where $\Delta$ is a fixed value, independent of $P$, with $|\Delta|\leq c\|C-C_{k}\|_{F}^{2}$ for some constant $c$. Since the same $\Delta$ distortion applies to all $P$, and since it is at most a constant times the true optimum, a near optimal low-rank approximation for $R$ still translates to a near optimal approximation for $C$.

At this point $R$ is a small matrix, and we can run any $O(\operatorname{nnz}(R))$ time algorithm to find a good low-rank factorization $EF^{T}$ to it, where $F^{T}$ is $k\times\tilde{O}(\sqrt{nk}/\mathop{\mathrm{missing}}{poly}(\epsilon))$. Since $R$ is a row PCP for $C$, by regressing the rows of $C$ to the span of $F$, we can obtain a near optimal low-rank approximation to $C$. We can solve this multi-response regression approximately in sublinear time via standard sampling techniques. Approximately regressing $A$ to the span of this approximation using similar techniques yields our final result. The total runtime is dominated by the input-sparsity low-rank approximation of $R$ requiring $O(\operatorname{nnz}(R))=\tilde{O}(nk/\mathop{\mathrm{missing}}{poly}(\epsilon))$ time.

To improve $\epsilon$ dependencies in our final runtime, achieving sample complexity $\tilde{O}\left(\frac{nk}{\epsilon^{2.5}}\right)$, we modify this approach somewhat, showing that $R$ actually satisfies a stronger spectral norm PCP property for $C$. This property lets us find a low-rank span $Z$ with $\|C-CZZ^{T}\|_{2}^{2}\leq\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}$, from which, through a series of approximate regression steps, we can extract a low-rank approximation to $A$ satisfying (1). This stronger spectral guarantee also lies at the core of our extensions to near optimal spectral norm low-rank approximation (Theorem 25), ridge regression (Theorem 27), and low-rank approximation where $B$ is restricted to be PSD (Theorem 12).

1.4 Some Further Intuition on Error Guarantees

Observe that in computing a low-rank approximation of $A$, we read just $\tilde{O}(n\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon))$ entries of the matrix, which is, up to lower order terms, the same number of entries (corresponding to column dot products of $A^{1/2}$) that we accessed to compute a low-rank approximation of $A^{1/2}$ in our description above. However, these sets of entries are very different. While low-rank approximation of $A^{1/2}$ looks at an $n\times\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ sized submatrix of $A$ together with the diagonal entries, our algorithm considers a carefully chosen $\sqrt{nk}\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)\times\sqrt{nk}\mathop{\mathrm{missing}}{poly}(\log n/\epsilon)$ submatrix together with the diagonal entries, which gives significantly more information about the spectrum of $A$.

As a simple example, consider $A$ with top eigenvalue $\lambda_{1}=\sqrt{n}$, and $\lambda_{i}=1$ for $i=2,...n$. $\|A^{1/2}\|_{F}^{2}=\sum_{i=1}^{n}\lambda_{i}=\sqrt{n}+n-1$ while $\|A^{1/2}-A^{1/2}_{1}\|_{F}^{2}=\sum_{i=2}^{n}\lambda_{i}=n-1$. So, $A^{1/2}$ has no good rank-$1$ approximation. Unless we set $\epsilon=O(1/\sqrt{n})$, a low-rank approximation algorithm for $A^{1/2}$ can learn nothing about $\lambda_{1}$ and still be near optimal. In contrast, $\|A\|_{F}^{2}=\sum_{i=1}^{n}\lambda_{i}^{2}=2n-1$ and $\|A-A_{1}\|_{F}^{2}=\sum_{i=2}^{n}\lambda_{i}^{2}=n-1$. So, even with $\epsilon=1/2$, any rank-$1$ approximation algorithm for $A$ must identify the presence of $\lambda_{1}$ and project this direction off the matrix. In this sense, our algorithm is able to obtain a much more accurate picture of $A$’s spectrum.

With incoherence assumptions, prior work on PSD low-rank approximation [GM13] obtains the bound $\|A-B\|_{*}\leq(1+\epsilon)\|A-A_{k}\|_{*}$ in sublinear time, where $\|M\|_{*}=\sum_{i=1}^{n}\sigma_{i}(M)$ is the nuclear norm of $M$. Recent work ([AGR16] in combination with [BW09]) gives $\|A-B\|_{F}\leq(k+1)\|A-A_{k}\|_{*}$ without the incoherence assumption. These nuclear norm bounds are closely related to approximation bounds for $A^{1/2}$ and it is not hard to see that neither require $\lambda_{1}$ to be detected in the example above, and so in this sense are weaker than our Frobenius norm bound.

A natural question if even stronger bounds are possible: e.g., can we compute $B$ with $\|A-B\|_{2}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{2}^{2}$ in sublinear time? We partially answer this question in Theorem 25. In $\tilde{O}(nk^{\omega}\mathop{\mathrm{missing}}{poly}(\log n/\epsilon))$ time, we can find $B$ satisfying $\|A-B\|_{2}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{2}^{2}+\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}$.

Significantly improving the above bound seems difficult: it is easy to see that a relative error spectral norm guarantee requires $\Omega(n^{2})$ time. Consider $A$ which is the identity except with $A_{i,j}=A_{j,i}=1$ for some uniform random pair $(i,j)$. Finding $(i,j)$ requires $\Omega(n^{2})$ queries to $A$. However, it is necessary to achieve a relative error spectral norm guarantee with $\epsilon<3$ since $\|A\|_{2}^{2}=4$ while $\|A-A_{1}\|_{2}^{2}=1$ where $A_{1}$ is all zeros with ones at its $(i,i)$, $(j,j)$, $(i,j)$, and $(j,i)$ entries.

A similar argument shows that relative error low-rank approximation in higher Schatten-$p$ norms, i.e., $\|A-B\|_{p}^{p}$ for $p>2$ requires superlinear dependence on $n$ (where $\|M\|_{p}^{p}=\sum_{i=1}^{n}\sigma_{i}^{p}(M)$.) We can set $A$ to be the identity but with an all ones block on a uniform random subset of $n^{1/p}$ indices. This block has associated eigenvalue $\lambda_{1}=n^{1/p}$ and so, since all other ($n-n^{1/p}$) eigenvalues of $A$ are $1$, $\|A\|_{p}^{p}=\Theta(n)$, and the block must be recovered to give a relative error approximation to $\|A-A_{1}\|_{p}^{p}$. However, as the block is placed uniformly at random and contains just $n^{2/p}$ entries, finding even a single entry requires $n^{2-2/p}$ queries to $A$ – superlinear for $p>2$.

1.5 Open Questions

While it is apparent that obtaining stronger error guarantees than (1) may require increased runtime, understanding exactly what can be achieved in sublinear time is an interesting direction for future work. We also note that it is still unknown how to compute a number of basic properties of PSD matrices in sublinear time. For example, while we can output $B$ satisfying $\|A-B\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$, surprisingly it is not clear how to actually estimate the value $\|A-A_{k}\|_{F}^{2}$ to within a $(1\pm\epsilon)$ factor. This can be achieved in $n^{3/2}\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time using our PCP techniques. However, obtaining linear runtime in $n$ is open. Estimating $\|A-A_{k}\|_{F}^{2}$ seems strongly connected to estimating other important quantities such as the statistical dimension of $A$ for ridge regression (see Theorem 27) which we do not know how to do in $o(n^{3/2})$ time.

Finally, an open question is if these techniques can be generalized to a broader class of matrices. As discussed, in the matrix completion literature, much attention has focused on incoherent low-rank matrices [CR09] which can be approximated with uniform sampling. PSD matrices are not incoherent in general, which is highlighted by the fact that our sampling schemes are far from uniform and very adaptive to previously seen matrix entries. However, perhaps there is some other parameter (maybe relating to a measure of diagonal dominance) which characterizes when low-rank approximation can be performed with just a small number of adaptive accesses to $A$.

1.6 Paper Outline

Section 2: Ridge Leverage Score Sampling. We show that the ridge leverage scores of $A$ are within an $O(\sqrt{n/k})$ factor of those of $A^{1/2}$, letting us use the fast ridge leverage score sampling algorithm of [MM16] to sample $\tilde{O}(\sqrt{nk}/\epsilon^{2})$ columns of $A$ that form a column PCP of the matrix.

Section 3: Row Sampling. We discuss how to further accelerate our algorithm by obtaining a row PCP for our column sample, allowing us to achieve runtime linear in $n$.

Section 4: Full Algorithm. We use the primitives in the previous sections along with approximate regression techniques to give our full sublinear time low-rank approximation algorithm.

Section 5: Lower Bounds. We show that our algorithm is nearly optimal – any relative error low-rank approximation algorithm must read $\Omega(nk/\epsilon)$ entries of $A$.

Section 6: Spectral Norm Bounds. We modify the algorithm of Section 4 to give a tighter approximation in the spectral norm and discuss applications to sublinear time ridge regression.

2 Ridge Leverage Score Sampling

Our main algorithmic tool will be ridge leverage score sampling, which is used to identify a small subset of columns of $A$ that span a good low-rank approximation of the matrix. Following the definition of [CMM17], the rank-$k$ ridge leverage scores of any matrix $A$ are given by:

Definition 1** (Ridge Leverage Scores).**

For any $A\in\mathbb{R}^{n\times d}$, letting $a_{i}\in\mathbb{R}^{n}$ be the $i^{th}$ column of $A$, the $i^{th}$ rank-$k$ column ridge leverage score of $A$ is:

[TABLE]

Above $I$ is the appropriately sized identity matrix and $M^{+}$ denotes the matrix pseudoinverse, equivalent to the inverse unless $\|A-A_{k}\|_{F}^{2}=0$ and $A$ is singular. Analogous scores can be defined for the rows of $A$ by simply transposing the matrix. It is not hard to see that $0<\tau^{k}_{i}(A)<1$ for all $i$. Since we use these scores as sampling probabilities, it is critical that the sum of scores, and hence the size of the subsets we sample, is not too large. We have the following (see Appendix B):

Lemma 2** (Sum of Ridge Leverage Scores).**

For any $A\in\mathbb{R}^{n\times d}$, $\sum_{i=1}^{d}\tau^{k}_{i}(A)\leq 2k$.

Intuitively, the ridge leverage scores are similar to the standard leverage scores of $A$, which are given by $a_{i}^{T}(AA^{T})^{+}a_{i}$. By writing $A=U\Sigma V^{T}$ in its SVD, one sees that standard leverage scores are just the squared column norms of $V^{T}$. Sampling columns by ridge leverage scores yields a spectral approximation to the matrix. The addition of the weighted identity (or ‘ridge’) $\frac{\|A-A_{k}\|_{F}^{2}}{k}I$ ‘dampens’ contributions from smaller singular directions of $A$, decreasing the sum of the scores and allowing us to sample fewer columns. At the same time, it introduces error dependent on the size of the tail $\|A-A_{k}\|_{F}^{2}$, ultimately giving an approximation from which it is possible to output a near optimal low-rank approximation to the original matrix. Specifically, sampling by ridge leverage scores yields a projection-cost preserving sketch (PCP) of $A$:

Lemma 3** (Theorem 6 of [CMM17]).**

For any $A\in\mathbb{R}^{n\times d}$, for $i\in\{1,\ldots,d\}$, let $\tilde{\tau}_{i}^{k}\geq\tau_{i}^{k}(A)$ be an overestimate for the $i^{th}$ rank-$k$ ridge leverage score. Let $p_{i}=\frac{\tilde{\tau}^{k}_{i}}{\sum_{i}\tilde{\tau}^{k}_{i}}$ and $t=\frac{c\log(k/\delta)}{\epsilon^{2}}\sum_{i}\tilde{\tau}^{k}_{i}$ for any $\epsilon<1$ and sufficiently large constant $c$. Construct $C$ by sampling $t$ columns of $A$, each set to $\frac{1}{\sqrt{tp_{i}}}a_{i}$ with probability $p_{i}$. With probability $1-\delta$, for any rank-$k$ orthogonal projection $P\in\mathbb{R}^{n\times n}$,

[TABLE]

We refer to $C$ as an $(\epsilon,k)$-column PCP of $A$.

Since the ‘cost’ $\|A-PA\|_{F}^{2}$ of any rank-$k$ projection of $A$ is preserved by $C$, any near-optimal low-rank approximation of $C$ yields a near optimal low-rank approximation of $A$. Further, $C$ is much smaller than $A$, so such a low-rank approximation can be computed quickly. The difficulty is in computing the approximate leverage scores. To do this, we use the main result from [MM16]:

Lemma 4** (Corollary of Theorem 20 of [MM16]).**

There is an algorithm that given any PSD matrix $A\in\mathbb{R}^{n\times n}$, runs in $O(n(k\log(k/\delta))^{\omega-1})$ time, accesses $O(nk\log(k/\delta))$ entries of $A$, and returns for each $i\in[1,..,n]$, $\tilde{\tau}^{k}_{i}(A^{1/2})$ such that with probability $1-\delta$, for all $i$:

[TABLE]

Proof.

Theorem 20 of [MM16] shows that by using a recursive ridge leverage score sampling algorithm, it is possible to return (with probability $1-\delta$) a sampling matrix $S\in\mathbb{R}^{n\times s}$ with $s=O\left(k\log(k/\delta)\right)$ such that, letting $\lambda=\frac{1}{k}\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}$:

[TABLE]

where $M\preceq N$ indicates $x^{T}Mx\leq x^{T}Nx$ for all $x$. If we set $\tilde{\tau}^{k}_{i}(A^{1/2})=2\cdot x_{i}^{T}\left(A^{1/2}SS^{T}A^{1/2}+\lambda I\right)^{+}x_{i}$, where $x_{i}$ is the $i^{th}$ column of $A^{1/2}$ we have the desired bound. Of course, we cannot directly compute this value without factoring $A$ to form $A^{1/2}$. However, as shown in Lemma 6 of [MM16]:

[TABLE]

Computing $(S^{T}AS+\lambda I)^{-1}$ requires $O(s^{2})=O((k\log(k/\delta))^{2})$ accesses to $A$ and $O(s^{\omega})=O((k\log(k/\delta))^{\omega})$ time. Computing all $n$ diagonal entries of $AS(S^{T}AS+\lambda I)^{-1}S^{T}A$ then requires $O(nk\log(k/\delta))$ accesses to $A$ and $O(n(k\log(k/\delta))^{\omega-1})$ time. With these entries in hand we can simply subtract from the diagonal entries of $A$ and rescale to give the final leverage score approximation. Critically, this calculation always reads all diagonal entries of $A$, allowing it to identify rows containing large off diagonal entries and skirt the $\operatorname{nnz}(A)$ time lower bound for general matrices.

Note that the stated runtime in [MM16] for outputting $S$ is $\tilde{O}(nk)$ accesses to $A$ (kernel evaluations in the language of [MM16]) and $\tilde{O}(nk^{2})$ runtime. However this runtime is improved to $\tilde{O}(nk^{\omega-1})$ using fast matrix multiplication. ∎

In order to apply Lemmas 3 and 4 to low-rank approximation of $A$, we now show that the ridge leverage scores of $A^{1/2}$ coarsely approximate those of $A$:

Lemma 5** (Ridge Leverage Score Bound).**

For any PSD matrix $A\in\mathbb{R}^{n\times n}$:

[TABLE]

Proof.

We write $A^{1/2}$ in its eigendecomposition $A^{1/2}=U\Lambda^{1/2}U^{T}$, where $\Lambda_{i,i}=\lambda_{i}$ is the $i^{th}$ eigenvalue of $A$. Letting $x_{i}$ denote the $i^{th}$ column of $A^{1/2}$ we have:

[TABLE]

where $\bar{\Lambda}_{i,i}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\frac{1}{\lambda_{i}+\frac{1}{k}\sum_{j=k+1}^{n}\lambda_{j}}.$ We can similarly write:

[TABLE]

where $\hat{\Lambda}_{i,i}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\frac{\lambda_{i}}{\lambda_{i}^{2}+\frac{1}{k}\sum_{j=k+1}^{n}\lambda_{j}^{2}}.$ Showing $\hat{\Lambda}\preceq 2\sqrt{\frac{n}{k}}\cdot\bar{\Lambda}$ is enough to give the lemma. Specifically we must show, for all $i$, $\hat{\Lambda}_{i,i}\leq 2\sqrt{\frac{n}{k}}\cdot\bar{\Lambda}_{i,i}$ which after cross-multiplying is equivalent to:

[TABLE]

First consider the relatively large eigenvalues. Say we have $\frac{1}{k}\sum_{j=k+1}^{n}\lambda_{j}\leq\sqrt{\frac{n}{k}}\lambda_{i}$. Then:

[TABLE]

which gives (3). Next consider small eigenvalues with $\frac{1}{k}\sum_{j=k+1}^{n}\lambda_{j}\geq\sqrt{\frac{n}{k}}\lambda_{i}$. In this case:

[TABLE]

which gives (3), completing the proof. ∎

Combining Lemmas 2, 3, 4, 5 we have:

Corollary 6** (Fast PSD Ridge Leverage Score Sampling).**

There is an algorithm that given any PSD matrix $A\in\mathbb{R}^{n\times n}$ runs in $\tilde{O}(nk^{\omega-1})$ time, accesses $\tilde{O}(nk)$ entries of $A$, and with prob. $1-\delta$ outputs a weighted sampling matrix $S_{1}\in\mathbb{R}^{n\times\tilde{O}\left(\frac{\sqrt{nk}}{\epsilon^{2}}\right)}$ such that $AS_{1}$ is an $(\epsilon$,$k)$-column PCP of $A$.

Proof.

By Lemma 4 we can compute constant factor approximations to the ridge leverage scores of $A^{1/2}$ in time $\tilde{O}(nk^{\omega-1})$. Applying Lemma 5, if we scale these scores up by $2\sqrt{n/k}$ they will be overestimates of the ridge leverage scores of $A$. If we set $t=O\left(\frac{\log(k/\delta)}{\epsilon^{2}}\cdot\sum\tilde{\tau}_{i}^{k}\right)$, and generate $S_{1}$ by sampling $t$ columns of $A$ with probabilities proportional to these estimated scores, by Lemma 3, $AS_{1}$ will be an $(\epsilon,k)$-column PCP of $A$ with probability $1-\delta$. By Lemma 2, $\sum_{i=1}^{n}\tau_{i}^{k}(A^{1/2})\leq 2k$. So we have $t=\tilde{O}(\sum\tilde{\tau}_{i}^{k}/\epsilon^{2})=\tilde{O}(\sqrt{nk}/\epsilon^{2})$. ∎

Forming $AS_{1}$ requires reading just $\tilde{O}(n^{3/2}\sqrt{k}/\epsilon^{2})$ entries of $A$. At this point, we could employ any input sparsity time algorithm to find a near optimal rank-$k$ projection $P$ for approximating $AS_{1}$ in $O(\operatorname{nnz}(AS_{1}))+n\mathop{\mathrm{missing}}{poly}(k/\epsilon)=n^{3/2}\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. This would in turn yield a near optimal low-rank approximation of $A$. However, as we will see in the next section, by further sampling the rows of $AS_{1}$, we can significantly improve this runtime.

3 Row Sampling

To achieve near linear dependence on $n$, we sample roughly $\sqrt{nk}$ rows from $AS_{1}$, producing an even smaller matrix $S_{2}^{T}AS_{1}$, which we can afford to fully read and from which we can form a near optimal low-rank approximation to $AS_{1}$ and consequently to $A$. However, sampling $AS_{1}$ is challenging: we cannot employ input sparsity time methods as we cannot afford to read the full matrix, and since it is no longer PSD, we cannot apply the same approach we used for $A$, approximating the ridge leverage scores with those of $A^{1/2}$.

Rewriting Definition 1 using the SVD $AS_{1}=U\Sigma V^{T}$ (and transposing $AS_{1}$ to give row instead of column scores) we see that the row ridge leverage scores of $AS_{1}$ are the diagonal entries of:

[TABLE]

where $\bar{\Sigma}_{i,i}=\frac{\Sigma_{i,i}^{2}}{\Sigma_{i,i}^{2}+\frac{\|AS_{1}-(AS_{1})_{k}\|_{F}^{2}}{k}}$. That is, the row ridge leverage scores depend only on the column span $U$ of $AS_{1}$ and its spectrum. Since $AS_{1}$ is a column PCP of $A$ this gives hope that the two matrices have similar row ridge leverage scores.

Unfortunately, this is not the case. It is possible to have rows in $AS_{1}$ with ridge leverage scores significantly higher than in $A$. Thus, even if we knew the ridge leverage scores of $A$, we would have to scale them up significantly to sample from $AS_{1}$. As an example, consider $A$ with relatively uniform ridge leverage scores: $\tau_{i}(A)\approx k/n$ for all $i$. When a column is selected to be included in $AS_{1}$ it will be reweighted by roughly a factor of $\sqrt{n/k}$. Now, append a number of rows to $A$ each with very small norm and just a containing single non-zero entry. These rows will have little effect on the ridge leverage scores if their norms are small enough. However, if the column corresponding to the nonzero in a row is selected, the row will appear in $AS_{1}$ with $\sqrt{n/k}$ times the weight that it appears in $A$, and its ridge leverage score will be roughly a factor $n/k$ times higher.

Fortunately, we are still able to show that sampling the rows of $AS_{1}$ by the rank $k^{\prime}=O(k/\epsilon)$ leverage scores of $A^{1/2}$ scaled up by a $\sqrt{n/k^{\prime}}$ factor yields a row PCP for this matrix. Our proof works not with the ridge scores of $AS_{1}$ but with the standard leverage scores of a near optimal low-rank approximation to this matrix – specifically the approximation given by projecting onto the top eigenvectors of $A$. We have:

Lemma 7** (Row PCP).**

For any PSD $A\in\mathbb{R}^{n\times n}$ and $\epsilon\leq 1$ let $k^{\prime}=\lceil ck/\epsilon\rceil$ and let $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})\geq\tau_{i}^{k^{\prime}}(A^{1/2})$ be an overestimate for the $i^{th}$ rank-$k^{\prime}$ ridge leverage score of $A^{1/2}$. Let $\tilde{\ell}_{i}=\sqrt{\frac{16n\epsilon}{k}}\cdot\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})$, $p_{i}=\frac{\tilde{\ell}_{i}}{\sum_{i}\tilde{\ell}_{i}}$, and $t=\frac{c^{\prime}\log n}{\epsilon^{2}}\sum_{i}\tilde{\ell}_{i}$. Construct weighted sampling matrices $S_{1},S_{2}\in\mathbb{R}^{n\times t}$ each whose $j^{th}$ column is set to $\frac{1}{\sqrt{tp_{i}}}e_{i}$ with probability $p_{i}$.

For sufficiently large constants $c,c^{\prime}$, with probability $\frac{99}{100}$, letting $\tilde{A}=S_{2}^{T}AS_{1}$, for any rank-$k$ orthogonal projection $P\in\mathbb{R}^{t\times t}$:

[TABLE]

for some fixed $\Delta$ (independent of $P$) with $|\Delta|\leq 600\|A-A_{k}\|_{F}^{2}$. We refer to $\tilde{A}$ as an $(\epsilon,k)$-row PCP of $AS_{1}$.

Note that by Lemma 2, $\sum_{i}\tau_{i}^{k^{\prime}}(A^{1/2})=O(k/\epsilon)$. So if $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})$ is a constant factor approximation to $\tau_{i}^{k^{\prime}}(A^{1/2})$, $t=O\left(\frac{\sqrt{nk}\log n}{\epsilon^{2.5}}\right)$. Also note that the Lemma requires both $S_{1}$ and $S_{2}$ to be sampled using the rank $k^{\prime}$ ridge scores. If we sample $S_{1}$ using a sum of the rank-$k$ and rank-$k^{\prime}$ ridge scores (appropriately scaled) Lemma 7 and Lemma 3 will hold simultaneously.

By applying an input sparsity time low-rank approximation algorithm to $\tilde{A}$ (which has just $\tilde{O}\left(\frac{nk}{\epsilon^{5}}\right)$ entries) we can find a near optimal low-rank approximation of $AS_{1}$, and thus for $A$. However, in our final algorithm, we will take a somewhat different approach. We are able to show that using appropriate sampling probabilities, we can in fact sample $\tilde{A}$ which is a projection-cost preserving sketch of $AS_{1}$ for spectral norm error. As we will see, recovering a near optimal spectral norm low-rank approximation to $AS_{1}$ suffices to recover a near optimal Frobenius norm approximation to $A$, and allows us to improve $\epsilon$ dependencies in our final runtime.

Lemma 8** (Spectral Norm Row PCP).**

For any PSD $A\in\mathbb{R}^{n\times n}$, and $\epsilon<1$ let $k^{\prime}=\lceil ck/\epsilon^{2}\rceil$ and $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})\geq\tau_{i}^{k^{\prime}}(A^{1/2})$ be an overestimate for the $i^{th}$ rank-$k^{\prime}$ ridge leverage score of $A^{1/2}$. Let $\tilde{\ell}_{i}=4\epsilon\sqrt{\frac{n}{k}}\tau_{i}^{k^{\prime}}(A^{1/2})$, $p_{i}=\frac{\tilde{\ell}_{i}}{\sum_{i}\tilde{\ell}_{i}}$, and $t=\frac{c^{\prime}\log n}{\epsilon^{2}}\cdot\sum_{i}\tilde{\ell}_{i}$. Construct weighted sampling matrices $S_{1},S_{2}\in\mathbb{R}^{n\times t}$, each whose $j^{th}$ column is set to $\frac{1}{\sqrt{tp_{i}}}e_{i}$ with probability $p_{i}$.

For sufficiently large constants $c,c^{\prime}$, with high probability (i.e. probability $\geq 1-1/n^{d}$ for some large constant $d$), letting $\tilde{A}=S_{2}^{T}AS_{1}$, for any orthogonal projection $P\in\mathbb{R}^{t\times t}$:

[TABLE]

We refer to $\tilde{A}$ as an $(\epsilon,k)$-spectral PCP of $AS_{1}$.

Note that if $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})$ is a constant factor approximation to $\tau_{i}^{k^{\prime}}(A^{1/2})$, $t=O\left(\frac{\sqrt{nk}\log n}{\epsilon^{3}}\right)$. We defer the proofs of Lemmas 7 and 8 to Appendix B.

4 Full Low-Rank Approximation Algorithm

We are finally ready to give our main algorithm for relative error low-rank approximation of PSD matrices in $\tilde{O}(n\mathop{\mathrm{missing}}{poly}(k/\epsilon))$ time, Algorithm 1. We set $k_{1}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\lceil ck/\epsilon\rceil$ and estimate the both the rank-$k$ and rank-$c^{\prime}k_{1}$ ridge leverage scores of $A^{1/2}$ using the algorithm of [MM16] (Step 1). If $c,c^{\prime}$ are sufficiently large, sampling by the sum of these scores (Steps 2-3) ensures that $AS_{1}$ is an $(\epsilon,k)$-column PCP for $A$ and simultaneously, by applying Lemmas 7 and 8 that $\tilde{A}=S_{2}^{T}AS_{1}$ is a row PCP in both spectral and Frobenius norm with rank $k_{1}$ and error $\epsilon=1/2$ for $AS_{1}$

In conjunction, these guarantees ensure that we can apply an input sparsity time algorithm to $\tilde{A}$ (Step 4) to find a rank-$k_{1}$ $Z$ satisfying $\|AS_{1}-AZZ^{T}\|_{2}^{2}=O(\|AS_{1}-(AS_{1})_{k_{1}}\|_{2}^{2}+\frac{1}{k_{1}}\|A-A_{k}\|_{F}^{2})=O\left(\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}\right)$, where the final bound holds since $k_{1}=\Theta(k/\epsilon)$. It is not hard to show that, due to this strong spectral norm bound, projecting $AS_{1}$ to $Z$ and taking the best rank-$k$ approximation in the span gives a near optimal Frobenius norm low-rank approximation to $AS_{1}$ and hence $A$.

We can still not afford to read $AS_{1}$ in its entirety, so we employ a number of standard leverage score sampling techniques to perform this projection approximately. In Step 5, we sample $\tilde{O}(k/\epsilon^{2})$ columns of $AS_{1}$ using the leverage scores of $Z$ (its row norms since it is an orthonormal matrix) to form $AS_{1}S_{3}$. We argue that there is a good rank-$k$ approximation to $AS_{1}$ lying in both the column span of $AS_{1}S_{3}$ and the row span of $Z^{T}$. In Step 6 we find a near optimal such approximation by further sampling $\tilde{O}(k/\epsilon^{4})$ rows $AS_{1}$ by the leverage scores of $AS_{1}S_{3}$ (the row norms of $V$, an orthonormal basis for its span), and computing the best rank-$k$ approximation to the sampled matrix falling in the column span of $AS_{1}S_{3}$ and the row span of $Z^{T}$.

Finally, in Step 7 we approximately project $A$ itself to the span of this rank-$k$ approximation by first sampling by the leverage scores of the approximation (the row norms of $Q$) and projecting.

Algorithm 1

PSD Low-Rank Approximation

Let $k_{1}=\lceil ck/\epsilon\rceil$. For all $i\in[1,..,n]$ compute $\tilde{\tau}_{i}^{k}(A^{1/2})$ and $\tilde{\tau}_{i}^{c^{\prime}k_{1}}(A^{1/2})$ which are constant factor approximations to the ridge leverage scores $\tau_{i}^{k}(A^{1/2})$ and $\tau_{i}^{c^{\prime}k_{1}}(A^{1/2})$ respectively. 2. 2.

Set $\ell^{(1)}_{i}=\sqrt{\frac{n}{k}}\tilde{\tau}_{i}^{k}(A^{1/2})+\sqrt{\frac{n\epsilon^{4}}{k_{1}}}\tilde{\tau}_{i}^{c^{\prime}k_{1}}(A^{1/2})$ and $\ell^{(2)}_{i}=\sqrt{\frac{n}{k_{1}}}\tilde{\tau}_{i}^{c^{\prime}k_{1}}(A^{1/2})$. Set $p_{i}^{(1)}=\frac{\ell_{i}^{(1)}}{\sum_{i}\ell_{i}^{(1)}}$ and $p_{i}^{(2)}=\frac{\ell_{i}^{(2)}}{\sum_{i}\ell_{i}^{(2)}}$. 3. 3.

Set $t_{1}=\frac{c_{1}\log n}{\epsilon^{2}}\sum_{i}\ell_{i}^{(1)}$ and $t_{2}=c_{2}\log n\sum_{i}\ell_{i}^{(2)}$. Sample $S_{1}\in\mathbb{R}^{n\times t_{1}}$ whose $j^{th}$ column is set to $\frac{1}{\sqrt{tp^{(1)}_{i}}}e_{i}$ with probability $p^{(1)}_{i}$. Sample $S_{2}\in\mathbb{R}^{n\times t_{2}}$ analogously using $p^{(2)}_{i}$. 4. 4.

Let $\tilde{A}=S_{2}^{T}AS_{1}$, and use an input sparsity time algorithm to compute orthonormal $Z\in\mathbb{R}^{t_{1}\times k_{1}}$ satisfying the spectral guarantee $\|\tilde{A}-\tilde{A}ZZ^{T}\|_{2}^{2}\leq 2\|\tilde{A}-\tilde{A}_{k_{1}}\|_{2}^{2}+\frac{2}{k_{1}}\|\tilde{A}-\tilde{A}_{k_{1}}\|_{F}^{2}$. 5. 5.

Let $t_{3}=c_{3}\left(\frac{k\log(k/\epsilon)}{\epsilon}+\frac{k}{\epsilon^{2}}\right)$, set $p_{i}^{(3)}=\frac{\|z_{i}\|_{2}^{2}}{\|Z\|_{F}^{2}}$, and sample $S_{3}\in\mathbb{R}^{t_{1}\times t_{3}}$ where the $j^{th}$ column is set to $\frac{1}{\sqrt{t_{3}p_{i}^{(3)}}}e_{i}$ with probability $p_{i}^{(3)}$. Compute $V\in\mathbb{R}^{n\times t_{3}}$ which is an orthogonal basis for the column span of $AS_{1}S_{3}$. 6. 6.

Let $p_{i}^{(4)}=\frac{\|v_{i}\|_{2}^{2}}{\|V\|_{F}^{2}}$ and $t_{4}=c_{4}\left(\frac{t_{3}\log t_{3}}{\epsilon^{2}}\right)$. Sample $S_{4}\in\mathbb{R}^{n\times t_{4}}$ where the $j^{th}$ column is set to $\frac{1}{\sqrt{t_{4}p_{i}^{(4)}}}e_{i}$ with probability $p_{i}^{(4)}$. Compute $W\in\mathbb{R}^{t_{3}\times t_{1}}$ satisfying:

[TABLE] 7. 7.

Compute an orthogonal basis $Q\in\mathbb{R}^{n\times k}$ for the column span of $AS_{1}S_{3}W$. Let $t_{5}=c_{5}\left(k\log k+\frac{k}{\epsilon}\right)$, set $p_{i}^{(5)}=\frac{\|q_{i}\|_{2}^{2}}{\|Q\|_{F}^{2}}$, where $q_{i}$ is the $i^{th}$ row of $Q$. Sample $S_{5}\in\mathbb{R}^{n\times t_{5}}$ where the $j^{th}$ column is set of $\frac{1}{\sqrt{t_{5}p_{i}^{(5)}}}e_{i}$ with probability $p_{i}^{(5)}$. Solve:

[TABLE] 8. 8.

Return $Q,N\in\mathbb{R}^{n\times k}$.

Theorem 9** (Sublinear Time Low-Rank Approximation).**

Given any PSD $A\in\mathbb{R}^{n\times n}$, for sufficiently large constants $c,c^{\prime},c_{1},c_{2},c_{3},c_{4},c_{5}$, for any $\epsilon<1$, Algorithm 1 accesses $O(\frac{n\cdot k\log^{2}n}{\epsilon^{2.5}}+\sqrt{n}k^{1.5}\cdot\log^{2}n\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon))$ entries of $A$, runs in $\tilde{O}\left(\frac{nk^{\omega-1}}{\epsilon^{2(\omega-1)}}+\sqrt{n}k^{\omega-.5}\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon)\right)$ time, and with probability at least $9/10$ outputs $M,N\in\mathbb{R}^{n\times k}$ with:

[TABLE]

For simplicity, we will show $\|A-MN^{T}\|_{F}^{2}\leq(1+O(\epsilon))\|A-A_{k}\|_{F}^{2}$. By scaling up the constants $c,c^{\prime},c_{1},c_{2},c_{3},c_{4},c_{5}$ we can then scale down $\epsilon$, achieving the claimed result. We first show:

Lemma 10**.**

For sufficiently large constants $c,c^{\prime},c_{1},c_{2}$, with probability $98/100$, Steps 1-4 of Algorithm 1 produce $Z\in\mathbb{R}^{t_{1}\times k_{1}}$ satisfying:

[TABLE]

Proof.

In Step 3, $S_{1}$ is sampled using probabilities proportional to the scores $\ell_{i}^{(1)}=\sqrt{\frac{n}{k}}\tilde{\tau}_{i}^{k}(A^{1/2})+\sqrt{\frac{n\epsilon^{4}}{k_{1}}}\tilde{\tau}_{i}^{c^{\prime}k_{1}}(A^{1/2})$. The first part of this score ensures that by Lemma 3, with high probability $AS_{1}$ is an $(\epsilon,k)$-column PCP of $A$. The second half, in combination with the construction of $S_{2}$ ensures via Lemmas 7 and 8 that if $c^{\prime}$ is large enough, $\tilde{A}=S_{2}^{T}AS_{1}$ is both a $(1/2,k_{1})$-row PCP of $AS_{1}$ for Frobenius norm error with probability $99/100$ and a $(1/2,k_{1})$-spectral PCP of $AS_{1}$ with high probability. Note that in $\ell_{i}^{(1)}$ the scores corresponding to $\tilde{\tau}_{i}^{c^{\prime}k_{1}}(A^{1/2})$ are scaled down by an $\epsilon^{2}$ factor from what is required to apply Lemmas 7 and 8 with rank $k_{1}$ and error $1/2$. However, $t_{1}$ is oversampled by a $1/\epsilon^{2}$ factor, balancing this scaling. By a union bound, all three PCP properties hold simultaneously with probability $98/100$.

By the Frobenius error PCP property of Lemma 7, letting $U_{k_{1}}$ contain the top $k_{1}$ row singular vectors of $AS_{1}$ we have:

[TABLE]

by the fact that $|\Delta|\leq 600\|A-A_{k}\|_{F}^{2}$ and $\|AS_{1}-(AS_{1})_{k_{1}}\|_{F}^{2}\leq\|AS_{1}-(AS_{1})_{k}\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}\leq 2\|A-A_{k}\|_{F}^{2}$ since $AS_{1}$ is an $(\epsilon,k)$-column PCP of $A$. Further, via Lemma 8, for any rank-$k_{1}$ projection $P$:

[TABLE]

By this bound, we have

[TABLE]

where the final bound follows because if we set $c\geq 2$, $k_{1}\geq\lceil 2k/\epsilon\rceil$ so $\|AS_{1}-(AS_{1})_{k_{1}}\|_{2}^{2}\leq\frac{\epsilon}{k}\|AS_{1}-(AS_{1})_{k}\|_{F}^{2}\leq\frac{2\epsilon}{k}\|A-A_{k}\|_{F}^{2}$. With (4) and (6) in place, applying (5) again, if we compute $Z$ satisfying the guarantee in Step 4 we have:

[TABLE]

which gives the lemma after adjusting constants on $\epsilon$ by making $c$, $c^{\prime}$, $c_{1}$ and $c_{2}$ sufficiently large. ∎

Lemma 10 ensures that the rank-$k$ matrix $W$ computed in Step 6 gives a near optimal low-rank approximation of $AS_{1}$. Specifically we have:

Lemma 11**.**

With probability $95/100$, for sufficiently large $c,c^{\prime},c_{1},c_{2},c_{3},c_{4}$, Step 5 of Algorithm 1 produces $W$ such that, letting $Q\in\mathbb{R}^{n\times k}$ be an orthonormal basis for the column span of $AS_{1}S_{3}W$:

[TABLE]

Proof.

We first consider the optimization problem:

[TABLE]

We have:

[TABLE]

where the second to last step uses that $\|AS_{1}(I-ZZ^{T})\|_{2}^{2}\leq\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}$ by Lemma 10 and the last step uses that $\|AS_{1}-(AS_{1})_{k}\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$ since $AS_{1}$ is an $(\epsilon,k)$-column PCP of $A$.

Note that since it is rank-$k$ we can write $M^{*}=Y^{*}N^{*}$ where $Y^{*}\in\mathbb{R}^{n\times k}$ and $N^{*}\in\mathbb{R}^{k\times k_{1}}$ has orthonormal rows. $Y^{*}$ is the solution to the unconstrained low-rank approximation problem:

[TABLE]

$N^{*}Z^{T}$ has orthonormal rows, so its column norms are its leverage scores. $S_{3}$ is sampled using the column norms of $Z^{T}$, which upper bound those of $N^{*}Z^{T}$. Since $\|Z\|_{F}^{2}=\lceil ck/\epsilon\rceil$, $t_{3}=\Theta\left(\|Z\|_{F}^{2}\log(\|Z\|_{F}^{2})+\|Z\|_{F}^{2}/\epsilon\right)$. It is thus well known (see e.g. Theorem 38 of [CW13], and [DMM08] for earlier work) that if we set

[TABLE]

we have with probability $99/100$:

[TABLE]

$\tilde{Y}$ can be computed in closed form as $\tilde{Y}=AS_{1}S_{3}(N^{*}Z^{T}S_{3})^{+}$, which is in the column span of $AS_{1}S_{3}$. Thus (8) demonstrates that there is some rank-$k$ $M$ in this span satisfying $\|MZ^{T}-AS_{1}\|_{F}^{2}\leq(1+\epsilon)\|M^{*}Z^{T}-AS_{1}\|_{F}^{2}\leq(1+4\epsilon)\|A-A_{k}\|_{F}^{2}$ by (4). Thus if we compute

[TABLE]

we will have

[TABLE]

We can compute $W^{*}$ approximately by sampling the rows of $AS_{1}S_{3}$ by their leverage scores. These are given by the row norms of the orthogonal basis $V$ as computed in Step 5 and sampled with to obtain $S_{4}$ in Step 6. For any $W$, by the Pythagorean theorem, since $V$ spans $AS_{1}S_{3}$,

[TABLE]

Similarly

[TABLE]

For the first term, since $S_{4}$ is sampled via the leverage scores of $AS_{1}S_{3}$ and $t_{4}=c_{4}\left(\frac{t_{3}\log t_{3}}{\epsilon^{2}}\right)$, $S_{4}$ gives a subspace embedding for the column span of $AS_{1}S_{3}$ with high probability and:

[TABLE]

For the cross term, again since $V$ spans the columns of $AS_{1}S_{3}$ we have:

[TABLE]

where the last bound follows with probability $99/100$ from a standard approximate matrix multiplication result [DKM06] since $S_{4}$ is sampled using the row norms of $V$. $\|(I-VV^{T})AS_{1}\|_{F}\leq\|AS_{1}S_{3}W^{*}Z^{T}-AS_{1}\|_{F}\leq\|AS_{1}S_{3}WZ^{T}-AS_{1}\|_{F}$ and so overall:

[TABLE]

Combining (10), (4), (12), and (13) for any $W$,

[TABLE]

where $\Delta=\|S_{4}^{T}(I-VV^{T})AS_{1}\|_{F}^{2}-\|(I-VV^{T})AS_{1}\|_{F}^{2}$ is fixed independent of $W$. Thus, if we set:

[TABLE]

with probability $95/100$ union bounding over all the probabilistic events above, and applying (9):

[TABLE]

which gives the lemma after adjusting constants on $\epsilon$ by making $c,c^{\prime},c_{1},c_{2},c_{3},c_{4}$ large enough. ∎

Proof of Theorem 9.

By Lemma 11 and since $AS_{1}$ is an $(\epsilon,k)$-column PCP of $A$, we have $\|A-QQ^{T}A\|_{F}^{2}\leq\frac{1+\epsilon}{1-\epsilon}\|A-A_{k}\|_{F}^{2}.$ In Step 7 we sample $c_{5}(k\log k+k/\epsilon)$ rows of Q by their standard leverage scores (their row norms) and so have with probability 99/100:

[TABLE]

Union bounding over failure probabilities and adjusting constants on $\epsilon$ yields the final claim.

It just remains to discuss Algorithm 1’s runtime and query complexity. In Step 3, forming $\tilde{A}=S_{2}^{T}AS_{1}$ requires reading $t_{1}\cdot t_{2}=O\left(\frac{nk\log^{2}n}{\epsilon^{2.5}}\right)$ entries of $A$. This dominates the access cost of Step 1, which requires $O(nk_{1}\log k_{1})=O\left(\frac{nk\log(k/\epsilon)}{\epsilon}\right)$ accesses by Lemma 4, forming $AS_{1}S_{3}$ in Step 5 which requires $O\left(\frac{nk\log(k/\epsilon)}{\epsilon}+\frac{nk}{\epsilon^{2}}\right)$ accesses, and forming $S_{5}^{T}A$ in Step 7 which requires $O\left(nk\log k+\frac{nk}{\epsilon}\right)$ accesses. Forming $S_{4}^{T}AS_{1}$ in Step 6 requires $t_{1}\cdot t_{4}=O\left(\frac{\sqrt{n}k^{1.5}\log(k/\epsilon)\log n}{\epsilon^{6}}\right)$ accesses, which when $n$ is large compared to $k$ and $1/\epsilon$ will be dominated by our linear in $n$ terms.

For time complexity, Step 1 requires $O(n(k_{1}\log k_{1})^{\omega-1})=\tilde{O}\left(\frac{nk^{\omega-1}}{\epsilon^{\omega-1}}\right)$ time. Computing $Z$ in Step 4 can be done using an input sparsity time algorithm for spectral norm error with rank $k_{1}$ and error parameter $\epsilon^{\prime}=\Theta(1)$. By Theorem 27 of [CEM*+*15] or using input sparsity time ridge leverage score sampling [CMM17] in conjunction with the spectral norm PCP result of Lemma 24, the total runtime required is $O(\operatorname{nnz}(\tilde{A}))+\tilde{O}(\frac{\sqrt{n}k^{\omega-1}}{\epsilon^{\omega-1}})$.

We can compute $V$ in Step 5 in $O(nt_{3}^{\omega-1})=\tilde{O}\left(\frac{nk^{\omega-1}}{\epsilon^{2(\omega-1)}}\right)$. In Step 6, we can compute $W$ by first multiplying $S_{4}^{T}AS_{1}$ by $Z$ and then multiplying by $(S_{4}^{T}AS_{1}S_{3})^{+}$ and taking the best rank-$k$ approximation of the result. The total runtime is $\tilde{O}\left(t_{1}\cdot k^{\omega-1}+k^{\omega}\right)\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon)=\tilde{O}(\sqrt{n}k^{\omega-.5})\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon)$. Finally, we can compute $Q$ in Step 7 by first computing a $t_{3}\times k$ span of the column space of $W$, multiplying this by $AS_{1}S_{3}$, and then taking an orthogonal basis for the result, requiring total time $\tilde{O}(k^{\omega}\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon))+O\left(\frac{nk^{\omega-1}\log(k/\epsilon)}{\epsilon^{2}}+nk^{\omega-1}\right)$. The final regression problem can be solved by forming the pseudoinverse of $S_{5}^{T}Q$ in $O(k^{\omega}(\log k+1/\epsilon))$ time and applying it to $S_{5}^{T}A$ in $O(nk^{\omega-1}(\log k+1/\epsilon))$ time. So overall our runtime cost with linear dependence on $n$ is dominated by the cost of Step 5. We additionally must add the $\tilde{O}(\sqrt{n}k^{\omega-.5})\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon)$ term from Step 6, which only dominates if $n$ is relatively small. ∎

In many applications it is desirable that the low-rank approximation to $A$ is also symmetric and positive semidefinite. We show in Appendix C that a modification to Algorithm 1 can satisfy this constraint also in $\tilde{O}(n\mathop{\mathrm{missing}}{poly}(k/\epsilon))$ time. The upshot is:

Theorem 12** (Sublinear Time Low-Rank Approximation – PSD Output).**

There is an algorithm that given any PSD $A\in\mathbb{R}^{n\times n}$, accesses $\tilde{O}\left(\frac{nk^{2}}{\epsilon^{2}}+\frac{nk}{\epsilon^{3}}\right)$ entries of $A$, runs in $\tilde{O}\left(\frac{nk^{\omega}}{\epsilon^{\omega}}+\frac{nk^{\omega-1}}{\epsilon^{3(\omega-1)}}\right)$ time and with probability at least $9/10$ outputs $M\in\mathbb{R}^{n\times k}$ with:

[TABLE]

5 Query Lower Bound

We now present our lower bound on the number of accesses to $A$ required to compute a near-optimal low-rank approximation, matching the query complexity of Algorithm 4 up to a $\tilde{O}(1/\epsilon^{1.5})$ factor.

Theorem 13**.**

Assume that $k,\epsilon$ are such that $nk/\epsilon=o(n^{2})$. Consider any (possibly randomized) algorithm $\mathcal{A}$ that, given any PSD $A\in\mathbb{R}^{n\times n}$, outputs a $(1+\epsilon)$-approximate rank-$k$ approximation to $A$ (in the Frobenius norm) with probability at least $2/3$. Then there must be some input ${A}$ of which $\mathcal{A}$ reads at least $\Omega(nk/\epsilon)$ positions, possibly adaptively, in expectation.

5.1 Lower Bound Approach

We prove Theorem 13 via Yao’s minimax principle [Yao77], proving a lower bound for randomized algorithms on worst-case inputs via a lower bound for deterministic algorithms on a hard input distribution. Specifically, we will draw the input $A$ from a distribution over binary matrices. $A$ has all $1$’s on its diagonal, along with $k$ randomly positioned (non-contiguous) blocks of all $1$’s, each of size $\sqrt{2\epsilon n/k}\times\sqrt{2\epsilon n/k}$. In other words, $A$ is the adjacency matrix (plus identity) of a graph with $k$ cliques of size $\sqrt{2\epsilon n/k}$, placed on random subsets of the vertices, with all other vertices isolated.

It is easy to see that every $A$ chosen from the above distribution is PSD since applying a permutation yields a block diagonal matrix, each of whose blocks is a PSD matrix (either a single $1$ entry or a rank-$1$ all $1$’s block). Additionally, for every ${A}$ chosen from the distribution, the optimal rank-$k$ approximation to $A$ projects off each of the $k$ blocks, achieving Frobenius norm error $\|A-A_{k}\|_{F}^{2}=n-k\sqrt{2\epsilon n/k}\approx n$. To match this up to a $1+\epsilon$ factor, any near-optimal rank-$k$ approximation must at least capture a constant fraction of the Frobenius norm mass in the blocks since this mass is $k\cdot 2\epsilon n/k=2\epsilon n$.

Doing so requires identifying at least a constant fraction of the blocks. However, since block positions are chosen uniformly at random, and since the diagonal entries of $A$ are identical and so convey no information about the positions, to identify a single block, any algorithm essentially must read arbitrary off-diagonal entries until it finds a $1$. There are $\approx n^{2}$ off-diagonal entries with just $2\epsilon n$ of them $1$’s, so identifying a first block requires $\Omega(n/\epsilon)$ queries to $A$ in expectation (over the random choice of ${A}$). Since the vast majority of vertices are isolated and not contained within a block, finding this first block does little to make finding future blocks easier. So overall, the algorithm must make $\Omega(nk/\epsilon)$ queries in expectation to find a constant fraction of the $k$ blocks and output a near-optimal low-rank approximation with good probability.

While the above intuition is the key idea behind the lower bound, a rigorous proof requires a number of additional steps and modifications, detailed in the remainder of this section. In Section 5.2 we introduce the notion of a primitive approximation to a matrix, which we will employ in our lower bound for low-rank approximation. In Section 5.3 we show that any deterministic low-rank approximation algorithm that succeeds with good probability on the input distribution described above can be used to give an algorithm that computes a primitive approximation with good probability on a matrix drawn from a related input distribution (Lemma 20). This reduction yields a lower bound for deterministic low-rank approximation algorithms (Theorem 23), which gives Theorem 13 after an application of Yao’s principle.

5.2 Primitive Approximation

We first define the notion of an $\epsilon$-primitive approximation to a matrix and establish some basic properties of these approximations.

Definition 14** ($\epsilon$-primitive Approximation).**

A matrix $A^{\prime}\in\mathbb{R}^{m\times m}$ is said to be $\epsilon$-primitive for $A\in\mathbb{R}^{m\times m}$ if the squared Frobenius norm of $A-A^{\prime}$, restricted to its off-diagonal entries is $\leq\epsilon m$. Note that $A^{\prime}$ is allowed to have any rank.

We also define a distribution on matrices with a randomly placed block of all one entries that will will later use in our hard input distribution construction.

Definition 15** (Random Block Matrix).**

For any $m,\epsilon$ with $m/\epsilon\leq m^{2}$, let $\mu(m,\epsilon)$ be the distribution on ${A}\in\mathbb{R}^{m\times m}$ defined as follows. We choose a uniformly random subset $S$ of $[m]$ where $|S|=\sqrt{16\epsilon m}$, where we assume for simplicity that $|S|$ is an integer. We generate a random matrix $A$ by setting for each $i\neq j\in S$, $A_{i,j}=1$. We then set $A_{i,i}=1$ for all $i$ and set all remaining entries of $A$ to equal [math].

Note that we associate a random subset $S$ with the sampling of a matrix $A$ according to $\mu(m,\epsilon)$. It is clear that any $A$ in the support of $\mu(m,\epsilon)$ is PSD. This is since any ${A}$ in the support of $\mu(m,\epsilon)$, after a permutation, is composed of an $|S|\times|S|$ all ones block and an $(m-|S|)\times(m-|S|)$ identity.

If ${A}^{\prime}$ is $\epsilon$-primitive for ${A}$ in the support of $\mu(m,\epsilon)$, it approximates ${A}$ up to error $\epsilon m$ on its off-diagonal entries. Let $R$ denote the set of off-diagonal entries of ${A}$ restricted to the intersection of the rows and columns indexed by $S$. The error on the entries of $R$ is at most $\epsilon m$. Further, restricted to these entries, ${A}$ is an all ones matrix. Thus, on the entries of $R$, both ${A}$ and ${A}^{\prime}$ look far from an identity matrix. Formally, we show:

Lemma 16**.**

If $A^{\prime}\in\mathbb{R}^{m\times m}$ is $\epsilon$-primitive for an $A\in\mathbb{R}^{m\times m}$ in the support of $\mu(m,\epsilon)$, then $A^{\prime}$ is not $\epsilon$-primitive for $I$.

Proof.

By definition, any $\epsilon$-primitive matrix ${A}^{\prime}$ for ${A}$ has squared Frobenius norm restricted to ${A}$’s off-diagonal entries of $\leq\epsilon m$. Let $R$ denote the set of off-diagonal entries in the intersection of the rows and columns of ${A}$ indexed by $S$. Assuming $|S|\geq 4$, which follows from the assumption in Definition 15 that $m/\epsilon\leq m^{2}$, $|R|=|S|^{2}-S\geq 12\epsilon m$.

Restricted to the entries in $R$, ${A}$ is an all ones matrix. Thus, on at least $8\epsilon m$ of these entries ${A}^{\prime}$ must have value at least $1/2$. Otherwise it would have greater than $|R|-4\epsilon m\geq 4\epsilon m$ entries with value $\leq 1/2$ and so squared Frobenius norm error on ${A}$’s off-diagonal entries greater than $\frac{1}{2^{2}}\cdot 4\epsilon m\geq\epsilon m$, contradicting the assumption that it is $\epsilon$-primitive for ${A}$.

Restricted to the entries in $R$, the identity matrix is all zero. Thus, ${A}^{\prime}$ has Frobenius norm error on these entries at least $\frac{1}{2^{2}}\cdot 8\epsilon m\geq 2\epsilon m$. Thus ${A}^{\prime}$ is not $\epsilon$-primitive for ${I}$, giving the lemma. ∎

Consider ${A}$ drawn from $\mu(m,\epsilon)$ with probability $1/2$ and set to ${I}$ with probability $1/2$. By Lemma 16, any (possibly randomized) algorithm $\mathcal{A}$ that returns an $\epsilon$-primitive approximation to ${A}$ with good probability can be used to distinguish with good probability whether ${A}$ is in the support of $\mu(m,\epsilon)$ or ${A}={I}$. This is because, by Lemma 16, the output of such an algorithm can only be correct either for some ${A}$ in the support of $\mu(m,\epsilon)$ or for the identity. We first define the distribution over matrices to which this result applies:

Definition 17**.**

For any $m,\epsilon$ with $m/\epsilon\leq m^{2}$, let $\gamma(m,\epsilon)$ be the distribution on ${A}\in\mathbb{R}^{m\times m}$ defined as follows. With probability $1/2$ draw ${A}$ from $\mu(m,\epsilon)$ (Definition 15). Otherwise, draw ${A}$ from $\nu(m)$, which is the distribution whose support is only the $m\times m$ identity matrix.

Note that, as with $\mu(m,\epsilon)$, we associate a random subset $S$ with $|S|=\sqrt{16\epsilon m}$ with the sampling of a matrix ${A}$ according to $\gamma(m,\epsilon)$. $S$ is not used in the case that $A$ is drawn from $\nu(m)$.

Formally, since $\mathcal{A}$ can distinguish whether ${A}$ is in the support of $\mu(m,\epsilon)$ or $\nu(m)$ (i.e., ${A}={I}$) with good probability we can prove that the distribution of $\mathcal{A}$’s access pattern (over randomness in the input and possible randomization in the algorithm) is significantly different when it is given input ${A}$ drawn from $\mu(m,\epsilon)$ than when it is given ${A}$ drawn from $\nu(m)$. Recall for distributions $\alpha$ and $\beta$ supported on elements $s$ of a finite set $S$, that the total variation distance $D_{TV}(\alpha,\beta)=\sum_{s\in S}|\alpha(s)-\beta(s)|$, where $\alpha(s)$ is the probability of $s$ in distribution $\alpha$.

Corollary 18**.**

Suppose that a (possibly randomized) algorithm $\mathcal{A}$, with probability at least $7/12$ over its random coin flips and random input $A\in\mathbb{R}^{n\times n}$ drawn from $\gamma(m,\epsilon)$, outputs an $\epsilon$-primitive matrix for $A$. Further, suppose that $\mathcal{A}$ reads at most $r$ positions of $A$, possibly adaptively.333That is, for any input ${A}$, in any random execution, $\mathcal{A}$ reads at most $r$ entries of ${A}$. Let $S$ be a random variable indicating the list of positions read and their corresponding values.444$S$ is determined by the random input ${A}$ and the random choices of $\mathcal{A}$. Since $\mathcal{A}$ reads at most $r$ positions of any input, we always have $|S|\leq r$. Let $L(\mu)$ denote the distribution of $S$ conditioned on $A\sim\mu(m,\epsilon)$, and let $L(\nu)$ denote the distribution of $S$ conditioned on ${A}\sim\nu(m)$.555Here and throughout, we let $A\sim\mu(m,\epsilon)$ denote the event that, when ${A}$ is drawn from the distribution $\gamma(m,\epsilon)$ (Definition 17), which is a mixture of the distributions $\mu(m,\epsilon)$ and $\nu(m)$, that ${A}$ is drawn from $\mu(m,\epsilon)$. ${A}\sim\nu(m)$ denotes the analogous event for $\nu(m)$. Then

[TABLE]

Proof.

By Lemma 16, if algorithm $\mathcal{A}$ succeeds then its output can be used to decide if $A\sim\mu(m,\epsilon)$ or if $A\sim\nu(m)$. The success probability of any such algorithm is well-known (see, e.g., Proposition 2.58 of [BY02]) to be at most $1/2+D_{TV}(L(\mu),L(\nu))/2$. Making this quantity at least $7/12$ and solving for $D_{TV}(L(\mu),L(\nu))$ proves the corollary. ∎

5.3 Lower Bound for Low-Rank Approximation

We now give a reduction, showing that any deterministic relative error $k$-rank approximation algorithm that succeeds with good probability on the hard input distribution described in Section 5.1 can be used to compute an $\epsilon$-primitive approximation to ${A}$ that is drawn from $\gamma(n/(2k),\epsilon)$ with good probability. We first formally define this input distribution.

Definition 19** (Hard Input Distribution – Low-Rank Approximation).**

Suppose $2nk/\epsilon\leq n^{2}$. Let $\gamma_{b}$ be the distribution on ${A}\in\mathbb{R}^{n\times n}$ determined as follows. We draw a uniformly random subset $S$ of $[n]$ where $|S|=n/2$, where we assume for simplicity that $|S|$ is an integer. We further partition $S$ into $k$ subsets $S^{1},S^{2},...,S^{k}$ chosen uniformly at random. For all $\ell\in[k]$, $|S^{\ell}|=n/(2k)$, which we also assume to be an integer.

Letting ${A}^{\ell}$ denote the entries of ${A}$ restricted to the intersection of the rows and columns indexed by $S^{\ell}$, we independently draw each ${A}^{\ell}$ from $\gamma(n/(2k),\epsilon)$ (Definition 17).666By the assumption that $2nk/\epsilon\leq n^{2}$ we have $\frac{n}{2k\epsilon}\leq\left(\frac{n}{2k}\right)^{2}$ and so this is a valid setting of the parameters for $\gamma(m,\epsilon)$. We then set ${A}_{i,i}=1$ for all $i$ and set all remaining entries of $A$ to equal [math].

Our reduction from $\epsilon$-primitive approximation to low-rank approximation is as follows:

Lemma 20** (Reduction from Primitive Approximation to Low-Rank Approximation).**

Suppose that $nk/\epsilon=o(n^{2})$ and that $\mathcal{A}$ is a deterministic algorithm that, with probability $\geq 2/3$ on random input ${A}\in\mathbb{R}^{n\times n}$ drawn from $\gamma_{b}$, outputs a $(1+\epsilon/26)$-approximate rank-$k$ approximation to ${A}$. Further suppose $\mathcal{A}$ reads at most $r$ positions of ${A}$, possibly adaptively.

Then there is a randomized algorithm $\mathcal{B}$ that, with probability $\geq 7/12$ over its random coin flips and random input ${B}$ drawn from $\gamma(n/(2k),\epsilon)$, outputs an $\epsilon$-primitive matrix for ${B}$. Further, letting $L(\mu),L(\nu)$ be as defined in Corollary 18 for $\mathcal{B}$,

[TABLE]

Proof.

Consider a randomized algorithm $\mathcal{B}$ that, given ${B}$ drawn from $\gamma(n/(2k),\epsilon)$ generates a random matrix ${A}^{n\times n}$ drawn from $\gamma_{b}$ as follows. Choose a uniformly random subset $S$ of $[n]$ with $|S|=n/2$. Partition $S$ into $k$ subsets $S^{1},S^{2},...,S^{k}$ chosen uniformly at random. Note that for $\ell\in[k]$ $|S^{\ell}|=n/(2k)$. Letting ${A}^{\ell}$ denote the entries of ${A}$ restricted to the intersection of the rows and columns indexed by $S^{\ell}$, set ${A}^{1}={B}$, and for $\ell=2,...,k$ independently draw ${A}^{\ell}$ from $\gamma(n/(2k),\epsilon)$. Set ${A}_{i,i}=1$ for all $i$ and set all remaining entries of $A$ equal to [math].

After generating ${A}$, $\mathcal{B}$ then applies $\mathcal{A}$ to ${A}$ to compute a rank-k approximation ${A}^{\prime}$. $\mathcal{B}$ then outputs $({A}^{\prime})^{1}$, the $n/(2k)\times n/(2k)$ submatrix of ${A}^{\prime}$ corresponding to the intersection of the rows and columns indexed by $S^{1}$. We have the following:

Claim 21**.**

With probability $\geq 7/12$ over the random choices of $\mathcal{B}$ and over the random input ${B}$ drawn from $\gamma(n/(2k),\epsilon)$, $({A}^{\prime})^{1}$ is $\epsilon$-primitive for ${B}$.

Proof.

It is clear that ${A}$ generated by $\mathcal{B}$ is distributed according to $\gamma_{b}$ (Definition 19). By construction, for any $A\sim\gamma_{b}$ there is a rank-$k$ approximation of cost at most $n$; indeed this follows by choosing the best rank-$1$ solution for each ${A}^{\ell}$. Consequently if $\mathcal{A}$ succeeds on ${A}$, then its output $A^{\prime}$ satisfies $\|A-A^{\prime}\|_{F}^{2}\leq n+\epsilon n/26$.

Note that ${A}$ has ${A}_{i,i}=1$ for all $i\in[n]$. So the squared Frobenius norm cost of any rank-$k$ approximation restricted to the diagonal is at least $n-k$. Let $c_{\ell}$ be the squared Frobenius norm cost of ${A}^{\prime}$ restricted to the off diagonal entries in ${A}^{\ell}$. Note that $c_{\ell}$ is a random variable. Then,

[TABLE]

By averaging for at least a $11/12$ fraction of the blocks $i$,

[TABLE]

assuming $\epsilon n/(26k)\geq 8$, which holds if $\epsilon n/k=\omega(1)$, which follows from our assumption that $nk/\epsilon=o(n^{2})$.

It follows by symmetry of $\gamma_{b}$ with respect to the $k$ blocks that with probability at least $11/12$, if $\mathcal{A}$ succeeds on ${A}$, $c_{1}\leq\epsilon n/(2k)$. This gives that $({A}^{\prime})^{1}$ is $\epsilon$-primitive (Definition 14) for ${A}^{1}={B}$. This yields the claim after applying a union bound, since $\mathcal{A}$ succeeds with probability at least $2/3$. ∎

It remains to bound the total variation distance between $\mathcal{B}$’s access pattern when ${B}\sim\mu(n/(2k),\epsilon)$ and when ${B}\sim\nu(n/(2k))$. We have:

Claim 22**.**

For the algorithm $\mathcal{B}$ defined above

[TABLE]

Proof.

Let $\mathcal{W}$ denote the random variable that encompasses $\mathcal{B}$’s random choices in choosing the indices in $S^{2},...,S^{k}$ and in setting the entries in ${A}^{2},...,{A}^{k}$. Let $\Omega$ denote the set of all possible values of $\mathcal{W}$. Let $S$ denote the random subset with $|S|=\sqrt{n/(2k)\cdot\epsilon}$ associated with the drawing of ${B}$ from $\gamma(n/(2k),\epsilon)$ (see Definition 17). Let $\bar{S}$ denote the set of off-diagonal entries of ${A}$ in the intersection of the rows and columns indexed by $S$.

If ${B}\sim\mu(n/(2k),\epsilon)$ then all entries of $\bar{S}$ are $1$. If ${B}\sim\nu(n/(2k))$ then ${B}={I}$ and so all entries of $\bar{S}$ are [math]. Outside of the entries in $\bar{S}$, the entries of ${B}$ are identical in the two cases that ${B}\sim\mu(n/(2k),\epsilon)$ and ${B}\sim\nu(n/(2k))$ (in particular, they are all [math] off the diagonal and $1$ on the diagonal).

So, for any $w\in\Omega$, conditioned on $\mathcal{W}=w$, all entries outside $\bar{S}$ are fixed. Further conditioned on ${B}\sim\nu(n/(2k))$, all entries in $\bar{S}$ are [math]. So all entries of ${A}$ are fixed and $\mathcal{A}$ always reads the same sequence of entries $(i_{w,1},j_{w,1}),...,(i_{w,r},j_{w,r})$. Further, for any $w\in\Omega$, conditioned on $\mathcal{W}=w$, $S$ is a uniform random subset of $R=[n]\setminus(S^{2}\cup...\cup S^{k})$. $|R|=n-(k-1)\cdot n/(2k)\geq n/2$. So, for any $\ell\in[r]$ we have:

[TABLE]

By a union bound, letting $\mathcal{E}$ be the event that $\mathcal{A}$ reads $(i_{w,1},j_{w,1}),...,(i_{w,r},j_{w,r})$ and does not read any entry of $\bar{S}$ in its $r$ accesses to ${A}$ we have for any $w\in\Omega$:

[TABLE]

Thus, for any $w\in\Omega$, regardless of whether ${B}\sim\nu(n/(2k))$ or ${B}\sim\mu(n/(2k),\epsilon)$, $\mathcal{A}$ has access pattern $(i_{w,1},j_{w,1}),...,(i_{w,r},j_{w,r})$ to ${A}$ with probability $\geq 1-\frac{2\epsilon r}{nk}$. Correspondingly, $\mathcal{B}$ has a fixed access pattern to ${B}$ with probability $\geq 1-\frac{2\epsilon r}{nk}$. Thus,

[TABLE]

yielding the claim. ∎

In combination Claims 21 and 22 give Lemma 20. ∎

We can now use Lemma 20 to argue that if a deterministic low-rank approximation algorithm succeeding with good probability on a random input drawn from $\gamma_{b}$ accesses too few entries, then it can be used to give a primitive approximation algorithm violating Corollary 18.

Theorem 23** (Lower Bound for Deterministic Algorithms).**

Assume that $n,k,\epsilon$ are such that $nk/\epsilon=o(n^{2})$. Consider any deterministic algorithm $\mathcal{A}$ that, given random input ${A}$ drawn from $\gamma_{b}$, outputs a $(1+\epsilon)$-approximate rank-$k$ approximation to $A$ (in the Frobenius norm) with probability at least $2/3$. Further, suppose $\mathcal{A}$ reads at most $r$ positions of $A$, possibly adaptively. Then $r=\Omega(nk/\epsilon)$.

Proof.

Assume towards a contradiction that $r=o(nk/\epsilon)$. Then applying Lemma 20, $\mathcal{A}$ can be used to give a randomized algorithm $\mathcal{B}$ that with probability $\geq 7/12$ over its random coin flips and random input ${B}\in n/(2k)\times n/(2k)$ drawn from $\gamma(n/(2k),\epsilon)$, outputs an $\epsilon$-primitive matrix for ${B}$. Further, letting $L(\mu),L(\nu)$ be as defined in Corollary 18 for $\mathcal{B}$, by Lemma 20, $D_{TV}(L(\mu),L(\nu))\leq\frac{2\epsilon r}{nk}.$ For $r=o(nk/\epsilon)$, $\frac{2\epsilon r}{nk}=o(1)$, contradicting Corollary 18, and giving the theorem. ∎

Theorem 13 follows directly from applying Yao’s minimax principle ([Yao77], Theorem 3) to Theorem 23.

6 Spectral Norm Error Bounds

We conclude by showing how to modify Algorithm 1 to output a low-rank approximation $B$ achieving the spectral norm guarantee:

[TABLE]

This can be significantly stronger than the Frobenius guarantee (1) when $\|A-A_{k}\|_{F}^{2}$ is large, and, for example is critical in our application to ridge regression, discussed in Section 6.2.

It is not hard to see that since additive error in the Frobenius norm upper bounds additive error in the spectral norm (see e.g. Theorem 3.4 of [Gu14]) that for $B$ satisfying the Frobenius norm guarantee $\|A-B\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$, we immediately have the spectral bound $\|A-B\|_{2}^{2}\leq\|A-A_{k}\|_{2}^{2}+\epsilon\|A-A_{k}\|_{F}^{2}$. Thus, we can achieve 14 simply by running Algorithm 1 with error parameter $\epsilon/k$. However, this approach is suboptimal. Applying Theorem 9, our query complexity would be $\Theta\left(\frac{nk^{3.5}\log^{2}n}{\epsilon^{2.5}}\right)$. We improve this $k$ dependence significantly in Algorithm 2. Since (14) is often applied (see for example Section 6.2) with $k^{\prime}=k/\epsilon$ and $\epsilon=\Theta(1)$ to give $\|A-B\|_{2}^{2}\leq O\left(\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}\right)$, optimizing $k$ dependence is especially important.

We first give an extension of Lemma 3 to the spectral norm case. This lemma provides the column sampling analog to Lemma 8.

Lemma 24** (Spectral Norm PCP).**

For any $A\in\mathbb{R}^{n\times d}$, for $i\in\{1,\ldots,d\}$, let $\tilde{\tau}_{i}^{k}\geq\tau_{i}^{k}(A)$ be an overestimate for the $i^{th}$ rank-$k$ ridge leverage score. Let $p_{i}=\frac{\tilde{\tau}^{k}_{i}}{\sum_{i}\tilde{\tau}^{k}_{i}}$ and $t=\frac{c\log(k/\delta)}{\epsilon^{2}}\sum_{i}\tilde{\tau}^{k}_{i}$ for any $\epsilon<1$ and sufficiently large constant $c$. Construct $C$ by sampling $t$ columns of $A$, each set to $\frac{1}{\sqrt{tp_{i}}}a_{i}$ with probability $p_{i}$. With probability $1-\delta$, for any orthogonal projection $P\in\mathbb{R}^{n\times n}$,

[TABLE]

We refer to $C$ as an $(\epsilon,k)$-spectral PCP of $A$.

Proof.

This follows from Corollary 34. With probability $\geq 1-\delta$, sampling by the rank-$k$ ridge leverage scores gives $C$ satisfying:

[TABLE]

We can write for any $M$, $\|M\|_{2}^{2}=\max_{x:\|x\|_{2}^{2}=1}x^{T}Mx$. When $\|x\|_{2}^{2}=1$, $\|(I-P)x\|_{2}^{2}\leq 1$ so by (15) we have for any unit norm $x$:

[TABLE]

which gives $\|C-PC\|_{2}^{2}\leq\frac{1}{1-\epsilon}\|A-PA\|_{2}^{2}+\frac{\epsilon}{(1-\epsilon)k}\|A-A_{k}\|_{F}^{2}$. Similarly we have:

[TABLE]

which gives $\frac{1}{1+\epsilon}\|A-PA\|_{2}^{2}-\frac{\epsilon}{(1+\epsilon)k}\|A-A_{k}\|_{F}^{2}\leq\|C-PC\|_{2}^{2}$. The lemma follows from combining these upper and lower bounds after adjusting constant factors on $\epsilon$ (by making $c$ large enough). ∎

6.1 Spectral Norm Low-Rank Approximation Algorithm

We now use Lemma 24, along with its row sampling counterpart, Lemma 8, to give an algorithm (Algorithm 2) for computing a near optimal spectral norm low-rank approximation to $A$.

In Steps 1-3 we sample both rows and columns of $A$ via the rank $\Theta(k/\epsilon)$ ridge leverage scores of $A^{1/2}$, ensuring with high probability that $AS_{1}$ is an $(\epsilon,k)$-spectral PCP of $A$ and $\tilde{A}$ is in turn an $(\epsilon,k)$-spectral row PCP of $AS_{1}$. Thus, if we compute (using an input sparsity time algorithm) a span $Z$ which gives a near optimal spectral norm low-rank approximation to $\tilde{A}$ (Step 3), this span will also be nearly optimal for $AS_{1}$. Since we cannot afford to fully read $AS_{1}$, we approximately project it to $Z$ by further sampling its columns using $Z$’s leverage scores (Step 4). We use leverage score sampling again in Step 5 to approximately project $A$ to the span of the result. This yields our final approximation, using the fact that $AS_{1}$ is a spectral PCP for $A$.

Algorithm 2

PSD Low-Rank Approximation – Spectral Error

Let $k_{1}=\lceil ck/\epsilon^{2}\rceil$. For all $i\in[1,..,n]$ compute $\tilde{\tau}_{i}^{k_{1}}(A^{1/2})$ which is a constant factor approximation to the ridge leverage score $\tau_{i}^{k_{1}}(A^{1/2})$. 2. 2.

Set $\ell^{(1)}_{i}=4\epsilon\sqrt{\frac{n}{k}}\tilde{\tau}_{i}^{k_{1}}(A^{1/2})$. Set $p_{i}^{(1)}=\frac{\ell_{i}^{(1)}}{\sum_{i}\ell_{i}^{(1)}}$ and $t_{1}=\frac{c_{1}\log n}{\epsilon^{2}}\sum_{i}\ell_{i}^{(1)}$. Sample $S_{1},S_{2}\in\mathbb{R}^{n\times t_{1}}$ each whose $j^{th}$ column is set to $\frac{1}{\sqrt{tp^{(1)}_{i}}}e_{i}$ with probability $p^{(1)}_{i}$. 3. 3.

Let $\tilde{A}=S_{2}^{T}AS_{1}$, and use an input sparsity time algorithm to compute orthonormal $Z\in\mathbb{R}^{t_{1}\times k}$ satisfying the Frobenius guarantee $\|\tilde{A}-\tilde{A}ZZ^{T}\|_{F}^{2}\leq 2\|\tilde{A}-\tilde{A}_{k}\|_{F}^{2}$ along with the spectral guarantee $\|\tilde{A}-\tilde{A}ZZ^{T}\|_{2}^{2}\leq(1+\epsilon)\|\tilde{A}-\tilde{A}_{k}\|_{2}^{2}+\frac{\epsilon}{k}\|\tilde{A}-\tilde{A}_{k}\|_{F}^{2}$. 4. 4.

Let $t_{3}=c_{3}\left(k\log k+\frac{k^{2}}{\epsilon}\right)$, set $p_{i}^{(3)}=\frac{\|z_{i}\|_{2}^{2}}{\|Z\|_{F}^{2}}$, and sample $S_{3}\in\mathbb{R}^{t_{1}\times t_{3}}$ where the $j^{th}$ column is set to $\frac{1}{\sqrt{t_{3}p_{i}^{(3)}}}e_{i}$ with probability $p_{i}^{(3)}$. Solve:

[TABLE] 5. 5.

Compute an orthogonal basis $Q\in\mathbb{R}^{n\times k}$ for the column span of $M$. Let $t_{4}=c_{4}\left(k\log k+\frac{k^{2}}{\epsilon}\right)$, set $p_{i}^{(4)}=\frac{\|q_{i}\|_{2}^{2}}{\|Q\|_{F}^{2}}$, and sample $S_{4}\in\mathbb{R}^{n\times t_{4}}$ where the $j^{th}$ column is set to $\frac{1}{\sqrt{t_{4}p_{i}^{(4)}}}e_{i}$ with probability $p_{i}^{(4)}$. Solve:

[TABLE] 6. 6.

Return $Q,N\in\mathbb{R}^{n\times k}$.

Theorem 25** (Sublinear Time Low-Rank Approximation –Spectral Norm Error).**

Given any PSD $A\in\mathbb{R}^{n\times n}$, for sufficiently large constants $c,c_{1},c_{2},c_{3},c_{4}$, Algorithm 2 accesses $O(\frac{n\cdot k\log^{2}n}{\epsilon^{6}}+\frac{nk^{2}}{\epsilon})$ entries of $A$, runs in $\tilde{O}\left(\frac{nk^{\omega}}{\epsilon}+\frac{nk}{\epsilon^{6}}+(\sqrt{n}k^{\omega-1}+k^{\omega+1})\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon)\right)$ time and with probability at least $9/10$ outputs $M,N\in\mathbb{R}^{n\times k}$ with:

[TABLE]

Proof.

$\ell_{i}^{(1)}=4\epsilon\sqrt{\frac{n}{k}}\tilde{\tau}_{i}^{k_{1}}(A^{1/2})$ which by Lemma 5 is within a constant factor of upper bounding the rank $k_{1}=\lceil ck/\epsilon^{2}\rceil$ ridge leverage scores of $A$. As long as $c\geq 1$, these scores upper bound the rank-$k$ ridge leverage scores. So by Lemma 24, with high probability $AS_{1}$ an $(\epsilon,k)$-spectral PCP of $A$ as long as $c_{1}$ is set large enough. Additionally, by Lemma 8, with high probability $\tilde{A}$ is an $(\epsilon,k)$-spectral row PCP of $AS_{1}$.

Frobenius norm PCP bounds also hold. By Lemma 3, $AS_{1}$ is an $(\epsilon,k_{1})$-column PCP of $A$ with high probability and $\tilde{A}$ is an $(\epsilon,k_{1})$-row PCP for $AS_{1}$ with probability $99/100$. These bounds trivially give that $AS_{1}$ and $\tilde{A}$ are $(\epsilon,k)$ PCPs for $A$ and $AS_{1}$ respectively, which gives:

[TABLE]

for some constant $c_{5}$. By (16) and the fact that $\tilde{A}$ is an $(\epsilon,k)$-spectral PCP of $AS_{1}$, for $Z$ computed in Step 3 of Algorithm 2 we have:

[TABLE]

where we assume without loss of generality that $\epsilon<1/2$ as this can be achieved by scaling our constants sufficiently.

Now, by standard approximate regression, since $S_{3}$ is sampled by the leverage scores of $Z$, and since $t_{3}=\Theta\left(k\log k+k/\epsilon^{\prime}\right)$ for $\epsilon^{\prime}=\epsilon/k$, for $M$ computed in Step 4, with probability $99/100$:

[TABLE]

This Frobenius norm bound also implies a spectral norm bound. Specifically, we can write:

[TABLE]

so by (18) we must have $\|AS_{1}ZZ^{T}-MZ^{T}\|_{F}^{2}\leq\frac{\epsilon}{k}\|AS_{1}(I-ZZ^{T})\|_{F}^{2}$. This gives:

[TABLE]

Note that $\|AS_{1}(I-ZZ^{T})\|_{F}^{2}=O(\|A-A_{k}\|_{F}^{2})$ by the Frobenius norm guarantee required in Step 3, and the fact that $\tilde{A}$ and $AS_{1}$ are $(\epsilon,k)$ PCPs for $AS_{1}$ and $A$ respectively. So combining with (16) the above implies that for $Q$ spanning the columns of $M$,

[TABLE]

Since $AS_{1}$ is an $(\epsilon,k)$-spectral PCP for $A$ we also have $\|A-QQ^{T}A\|_{2}^{2}\leq(1+O(\epsilon))\|AS_{1}-(AS_{1})_{k}\|_{2}^{2}+O\left(\frac{\epsilon}{k}\right)\|A-A_{k}\|_{F}^{2}.$ The theorem follows by applying an identical approximate regression argument for $N$ computed in Step 5 to show that

[TABLE]

with probability $99/100$. Adjusting constants on $\epsilon$ and union bounding over failure probabilities yields the final bound.

It just remains to discuss runtime and sample complexity. Constructing $S_{2}^{T}AS_{1}$ in Step 3 requires reading $t_{1}^{2}=O\left(\frac{nk\log^{2}n}{\epsilon^{6}}\right)$ entries of $A$. Constructing $AS_{1}S_{3}$ in Step 4 and $AS_{4}$ in Step 5 both require reading $O\left(nk\log k+\frac{nk^{2}}{\epsilon}\right)$ entries.

For runtime, computing $Z$ in Step 3 requires $O(\operatorname{nnz}(\tilde{A}))+\tilde{O}(\sqrt{n}k^{\omega-1}\cdot\mathop{\mathrm{missing}}{poly}(1/\epsilon))$ time using an input sparsity time algorithm (e.g. by Theorem 27 of [CEM*+*15] or using input sparsity time ridge leverage score sampling [CMM17] in conjunction with the spectral norm PCP result of Lemma 24). Computing $M$ in Step 4 and $N$ in Step 5 both require computing the pseudoinverse of a $k\times O\left(k\log k+\frac{k^{2}}{\epsilon}\right)$ matrix in $\tilde{O}(k^{\omega+1}\mathop{\mathrm{missing}}{poly}(1/\epsilon))$ time, and then applying this to an $n\times O\left(k\log k+\frac{k^{2}}{\epsilon}\right)$ matrix requiring $\tilde{O}\left(\frac{nk^{\omega}}{\epsilon}\right)$ time. ∎

6.2 Sublinear Time Ridge Regression

We now demonstrate how Theorem 25 can be leveraged to give a sublinear time, relative error algorithm for approximately solving the ridge regression problem:

[TABLE]

for PSD $A$. We begin with a lemma showing that any approximation to $A$ with small spectral norm error can be used to approximately solve (19) up to relative error.

Lemma 26** (Ridge Regression via Spectral Norm Low-Rank Approximation).**

For any PSD $A\in\mathbb{R}^{n\times n}$, $y\in\mathbb{R}^{n}$, regularization parameter $\lambda\geq 0$ and $B$ with $\|A-B\|_{2}^{2}\leq\epsilon^{2}\lambda$, let $\tilde{x}\in\mathbb{R}^{n}$ be any vector satisfying:

[TABLE]

we have:

[TABLE]

Proof.

For any $x\in\mathbb{R}^{n}$ we have:

[TABLE]

Now $\|(B-A)x\|_{2}^{2}\leq\epsilon^{2}\lambda\|x\|_{2}^{2}$ and further

[TABLE]

So for any $x$, $\|Bx-y\|_{2}^{2}+\lambda\|x\|_{2}^{2}\in(1\pm 2\epsilon)\left(\|Ax-y\|_{2}^{2}+\lambda\|x\|_{2}^{2}\right)$ which gives the lemma since any nearly optimal $\tilde{x}$ for the ridge regression problem on $B$ will also be nearly optimal for $A$. ∎

Combining Lemma 26 with Theorem 25 gives the following:

Theorem 27** (Sublinear Time Ridge Regression).**

Given any PSD $A\in\mathbb{R}^{n\times n}$, regularization parameter $\lambda\geq 0$, $y\in\mathbb{R}^{n}$, and upper bound $\tilde{s}_{\lambda}$ on the statistical dimension $s_{\lambda}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\operatorname{tr}((A^{2}+\lambda I)^{-1}A^{2})$, there is an algorithm accessing $\tilde{O}\left(\frac{n\tilde{s}_{\lambda}^{2}}{\epsilon^{4}}\right)$ entries of $A$ and running in $\tilde{O}\left(\frac{n\tilde{s}_{\lambda}^{\omega}}{\epsilon^{2\omega}}\right)$ time, which outputs $\tilde{x}$ satisfying:

[TABLE]

When $\tilde{s}_{\lambda}\ll n$ as is often the case, the above significantly improves upon state-of-the-art input sparsity time runtimes for general matrices [ACW16].

Proof.

Let $k=\frac{c\tilde{s}_{\lambda}}{\epsilon^{2}}$ for sufficiently large constant $c$. We have:

[TABLE]

So $|\{i:\lambda_{i}^{2}(A)\geq\epsilon^{2}\lambda\}|\leq\frac{2s_{\lambda}}{\epsilon^{2}}$ and for large enough $c$ and $k=\frac{c\tilde{s}_{\lambda}}{\epsilon^{2}}\geq\frac{cs_{\lambda}}{\epsilon^{2}}$ we have $\|A-A_{k}\|_{2}^{2}\leq\frac{\epsilon^{2}\lambda}{2}$ and can run Algorithm 2 with error parameter $\epsilon=\Theta(1)$ to find $M,N\in\mathbb{R}^{n\times k}$ with $\|A-MN^{T}\|_{2}^{2}\leq\epsilon^{2}\lambda$. We can then apply Lemma 26 – solving $\tilde{x}=\min_{x\in\mathbb{R}^{n}}\|MN^{T}x-y\|_{2}^{2}+\lambda\|x\|_{2}^{2}$ exactly using an SVD in $O(nk^{\omega-1})$ time. $\tilde{x}$ will be a $(1+O(\epsilon))$ approximate solution for $A$, which, after adjusting constants on $\epsilon$, gives the lemma. The runtime follows from Theorem 25 with $k=c\tilde{s}_{\lambda}/\epsilon^{2}$ and $\epsilon^{\prime}=\Theta(1)$. The $O(nk^{\omega-1})$ regression cost is dominated by the cost of computing the low-rank approximation. ∎

Note that Theorem 27 ensures that if the $k\geq\frac{cs_{\lambda}}{\epsilon^{2}}$ for some constant $c$, $\tilde{x}$ is a good approximation to the ridge regression problem. Setting $k$ properly requires some knowledge of an upper bound $\tilde{s}_{\lambda}$ on $s_{\lambda}$. A constant factor approximation to $s_{\lambda}$ can be computed in $\tilde{O}(n^{3/2}\cdot\mathop{\mathrm{missing}}{poly}(s_{\lambda}))$ time using for example a column PCP as given by Lemma 3 and binary searching for an appropriate $k$ value.

An interesting open question is if $s_{\lambda}$ be be approximated more quickly – specifically with linear dependence on $n$. This question is closely related to if it is possible to estimate the cost $\|A-A_{k}\|_{F}^{2}$ in $\tilde{O}(n\cdot\mathop{\mathrm{missing}}{poly}(k))$, which surprisingly is also open.

Acknowledgements

The authors thank IBM Almaden where part of this work was done. David Woodruff also thanks the Simons Institute program on Machine Learning and the XDATA program of DARPA for support.

Appendix A Low-Rank Approximation of $A$ via Approximation of $A^{1/2}$

We first observe that a low-rank approximation for $A^{1/2}$ does not imply a good low-rank approximation for $A$. Intuitively, if $A$ has a large top singular value, the low-rank approximation for $A$ must capture the corresponding singular direction with significantly more accuracy than a good low-rank approximation for $A^{1/2}$, in which the singular value is relatively much smaller.

Theorem 28**.**

For any $k$, $\epsilon$ there exists a PSD matrix $A$ and a rank $k$ matrix $B$ such that $\|A^{1/2}-B\|_{F}^{2}\leq(1+\epsilon)\|A^{1/2}-A_{k}^{1/2}\|_{F}^{2}$ but for every matrix $C$ in the rowspan of $B$,

[TABLE]

Notably, if we set $B=A^{1/2}P$ for some rank $k$ orthogonal projection $P$, $AP$ can be an arbitrarily bad low-rank approximation of $A$.

Proof.

Let $A\in\mathbb{R}^{n\times n}$ be a diagonal matrix with $A_{i,i}=\alpha^{2}$ for $i=1,...,k$, $A_{k+1,k+1}=0$ and all other diagonal entries equal to $\beta^{2}$, where $\alpha>\beta>0$. Let $B$ be a rank $k$ matrix which has its last $n-k$ rows all zero. For $i=1,...,k$, let $B_{i,i}=A^{1/2}_{i,i}$ and $B_{1,k+2}=\sqrt{\epsilon(n-k-1)}\cdot\beta$. We have: $\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}=(n-k-1)\beta^{2}$ and $\|A^{1/2}-B\|_{F}^{2}=(1+\epsilon)(n-k-1)\beta^{2}$. Note that the first row $b_{1}$ aligns somewhat well with $a_{1}$, but as we will see, not well enough to give a good low-rank approximation for $A$ itself.

Let $C$ be the projection of $A$ onto the rowspan of $B$, which gives the optimal low-rank approximation to $A$ within this span. For $i=2,...,k$, $c_{i}=a_{i}$, since $A$ and $B$ match exactly on these rows up to a scaling. For $i>k$, $c_{i}=\vec{0}$. Finally, $c_{1}$ = $\frac{b_{1}}{\|b_{1}\|_{2}^{2}}\cdot\langle b_{1},a_{1}\rangle=b_{1}\cdot\left(\frac{\alpha^{3}}{\alpha^{2}+\epsilon(n-k-1)\beta^{2}}\right)$. Overall:

[TABLE]

By setting $\alpha\gg\beta$ we can make this approximation arbitrarily bad. Note that $\alpha^{2}/\beta^{2}=\lambda_{1}(A)/\lambda_{k+1}(A)$. This ratio will be large whenever $A$ is well approximated by a low-rank matrix. ∎

Despite the above example, we can show that for a projection $P$, if $A^{1/2}P$ is a very near optimal low-rank approximation of $A^{1/2}$ then $A^{1/2}PA^{1/2}$ is a relative error low-rank approximation of $A$:

Theorem 29**.**

Let $P\in\mathbb{R}^{n\times n}$ be an orthogonal projection matrix such that $\|A^{1/2}-A^{1/2}P\|_{F}^{2}\leq(1+\epsilon/\sqrt{n})\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}$. Then:

[TABLE]

Proof.

We can rewrite using the fact that $(I-P)=(I-P)^{2}$ since it is a projection:

[TABLE]

Let $\delta_{i}=\sigma_{i}(A^{1/2}(I-P))$ denote the $i^{th}$ singular value of $A^{1/2}(I-P)$. Let $\lambda_{i}$ be the $i^{th}$ eigenvalue of $A$. By the assumption that $P$ gives a near optimal low-rank approximation of $A^{1/2}$:

[TABLE]

Additionally, by Weyl’s monotonicity theorem (see e.g. Theorem 3.2 in [Gu14] and proof of Lemma 15 in [MM15]), for all $i$, $\delta_{i}\geq\lambda^{1/2}_{i+k}$. We thus have:

[TABLE]

If $\lambda_{k+1}\geq 1/\sqrt{n}\cdot\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}$ then $\left(\lambda_{k+1}+\epsilon/\sqrt{n}\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}\right)^{2}\leq(1+\epsilon)^{2}\lambda_{k+1}^{2}\leq(1+3\epsilon)\lambda_{k+1}^{2}$ and hence:

[TABLE]

Alternatively if $\lambda_{k+1}\leq 1/\sqrt{n}\cdot\|A^{1/2}-A^{1/2}_{k}\|_{F}^{2}$ then

[TABLE]

which also gives the theorem. ∎

A.1 PSD Low-Rank Approximation in $n^{1.69}\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ Time

We now combine Theorem 29 with the ridge leverage score sampling algorithm of [MM16] to give a sublinear time algorithm for low-rank approximation of $A$ reading $n^{3/2}\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ entries of the matrix and running in $n^{1.69}\cdot\mathop{\mathrm{missing}}{poly}(k/\epsilon)$ time. We note that this approach could also be used with adaptive sampling [DV06] or volume sampling techniques [AGR16], as outlined in the introduction.

Theorem 30**.**

There is an algorithm based off ridge leverage score sampling which, given PSD $A\in\mathbb{R}^{n\times n}$ with high probability outputs $M\in\mathbb{R}^{n\times k}$ with $\|A-MM^{T}\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$. The algorithm reads $\tilde{O}(n^{3/2}k/\epsilon)$ entries of $A$ and runs in $\tilde{O}\left(n^{(\omega+1)/2}\cdot(k/\epsilon)^{\omega-1}\right)$ time where $\omega<2.38$ is the exponent of matrix multiplication.

This follows from Lemma 4, adapted from Theorem 8 of [MM16], which shows that it is possible to estimate the ridge leverage scores of $A^{1/2}$ with $\tilde{O}(nk)$ accesses to $A$ and $O(nk^{\omega-1})$ time. We can use these scores to sample a set of rows from $A^{1/2}$ whose span contains a near optimal low-rank approximation. Specifically we have:

Lemma 31** (Theorem 7 of [CMM17]).**

For any $B\in\mathbb{R}^{n\times n}$, for $i\in\{1,\ldots,n\}$, let $\tilde{\tau}_{i}^{k}\geq\tau_{i}^{k}(B)$ be an overestimate for the $i^{th}$ rank-$k$ ridge leverage score of $B$. Let $p_{i}=\frac{\tilde{\tau}^{k}_{i}}{\sum_{i}\tilde{\tau}^{k}_{i}}$ and $t=c\left(\log k+\frac{\log(1/\delta)}{\epsilon}\right)\sum_{i}^{k}\tilde{\tau}^{k}_{i}$ for $\epsilon<1$ and some sufficiently large constant $c$. Construct $R$ by sampling $t$ rows of $B$, each set to row $b_{i}$ with probability $p_{i}$. With probability $1-\delta$, letting $P_{R}$ be the projection onto the rows of $R$,

[TABLE]

Note that $(BP_{R})_{k}$ can be written as a row projection of $B$ – onto the top $k$ singular vectors of $BP_{R}$. So, if we compute for each $i$, $\tilde{\tau}_{i}^{k}\geq\tau_{i}^{k}(A^{1/2})$ using Lemma 4, set $\epsilon^{\prime}=\epsilon/3\sqrt{n}$, and let $S$ be a sampling matrix selecting $\tilde{O}(\sum\tilde{\tau}^{k}_{i}/\epsilon^{\prime})=\tilde{O}(k\sqrt{n}/\epsilon)$ rows of $A^{1/2}$, then by Theorem 29, letting $P$ be the projection onto the rows of $SA^{1/2}$, we have $\|A-(A^{1/2}P)_{k}(A^{1/2}P)_{k}^{T}\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$.

Finally, $(A^{1/2}P)_{k}(A^{1/2}P)_{k}^{T}=(A^{1/2}PA^{1/2})_{k}=(AS(S^{T}AS)^{+}S^{T}A)_{k}$. We can compute a factorization of this matrix by computing $(S^{T}AS)^{+/2}$, then computing $AS(S^{T}AS)^{+/2}$ and taking the SVD of this matrix. Since $S$ has $\tilde{O}(k\sqrt{n}/\epsilon)$ columns, using fast matrix multiplication this requires time $\tilde{O}(n\cdot(k\sqrt{n}/\epsilon)^{\omega-1})=\tilde{O}(n^{(\omega+1)/2}\cdot(k/\epsilon)^{\omega-1})$ and $\tilde{O}(n^{3/2}k/\epsilon)$ accesses to $A$ (to read the entries of $AS$), giving Theorem 30.

Appendix B Additional Proofs for Main Algorithm

Lemma 2 (Sum of Ridge Leverage Scores).

For any $A\in\mathbb{R}^{n\times d}$, $\sum_{i=1}^{d}\tau_{i}(A)\leq 2k$.

Proof.

We rewrite Definition 1 using $A$’s singular value decomposition $A=U\Sigma V^{T}$.

[TABLE]

where $\bar{\Sigma}_{i,i}=\frac{1}{\sigma^{2}_{i}(A)+\frac{\|A-A_{k}\|_{F}^{2}}{k}}$. We then have:

[TABLE]

$(\Sigma^{2}\bar{\Sigma})_{i,i}=\frac{\sigma_{i}^{2}(A)}{\sigma^{2}_{i}(A)+\frac{\|A-A_{k}\|_{F}^{2}}{k}}$. For $i\leq k$ we simply upper bound this by $1$. So:

[TABLE]

∎

We first prove our row sampling PCP result for spectral norm error. We follow with the closely related proof Lemma 7 which gives Frobenius norm error.

Lemma 8 (Spectral Norm Row PCP).

For any PSD $A\in\mathbb{R}^{n\times n}$, and $\epsilon<1$ let $k^{\prime}=\lceil ck/\epsilon^{2}\rceil$ and $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})\geq\tau_{i}^{k^{\prime}}(A^{1/2})$ be an overestimate for the $i^{th}$ rank-$k^{\prime}$ ridge leverage score of $A^{1/2}$. Let $\tilde{\ell}_{i}=4\epsilon\sqrt{\frac{n}{k}}\tau_{i}^{k^{\prime}}(A^{1/2})$, $p_{i}=\frac{\tilde{\ell}_{i}}{\sum_{i}\tilde{\ell}_{i}}$, and $t=\frac{c^{\prime}\log n}{\epsilon^{2}}\cdot\sum_{i}\tilde{\ell}_{i}$. Construct weighted sampling matrices $S_{1},S_{2}\in\mathbb{R}^{n\times t}$, where the $j^{th}$ column is set to $\frac{1}{\sqrt{tp_{i}}}e_{i}$ with probability $p_{i}$.

For sufficiently large constants $c,c^{\prime}$, with high probability, letting $\tilde{A}=S_{2}^{T}AS_{1}$, for any orthogonal projection $P\in\mathbb{R}^{t\times t}$:

[TABLE]

We refer to $\tilde{A}$ as an $(\epsilon,k)$-spectral PCP of $AS_{1}$.

Note that if $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})$ is a constant factor approximation to $\tau_{i}^{k^{\prime}}(A^{1/2})$, $t=O\left(\frac{\sqrt{nk}\log n}{\epsilon^{3}}\right)$. We use ‘with high probability’ to mean with probability $\geq 1-1/n^{d}$ for some large constant $d$.

Proof.

For conciseness write $C\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}AS_{1}$ and write the eigendecomposition $A=U\Lambda U^{T}$ with $\lambda_{i}=\Lambda_{i,i}$. Applying the triangle inequality we have:

[TABLE]

Thus to show the Lemma it suffices to show:

[TABLE]

Let $m$ be the largest index with $\lambda_{m}^{2}\geq\frac{\epsilon^{2}}{k}\|A-A_{k}\|_{F}^{2}$. Let $U_{H}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}U_{m}$ contain the top $m$ ‘head’ eigenvectors of $A$ and let $U_{T}$ contain the remaining ‘tail’ eigenvectors. Let $C_{H}=U_{H}U_{H}^{T}C$ and $C_{T}=U_{T}U_{T}^{T}C$. $C_{H}+C_{T}=C$ and so by the triangle inequality:

[TABLE]

We bound each of the terms in the above sum separately. Specifically we show:

• •

Head Term: $\|(I-P)(C_{H}^{T}S_{2}S_{2}^{T}C_{H}-C_{H}^{T}C_{H})(I-P)\|_{2}\leq\epsilon\|C(I-P)\|_{2}^{2}$

• •

Tail Term: $\|(I-P)(C_{T}^{T}S_{2}S_{2}^{T}C_{T}-C_{T}^{T}C_{T})(I-P)\|_{2}\leq\frac{13\epsilon^{2}}{k}\|A-A_{k}\|_{F}^{2}$.

• •

Cross Term: $\|(I-P)C_{H}^{T}S_{2}S_{2}^{T}C_{T}(I-P)\|_{2}\leq\frac{10\epsilon}{k}\|A-A_{k}\|_{F}^{2}+4\epsilon\|C(I-P)\|_{F}^{2}$.

Combining these three bounds, after adjusting constant factors on $\epsilon$ by making the constants $c$ and $c^{\prime}$ in the rank parameter $k^{\prime}$ and sample size $t$ large enough, gives (20) and thus the lemma. For the remainder of the proof we thus fix $c=1$ so $k^{\prime}=\lceil k/\epsilon^{2}\rceil$.

Head Term:

We first show that the ridge scores of $A^{1/2}$ upper bound the standard leverage scores of $U_{m}$.

Lemma 32**.**

For any $p$ with $\lambda^{2}_{p}(A)\geq\frac{1}{k}\|A-A_{k}\|_{F}^{2}$ we have:

[TABLE]

where $\|(U_{p})_{i}\|_{2}^{2}$ is the $i^{th}$ row norm of $U_{p}$, whose columns are the top $p$ eigenvectors of $A$.

Proof.

If $U_{p}$ contains no eigenvectors, this is true vacuously as all row norms are [math]. Otherwise

[TABLE]

where $\hat{\Lambda}_{j,j}=\frac{\lambda_{j}^{2}}{\lambda_{j}^{2}+\frac{\|A-A_{k}\|_{F}^{2}}{k}}.$ We can then write:

[TABLE]

This gives the lemma combined with the fact that $\tau_{i}^{k}(A)\leq 2\sqrt{\frac{n}{k}}\tau_{i}^{k}(A^{1/2})$ by Lemma 5. ∎

Applying Lemma 32 with $k^{\prime}=\lceil k/\epsilon^{2}\rceil$, since $\lambda_{m}^{2}(A)\geq\frac{1}{k^{\prime}}\|A-A_{k^{\prime}}\|_{F}^{2}$, we have:

[TABLE]

Further, $t=\frac{c^{\prime}\log n}{\epsilon^{2}}\cdot\sum_{i}\tilde{\ell}_{i}$ so if we set $c^{\prime}$ large enough, we have by a standard matrix Chernoff bound (see Lemma 33) since $U_{H}$ is an orthogonal span for the columns of $C_{H}$ and hence and its row norms are the leverage scores of $C_{H}$, with high probability:

[TABLE]

This in turn gives:

[TABLE]

Tail Term:

We can loosely bound via the triangle inequality and the fact that $\|I-P\|_{2}\leq 1$:

[TABLE]

Since $k^{\prime}=\lceil k/\epsilon^{2}\rceil$:

[TABLE]

Thus for $S=S_{1}$ or $S=S_{2}$, for sufficiently large $c^{\prime}$ in our sample size $t=\frac{c^{\prime}\log n}{\epsilon^{2}}$, by a matrix Chernoff bound (Corollary 34), with high probability

[TABLE]

Which in turn implies

[TABLE]

and so, applying the AMGM inequality,

[TABLE]

Similarly we have:

[TABLE]

So overall, plugging back into (24) gives:

[TABLE]

Cross Term:

By submultiplicativity of the spectral norm:

[TABLE]

As shown above for our tail bound, $\|C_{T}^{T}S_{2}\|_{2}\leq\sqrt{\frac{10\epsilon^{2}}{k}}\|A-A_{k}\|_{F}$. Further, as shown in (22), $S_{2}$ gives a subspace embedding for $C_{H}$ so we have:

[TABLE]

where the last bound follows form the AMGM inequality.

Plugging our head (23), tail (25) and cross term (26) bounds in (B), and adjusting constants on $\epsilon$ by making $c$ and $c^{\prime}$ sufficiently large gives the lemma.

∎

Lemma 7 (Frobenius Norm Row PCP).

For any PSD $A\in\mathbb{R}^{n\times n}$ and $\epsilon\leq 1$ let $k^{\prime}=\lceil ck/\epsilon\rceil$ and let $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})\geq\tau_{i}^{k^{\prime}}(A^{1/2})$ be an overestimate for the $i^{th}$ rank-$k^{\prime}$ ridge leverage score of $A^{1/2}$. Let $\tilde{\ell}_{i}=\sqrt{\frac{16n\epsilon}{k}}\cdot\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})$, $p_{i}=\frac{\tilde{\ell}_{i}}{\sum_{i}\tilde{\ell}_{i}}$, and $t=\frac{c^{\prime}\log n}{\epsilon^{2}}\sum_{i}\tilde{\ell}_{i}$. Construct weighted sampling matrices $S_{1},S_{2}\in\mathbb{R}^{n\times t}$ each whose $j^{th}$ column is set to $\frac{1}{\sqrt{tp_{i}}}e_{i}$ with probability $p_{i}$.

For sufficiently large constants $c,c^{\prime}$, with probability $\frac{99}{100}$, letting $\tilde{A}=S_{2}^{T}AS_{1}$, for any rank-$k$ orthogonal projection $P\in\mathbb{R}^{t\times t}$:

[TABLE]

for some fixed $\Delta$ (independent of $P$) with $|\Delta|\leq 600\|A-A_{k}\|_{F}^{2}$.

Note that if $\tilde{\tau}_{i}^{k^{\prime}}(A^{1/2})$ is a constant factor approximation to $\tau_{i}^{k^{\prime}}(A^{1/2})$, $t=O\left(\frac{\sqrt{nk}\log n}{\epsilon^{2.5}}\right)$.

Proof.

The proof is similar to that of Lemma 8. Denote $C\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}AS_{1}$ and write the eigendecomposition $A=U\Lambda U^{T}$ with $\lambda_{i}=\Lambda_{i,i}$. Let $m$ be the largest index with $\lambda_{m}^{2}\geq\frac{\epsilon}{k}\|A-A_{k^{\prime}}\|_{F}^{2}.$ Let $U_{H}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}U_{m}$ contain the top $m$ ‘head’ eigenvectors of $A$ and let $U_{T}$ contain the remaining ‘tail’ eigenvectors. Let $C_{H}=U_{H}U_{H}^{T}C$ and $C_{T}=U_{T}U_{T}^{T}C$ and note that $C_{H}+C_{T}=C$.

By the Pythagorean theorem we have $\|C(I-P)\|_{F}^{2}=\|C_{H}(I-P)\|_{F}^{2}+\|C_{T}(I-P)\|_{F}^{2}$. Expanding using the identity $\|M\|_{F}^{2}=\operatorname{tr}(M^{T}M)$,

[TABLE]

We bound each of the terms in the above sum separately. Specifically we show:

• •

Head Term: $\|S_{2}^{T}C_{H}(I-P)\|_{F}^{2}\in(1\pm\epsilon)\|C_{H}(I-P)\|_{F}^{2}$

• •

Tail Term: $\|S_{2}^{T}C_{T}(I-P)\|_{F}^{2}+\Delta\in\|C_{T}(I-P)\|_{F}^{2}\pm\epsilon\|A-A_{k}\|_{F}^{2}$ where $|\Delta|\leq 600\|A-A_{k}\|_{F}^{2}$.

• •

Cross Term: $|\operatorname{tr}\left((I-P)C_{H}^{T}S_{2}S_{2}^{T}C_{T}(I-P)\right)|\leq\epsilon\|C(I-P)\|_{F}^{2}.$

It is not hard to see that combining these bounds gives the Lemma. $\|A-A_{k}\|_{F}^{2}\leq(1+3\epsilon)\|C(I-P)\|_{F}^{2}$ since $C$ is an $(\epsilon,k)$-column PCP of $A$ by Lemmas 3, 5, and the fact that rank-$k^{\prime}$ ridge scores upper bound the rank-$k$ scores (as long as $c>1$). Additionally, $\|A-A_{k}\|_{F}^{2}\leq\|A(I-P)\|_{F}^{2}$ for any rank-$k$ $P$. So plugging into (27) we have:

[TABLE]

This gives the lemma by adjusting constants on $\epsilon$ by making $c$ and $c^{\prime}$ large enough. For the remainder of the proof we thus fix $c=1$ and so $k^{\prime}=\lceil k/\epsilon\rceil$.

Head Term:

By Lemma 32 applied to $k^{\prime}=\lceil k/\epsilon\rceil$, since $\lambda_{m}^{2}\geq\frac{1}{k^{\prime}}\|A-A_{k^{\prime}}\|_{F}^{2}$:

[TABLE]

Further, $t=\frac{c^{\prime}\log n}{\epsilon^{2}}\cdot\sum_{i}\tilde{\ell}_{i}$ so if we set $c^{\prime}$ large enough, by a standard matrix Chernoff bound (see Lemma 33) since $U_{H}$ is an orthogonal span for the columns of $C_{H}$ and hence and its row norms are the leverage scores of $C_{H}$, with high probability:

[TABLE]

This gives the bound $\|S_{2}^{T}C_{H}(I-P)\|_{F}^{2}\in(1\pm\epsilon)\|C_{H}(I-P)\|_{F}^{2}$.

Tail Term: We want to show:

[TABLE]

where $|\Delta|\leq 600\|A-A_{k}\|_{F}^{2}$.

We again split using Pythagorean theorem, $\|S_{2}^{T}C_{T}(I-P)\|_{F}^{2}=\|S_{2}^{T}C_{T}\|_{F}^{2}-\|S_{2}^{T}C_{T}P\|_{F}^{2}$. We set $\Delta=\|C_{T}\|_{F}^{2}-\|S_{2}^{T}C_{T}\|_{F}^{2}$. We have $\operatorname*{\mathbb{E}}\|S_{2}C_{T}\|_{F}^{2}=\|C_{T}\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$, where the last bound follows since $C$ is an $(\epsilon,k)$-column PCP for $A$. Thus by a Markov bound, with probability $299/300$, $|\Delta|\leq(1+\epsilon)300\|A-A_{k}\|_{F}^{2}\leq 600\|A-A_{k}\|_{F}^{2}$.

Additionally, we have $\|C_{T}\|_{2}^{2}\leq\frac{10\epsilon}{k}\|A-A_{k}\|_{F}^{2}$ and $\|S_{2}^{T}C_{T}\|_{2}^{2}\leq\frac{10\epsilon}{k}\|A-A_{k}\|_{F}^{2}$ with high probability by an identical argument to that used for Lemma 8. This gives:

[TABLE]

which gives the main bound (28) after adjusting constants on $\epsilon$.

Cross Term: We want to show:

[TABLE]

We can write:

[TABLE]

where in the last step, inserting $(C^{T}C)^{+}(C^{T}C)$, which is the projection onto the row span of $C$ has no effect as the rows of $C_{H}=U_{H}U_{H}^{T}C$ already lie in this span. $\langle M,N\rangle=\operatorname{tr}(M(C^{T}C)^{+}N^{T})$ is a semi-inner product since $C^{T}C$ is positive semidefinite, so by Cauchy-Schwarz:

[TABLE]

Using the singular value decomposition $C=XSY^{T}$:

[TABLE]

As argued, by Lemma 32 we have $\tilde{\ell}_{i}\geq\|(U_{H})_{i}\|_{2}^{2}$, and so by a standard approximate matrix multiplication result [DKM06], with probability $299/300$ if we set the constant $c^{\prime}$ in our sample size $>c^{\prime\prime}$ for some fixed $c^{\prime\prime}$:

[TABLE]

where the last step follows from the fact that $\|C_{T}\|_{F}^{2}=O(\|A-A_{k}\|_{F}^{2})$ since $C$ is an $(\epsilon,k)$-column PCP for $A$. The final bound then follows from combining with (30) with the fact that $\|A-A_{k}\|_{F}^{2}\leq(1+\epsilon)\|C(I-P)\|_{F}^{2}$ and adjusting constants on $\epsilon$ by increasing the constant $c^{\prime}$ in the sample size $t$ and $c$ in the rank parameter $k^{\prime}=\lceil ck/\epsilon\rceil$.

The full lemma follows simply from noting that we can union bound over our failure probability for each term so all hold simultaneously with probability $99/100$. ∎

Lemma 33** (Leverage Score Sampling Matrix Chernoff).**

For any $A\in\mathbb{R}^{n\times d}$ for $i\in\{1,\ldots,d\}$, let $\ell_{i}(A)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}a_{i}^{T}(AA^{T})^{+}a_{i}$ be the $i^{th}$ column leverage score of $A$ and let $\tilde{\ell}_{i}\geq\ell_{i}(A)$ be an overestimate for this score. Let $p_{i}=\frac{\tilde{\ell}_{i}}{\sum_{i}\tilde{\ell}_{i}}$ and $t=\frac{c\log(d/\delta)}{\epsilon^{2}}\sum_{i}\tilde{\ell}_{i}$ for sufficiently large $c$. Construct $C$ by sampling $t$ columns of $A$, each set to $\frac{1}{\sqrt{tp_{i}}}a_{i}$ with probability $p_{i}$. With probability $1-\delta$:

[TABLE]

Proof.

Write the singular value decomposition $A=U\Sigma V^{T}$. Note that:

[TABLE]

Let $Y=\Sigma^{-1}U^{T}\left(CC^{T}-AA^{T}\right)U\Sigma^{-1}$. We can write:

[TABLE]

For each $j\in 1,\ldots,t$, $X_{j}$ is given by:

[TABLE]

$\operatorname*{\mathbb{E}}Y=0$ since $\operatorname*{\mathbb{E}}X_{j}=\sum_{i=1}^{d}p_{i}\left[\frac{1}{p_{i}}a_{i}a_{i}^{T}-AA^{T}\right]=0$. Furthermore, $CC^{T}=U\Sigma Y\Sigma U+AA^{T}$. Showing $\|Y\|_{2}\leq\epsilon$ gives $-\epsilon I\preceq Y\preceq\epsilon I$, which gives the lemma. We prove this bound using a matrix Bernstein inequality from [Tro15]. This inequality requires upper bounds on the spectral norm of each $X_{j}$ and on variance of $Y$. We first note that for any $i$, $\frac{1}{\ell_{i}(A)}a_{i}a_{i}^{T}\preceq AA^{T}.$ This follows from writing any $x$ in the column span of $A$ as $(AA^{T})^{+/2}y$ and then noting:

[TABLE]

since $(AA^{T})^{+/2}a_{i}a_{i}^{T}(AA^{T})^{+/2}$ is rank-$1$ and so has maximum eigenvalue $\operatorname{tr}\left((AA^{T})^{+/2}a_{i}a_{i}^{T}(AA^{T})^{+/2}\right)=\ell_{i}(A)$ by the cyclic property of trace. Further $x^{T}AA^{T}x=y^{T}(AA^{T})^{+/2}AA^{T}(AA^{T})^{+/2}y=\|y\|_{2}^{2}$, giving us the bound. We then have:

[TABLE]

And so $\frac{1}{tp_{i}}\Sigma^{-1}U^{T}a_{i}a_{i}^{T}U\Sigma^{-1}\preceq\frac{\epsilon^{2}}{c\log(d/\delta)\tilde{\ell}_{i}}\Sigma^{-1}U^{T}a_{i}a_{i}^{T}U\Sigma^{-1}\preceq\frac{\epsilon^{2}}{c\log(d/\delta)}I$. Additionally,

[TABLE]

as long as $\tilde{\ell}_{i}\leq 1$ which we may as well enforce since $\ell_{i}(A)\leq 1$. Overall this gives $\|X_{j}\|_{2}\leq\frac{\epsilon^{2}}{c\log(d/\delta)}$. Next we bound the variance of $Y$.

[TABLE]

By Theorem 7.3.1 of [Tro15], for $\epsilon<1$,

[TABLE]

Which gives $\Pr\left[\|Y\|_{2}\geq\epsilon\right]\leq 4de^{-\frac{c\log(d/\delta)}{4}}\leq\delta$ if we choose $c$ large enough. ∎

The above Lemma also gives an easy Corollary for the approximation obtained when sampling with ridge leverage scores:

Corollary 34** (Ridge Leverage Scores Sampling Matrix Chernoff).**

For any $A\in\mathbb{R}^{n\times d}$ for $i\in\{1,\ldots,d\}$, let $\tau^{\lambda}_{i}(A)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}a_{i}^{T}(AA^{T}+\lambda I)^{+}a_{i}$ be the $i^{th}$ $\lambda$-ridge leverage score of $A$ and let $\tilde{\tau}_{i}^{\lambda}\geq\tau_{i}^{\lambda}(A)$ be an overestimate for this score. Let $p_{i}=\frac{\tilde{\tau}_{i}^{\lambda}}{\sum_{i}\tilde{\tau}_{i}^{\lambda}}$ and $t=\frac{c\log(d/\delta)}{\epsilon^{2}}\sum_{i}\tilde{\tau}_{i}^{\lambda}$ for sufficiently large $c$. Construct $C$ by sampling $t$ columns of $A$, each set to $\frac{1}{\sqrt{tp_{i}}}a_{i}$ with probability $p_{i}$. With probability $1-\delta$:

[TABLE]

Proof.

We can instantiate Lemma 33 with $[A,\sqrt{\lambda I}]$ setting $\tilde{\ell}_{i}=\tilde{\tau}_{i}^{\lambda}$. We simply fix the columns of the identity to appear in our sample. This only decreases variance, all calculations go through, and with probability $\delta$:

[TABLE]

which gives the desired bound if we subtract $\lambda I$ from all sides. ∎

Appendix C Outputting a PSD Low-Rank Approximation

In this section, we prove Theorem 12, showing to to efficiently output a low-rank approximation to $A$ under the restriction that the approximation is also PSD. We start with the following lemma, which shows that as long as we have a good low rank subspace for approximating $A$ in the spectral norm (computable using Theorem 25), we can quickly find a near optimal PSD low-rank approximation:

Lemma 35**.**

Given PSD $A\in\mathbb{R}^{n\times n}$ and orthonormal basis $Z\in R^{n\times m}$ with $\|A-AZZ^{T}\|_{2}^{2}\leq\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}$, there is an algorithm which accesses $O\left(\frac{nm\log m}{\epsilon^{2}}\right)$ entries of $A$, runs in $\tilde{O}\left(\frac{nm^{\omega-1}}{\epsilon^{2(\omega-1)}}\right)$ time, and with probability $99/100$ outputs $M\in\mathbb{R}^{n\times k}$ satisfying $\|A-MM^{T}\|_{F}^{2}\leq(1+\epsilon)\|A-A_{k}\|_{F}^{2}$.

Proof.

By Lemma 10 of [CW17] for any basis $Z\in\mathbb{R}^{n\times m}$ with $\|A-AZZ^{T}\|_{2}^{2}\leq\frac{\epsilon}{k}\|A-A_{k}\|_{F}^{2}$,

[TABLE]

As in Algorithm 1 we can find a near optimal $X$ by further sampling $Z$ using its leverage scores. If we sample $t_{1}=O\left(\frac{m\log m}{\epsilon^{2}}\right)$ rows of $Z$ by their leverage scores (their norms since $Z$ is orthonormal) to form $S_{1}\in R^{n\times t_{1}}$, by Theorem 39 of [CW13], we will have an affine embedding of $Z$. Specifically, letting $B^{*}=\operatorname*{arg\,min}_{B}\|A-ZB\|_{F}^{2}$ and $E^{*}=A-ZB^{*}$, for any $B$ we have:

[TABLE]

Note that this is similar to the embedding property used in the proof of Theorem 9 to show that $W$ computed in Step 6 of Algorithm 1 gave a near optimal low-rank approximation to $AS_{1}$.

By a Markov bound, since $\operatorname*{\mathbb{E}}\|S_{1}^{T}E^{*}\|_{F}^{2}=\|E^{*}\|_{F}^{2}$, with probability $99/100$, $\left|\|E^{*}\|_{F}^{2}-\|S_{1}^{T}E^{*}\|_{F}^{2}\right|\leq 100\|E^{*}\|_{F}^{2}=O(1)\|A-A_{k}\|_{F}^{2}$. This guarantees that a $(1+\epsilon)$ approximation to the sketched problem gives a $(1+O(\epsilon))$ approximation to the original. That is, for any PSD $\tilde{X}$ with $\operatorname{rank}(\tilde{X})=k$, and

[TABLE]

we have $\|A-Z\tilde{X}Z^{T}\|_{F}^{2}\leq(1+O(\epsilon))\|A-A_{k}\|_{F}^{2}.$ Following [CW17], we write $S_{1}^{T}Z$ in its SVD $S_{1}^{T}Z=U_{z}\Sigma_{z}V_{z}^{T}$. Since $S_{1}ZXZ^{T}$ falls in the column span of $S_{1}Z$ and the row span of $Z^{T}$, we write $\delta\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\|(I-U_{z}U_{z}^{T})S_{1}^{T}A\|_{F}^{2}+\|U_{z}U_{z}^{T}S_{1}^{T}A(I-ZZ^{T})\|_{F}^{2}$ and by Pythagorean theorem have:

[TABLE]

Since $S_{1}$ is sampled via $Z$’s leverage scores, $(1-\epsilon)I\preceq S_{1}^{T}Z\preceq(1+\epsilon)I$ and so:

[TABLE]

Finally, following [CW17], letting $B=V_{z}^{T}\Sigma_{z}^{-1}U_{z}^{T}S_{1}^{T}AZ$, we have

[TABLE]

where $N_{k,+}$ has all but the top $k$ positive eigenvalues of $N$ set to [math]. We output $M=Z\tilde{X}^{1/2}$. Overall, the above algorithm requires accessing $O\left(\frac{nm\log m}{\epsilon^{2}}\right)$ entries of $A$ (the entries of $AS_{1}$) and has runtime $\tilde{O}\left(\frac{nm^{\omega-1}}{\epsilon^{2(\omega-1)}}\right)$ giving the lemma. ∎

We can obtain $Z$ with rank $m=\Theta(k/\epsilon)$ by applying Theorem 25 with rank $k^{\prime}=\Theta(k/\epsilon)$ and error parameter $\epsilon^{\prime}=\Theta(1)$. Combined with Lemma 35 this yields:

Theorem 12 (Sublinear Time Low-Rank Approximation – PSD Output).

There is an algorithm that given any PSD $A\in\mathbb{R}^{n\times n}$, accesses $\tilde{O}\left(\frac{nk^{2}}{\epsilon^{2}}+\frac{nk}{\epsilon^{3}}\right)$ entries of $A$, runs in $\tilde{O}\left(\frac{nk^{\omega}}{\epsilon^{\omega}}+\frac{nk^{\omega-1}}{\epsilon^{3(\omega-1)}}\right)$ time and with probability at least $9/10$ outputs $M\in\mathbb{R}^{n\times k}$ with:

[TABLE]