Input Sparsity Time Low-Rank Approximation
via Ridge Leverage Score Sampling
Michael B. Cohen

Cameron Musco

Christopher Musco∗

arXiv:1511.07263v2 [] 6 Oct 2016

Massachusetts Institute of Technology, EECS
Cambridge, MA 02139, USA
Email: {micohen,cnmusco,cpmusco}@mit.edu
October 10, 2016

Abstract
We present a new algorithm for finding a near optimal low-rank approximation of a matrix
A in O(nnz(A)) time. Our method is based on a recursive sampling scheme for computing a
representative subset of A’s columns, which is then used to find a low-rank approximation.
This approach differs substantially from prior O(nnz(A)) time algorithms, which are all
based on fast Johnson-Lindenstrauss random projections. It matches the guarantees of these
methods while offering a number of advantages.
Not only are sampling algorithms faster for sparse and structured data, but they can also
be applied in settings where random projections cannot. For example, we give new single-pass
streaming algorithms for the column subset selection and projection-cost preserving sample problems. Our method has also been used to give the fastest algorithms for provably approximating
kernel matrices [MM16].

∗

Part of this work was completed while the author interned at Yahoo Labs, NYC.

1

Introduction

Low-rank approximation is a fundamental task in statistics, machine learning, and computational
science. The goal is to find a rank k matrix that is as close as possible to an arbitrary input matrix
A ∈ Rn×d , with distance typically measured using the spectral or Frobenius norms.
Traditionally, the problem is solved using the singular value decomposition (SVD), which takes
O(nd2 ) time to compute. This high cost can be reduced using iterative algorithms like the power
method or Krylov methods, which require just O(nnz(A) · k) time per iteration, where nnz(A) is
the number of non-zero entries in A1 .
More recently, the cost of low-rank approximation has been reduced even further using sketching
methods based on Johnson-Lindenstrauss random projection [Sar06]. Remarkably, so-called “sparse
random projections” [CW13, MM13, NN13, BDN15, Coh16] give algorithms that run in time2 :
O(nnz(A)) + Õ (n · poly(k, ǫ)) .
These methods output a low-rank approximation within a (1 + ǫ) factor of optimal when error is
measured using the Frobenius norm. They are typically referred to as running in “input sparsity
time” since the O(nnz(A)) term is considered to dominate the runtime.
Input sparsity time algorithms for low-rank approximation were an important theoretical achievement and have also been influential in practice. Implementations are now available in a variety of
languages and machine learning libraries [Liu14, Oka10, IBM14, PVG+ 11, HRZ+ 09, VM15].

1.1

Our Contributions

We give an entirely different approach to obtaining input sparsity time algorithms for low-rank
approximation. Random projection methods are based on multiplying A by a sparse random
matrix Π ∈ Rd×poly(k,ǫ) to form a smaller matrix AΠ that contains enough information about A
to compute a near optimal low-rank approximation. Our techniques on the other hand are based
on sampling Õ(k/ǫ) columns from A, and computing a low-rank approximation using this sample.
Sampling itself is simple and extremely efficient. However, to obtain a good approximation
to A, columns must be sampled with non-uniform probabilities, carefully chosen to reflect their
relative importance. It is known that variations on the standard statistical leverage scores give
probabilities that are provably sufficient for low-rank approximation [Sar06, DMM08, CEM+ 15].
Unfortunately, computing any of these previously studied “low-rank leverage scores” is as difficult as low-rank approximation itself, so sampling did not yield fast algorithms3 .
We address this issue for the first time by introducing new importance sampling probabilities
which can be approximated efficiently using a simple recursive algorithm. In particular, we adapt
the so-called ridge leverage scores to low-rank matrix approximation. These scores have been used
as sampling probabilities in the context of linear regression and spectral approximation [LMP13,
KLM+ 14, AM15] but never for low-rank approximation.
By showing that ridge leverage scores display a unique monotonicity property under perturbations to A, we are able to prove that, unlike any prior low-rank leverage scores, they can be approximated using a relatively large uniform subsample of A’s columns. While too large to use directly,
the size of this subsample can be reduced recursively to give an overall fast algorithm. This approach
1

The number of iterations depends on the accuracy ǫ and/or spectral gap conditions. See [MM15] for an overview.
Õ(·) hides logarithmic factors, including a failure probability dependence.
3
ℓ2 norm sampling does yield very fast algorithms [FKV04, DKM06a, BJS15], but cannot give relative error
guarantees matching those of random projection or leverage score methods without additional assumptions on A.
2

1

resembles work on recursive methods for computing standard leverage scores, which were recently
used to give the first O(nnz(A)) time sampling algorithms for linear regression [LMP13, CLM+ 15].
Our main algorithmic result, which nearly matches the state-of-the-art in [NN13] follows:
Theorem 1. For any
 1], there exists a recursive column sampling algorithm that, in time
 θ1+θ∈ 2(0,

n
k
−1
, returns Z ∈ Rn×k satisfying:
O θ nnz(A) + Õ
ǫ4
kA − ZZT Ak2F ≤ (1 + ǫ)kA − Ak k2F .

(1)

Here Ak is the best rank k approximation to A. To prove Theorem 1, we show how to compute
2
a sampling matrix S such that AS ∈ Rn×Õ(k/ǫ ) satisfies a projection-cost preservation guarantee
(formalized in Section 2). This property ensures that it is possible to extract a near optimal lowrank approximation from the sample. It also allows AS to be used to approximately solve a broad
class of constrained low-rank approximation problems, including k-means clustering [CEM+ 15].
With a slightly smaller sample, we also prove that AS satisfies a standard (1 + ǫ) error column subset selection guarantee. Ridge leverage score sampling is the first algorithm, efficient or
otherwise, that obtains both of these important approximation goals simultaneously.

1.2

Why Sampling?

Besides the obvious goal of obtaining alternative state-of-the-art algorithms for low-rank approximation, we are interested in sampling methods for a few specific reasons:
Sampling maintains matrix sparsity and structure.
Without additional assumptions on A, our recursive sampling algorithms essentially match random
projection methods. However, they have the potential to run faster for sparse or structured data.
Random projection linearly combines all columns in A to form AΠ ∈ Rn×poly(k,ǫ), so this sketched
matrix is usually dense and unstructured. On the other hand, AS will remain sparse or structured
if A is sparse or structured, in which case it can be faster to post-process. Potential gains are
especially important when the Õ(n · poly(k, ǫ)) runtime term is not dominated by O(nnz(A)).
We note that the ability to maintain sparsity and structure motivated similar work on recursive sampling algorithms for fast linear regression [CLM+ 15]. While these techniques only match
random projection for general matrices, they have been important ingredients in designing faster algorithms for highly structured Laplacian and SDD matrices [LPS15, KLP+ 16, JK16]. We hope our
sampling methods for low-rank approximation will provide a foundation for similar contributions.
Sampling techniques lead to natural streaming algorithms.
In data analysis, sampling itself is often the primary goal. The idea is to select a subset of columns
from A whose span contains a good low-rank approximation for the matrix and hence represents
important or influential features [PZB+ 07, MD09].
Computing this subset in a streaming setting is of both theoretical and practical interest [Str14].
Unfortunately, while random projection methods adapt naturally to data streams [CW09], importance sampling is more difficult. The leverage score of one column depends on every other column,
including those that have not yet appeared in the stream. While random projections can be used
to approximate leverage scores, this approach inherently requires two passes over the data.
Fortunately, the same techniques used in our recursive algorithms apply naturally in the streaming setting. We can compute coarse approximations to the ridge leverage scores using just a small
number of columns and refine these approximations as the stream is revealed. By rejection sampling
columns as the probabilities are adjusted, we obtain the first space efficient single-pass streaming
algorithms for both column subset selection and projection-cost preserving sampling (Section 6).
2

Sampling offers additional flexibility for non-standard matrices.
In recent follow up work, the techniques in this paper are adapted to give the most efficient,
provably accurate algorithms for kernel matrix approximation [MM16]. In nearly all settings this
well-studied problem cannot be solved efficiently by random projection methods.
The goal in kernel approximation is to replace an n × n positive semidefinite kernel matrix
K with a low-rank approximation that takes less space to represent [AMS01, WS01, FS02, MD05,
RR07, BW09a, BW09b, Bac13, GM13, HI15, LJS16]. However, unlike in the standard low-rank approximation problem, K is not represented explicitly. Its entries can only be accessed by evaluating
a “kernel function” between each pair of the n points in a data set.
Sketching K using random projection requires computing the full matrix first, using Θ(n2 )
kernel evaluations. On the other hand, with recursive ridge leverage score sampling this is not
necessary – it is possible to compute entries of K ‘on the fly’, only when they are required to
compute ridge scores with respect to a subsample. [MM16] shows that this technique gives the first
provable algorithms for approximating kernel matrices that only require time linear in n. In other
words, the methods only evaluate a tiny fraction of the dot products required to build K. Notably
they do not require any coherence or regularity assumptions to achieve this runtime.
Aside from kernel approximation, we note that in [BJS15] the authors present a low-rank
approximation algorithm based on elementwise sampling that they show can be applied to the
product of two matrices without ever forming this product explicitly. This result again highlights
the flexibility of sampling-based methods for non-standard matrices. Without access to an efficiently
computable leverage score distribution for elementwise sampling, [BJS15] applies an approximation
based on ℓ2 sampling. This approximation only performs well under additional assumptions on A
and an interesting open question is to see if our techniques can be adapted to their framework in
order to eliminate such assumptions.

1.3

Techniques and Paper Layout

Sampling Bounds (Sections 2, 3): In Section 2 we review technical background and introduce
ridge leverage scores. In Section 3 we prove that sampling by ridge leverage scores gives solutions
to the projection-cost preserving sketch and column subset selection problems. These sections do
not address algorithmic considerations.
While the proofs are technical, we reduce both problems to a simple “additive-multiplicative
spectral guarantee,” which resembles the ubiquitous subspace embedding guarantee [Sar06]. This
approach greatly simplifies prior work on low-rank approximation bounds for sampling methods
[CEM+ 15] and we hope that it will prove generally useful in studying future sketching methods.
Ridge Leverage Score Monotonicity (Section 4): In Section 4 we prove a basic theorem
regarding the stability of ridge leverage scores. Specifically, we show that the ridge leverage score
of a column cannot increase if another column is added to the matrix. This fact, which do not
hold for any prior “low-rank leverage scores”, is essential in proving the correctness of our recursive
sampling procedure and streaming algorithms.
Recursive Sampling Algorithm (Section 5): In Section 5, we describe and prove the correctness of our main sampling algorithm. We show how a careful implementation of the algorithm
gives O(nnz(A)) running time for computing ridge leverage scores, and accordingly for solving the
low-rank approximation problem.
Application to Streaming (Section 6): We conclude with an application of our results to
low-rank sampling algorithms for single-pass column streams that are only possible thanks to the
stability result of Section 4.
3

2
2.1

Technical Background
Low-rank Approximation

Using the singular value decomposition (SVD), any rank r matrix A ∈ Rn×d can be factored as
A = UΣVT . U ∈ Rn×r and V ∈ Rd×r are orthonormal matrices whose columns are the left and
right singular vectors of A. Σ is a diagonal matrix containing A’s singular values σ1 ≥ σ2 ≥ . . . ≥
σr > 0 in decreasing order from top left to bottom right. When quality is measured with respect
to the Frobenius norm, the best low-rank approximation for A is given by Ak = Uk Σk VkT where
Uk ∈ Rn×k , Vk ∈ Rd×k , and Σk ∈ Rk×k contain the just the first k components of U, V, and Σ
respectively. In other words,
Ak =

arg min kA − BkF .

B:rank(B)≤k

Since U has orthonormal columns, we can rewrite Ak = Uk UTk A. That is, the best rank k
approximation can be found by projecting A onto the span of its top k singular vectors. Throughout,
we will use the shorthand A\k to denote the residual A − Ak . U\k ∈ Rn×r−k , V\k ∈ Rd×r−k , and
Σ\k ∈ Rr−k×r−k denote U, V, and Σ restricted to just their last k components.
When solving the low-rank approximation problem approximately, our goal is to find an orthonormal span Z ∈ Rn×k satisfying kA − ZZT AkF ≤ (1 + ǫ)kA − Uk UTk AkF .

2.2

Sketching Algorithms

Like many randomized linear algebra routines, our low-rank approximation algorithms are based
on “linear sketching”. Sketching algorithms use a typically randomized procedure to compress A ∈
′
Rd×n into an approximation (or “sketch”) C ∈ Rd ×n with many fewer columns (d′ ≪ d). Random
projection algorithms construct C by forming d′ random linear combinations of the columns in A.
Random sampling algorithms construct C by selecting and possibly reweighting a d′ columns in A.
After compression, a post-processing routine, which is often deterministic, is used to solve the
original linear algebra problem with just the information contained in C. In our case, the postprocessing step needs to extract an approximation to the span of A’s top left singular vectors. If
C is much smaller than A, the cost of post-processing is typically considered a low-order term in
comparison to the cost of computing the sketch to begin with.
When analyzing sketching algorithms it is common to separate the post-processing step from the
dimensionality reduction step. Known post-processing routines give good approximate solutions
to linear algebra problems under the condition that C satisfies certain approximation properties
with respect to A. The challenge then becomes proving that a specific dimensionality reduction
algorithm produces a sketch satisfying these required guarantees.

2.3

Sampling Guarantees for Low-rank Approximation

For low-rank approximation, most algorithms aim for one of two standard approximation guarantees, which we describe below. Since we will be focusing on sampling methods, from now on we
assume that C is a subset of A’s columns.
′

Definition 2 (Rank k Column Subset Selection). For d′ < d, a subset of A’s columns C ∈ Rn×d
′
is a (1 + ǫ) factor column subset selection if there exists a rank k matrix Q ∈ Rd ×d with
kA − CQk2F ≤ (1 + ǫ)kA − Ak k2F .
4

(2)

In other words, the column span C contains a good rank k approximation for A. Algorithmically,
we can recover this low-rank approximation via projection to the column subset [Sar06, CW13].
Beyond sketching for low-rank approximation, the column subset selection guarantee is used as
a metric in feature selection for high dimensional datasets [PZB+ 07, MD09]. With columns of A
interpreted as features and rows as data points, (2) ensures that we select d′ features that span the
feature space nearly as well as the top k principal components. The guarantee is also important in
algorithms for CUR matrix decomposition [DKM06b, MD09, BW14] and Nyström approximation
[WS01, MD05, BW09b, BW09a, GM13, HI15, MM16].
In addition to Definition 2, we consider a stronger guarantee for weighted column selection,
which has a broader range of algorithmic applications:
Definition 3 (Rank k Projection-Cost Preserving Sample). For d′ < d, a subset of rescaled columns
′
C ∈ Rn×d is a (1 + ǫ) projection-cost preserving sample if, for all rank k orthogonal projection
matrices X ∈ Rn×n ,
(1 − ǫ)kA − XAk2F ≤ kC − XCk2F ≤ (1 + ǫ)kA − XAk2F

(3)

Definition 3 is formalized in two recent papers [FSS13, CEM+ 15], though it appears implicitly
in prior work [DFK+ 04, BZMD15]. It ensures that C approximates the cost of any rank k column
projection of A. C can thus be used as a direct surrogate of A to solve low-rank projection problems.
Specifically, it’s not hard to see that if we use a post-processing algorithm that sets Z equal to the
top k left singular vectors of C, it will hold that kA−ZZT AkF ≤ (1+ǫ)kA−Uk UTk AkF [CEM+ 15].
Definition 3 also ensures that C can be used in approximately solving a variety of constrained
low-rank approximation problems, including k-means clustering of A’s rows (see [CEM+ 15]).

2.4

Leverage Scores

It is well known that sketches satisfying Definitions 2 and 3 can be constructed via importance
sampling routines which select columns using carefully chosen, non-uniform probabilities. Many of
these probabilities are modifications on traditional “statistical leverage scores”.
The statistical leverage score of the ith column ai of A is defined as4 :
def

τi = aTi (AAT )+ ai .

(4)

τi measures how important ai is in composing the range of A. It is maximized at 1 when ai is linearly
independent from A’s other columns and decreases when many other columns approximately align
with ai or when kai k2 is small.
Leverage scores are used in fast sketching algorithms for linear regression and matrix preconditioning [DMM06, Sar06, CLM+ 15]. They have also been applied to convex optimization [LSW15],
linear programming [LS14, LS15], matrix completion [CBSW15], multi-label classification [BK13],
and graph sparsification, where they are known as effective resistances [SS11].

2.5

Existing Low-rank Leverage Scores

For low-rank approximation problems, leverage scores need to be modified to only capture how
important each column ai is in composing the top few singular directions of A’s range.
4

+ denotes the Moore-Penrose pseudoinverse of a matrix. When AAT is full rank (AAT )+ = (AAT )−1

5

d′

In particular, it is known that a sketch C satisfying Definition 2 can be constructed by sampling
= O(k log k + k/ǫ) columns according to the so-called rank k subspace scores [Sar06, DMM08]:
(k) def

ith rank k subspace score:

ssi

= aTi (Ak ATk )+ ai .

(5)

These scores are exactly equivalent to standard leverage scores computed with respect to Ak , an
optimal low-rank approximation for A. The stronger projection-cost preservation guarantee of
Definition 3 can be achieved by sampling O(k log k/ǫ2 ) columns using a related, but somewhat
more complex, leverage score modification [CEM+ 15].

2.6

Ridge Leverage Scores

Notably, prior low-rank leverage scores are defined in terms of Ak , which is not always unique and
regardless can be sensitive to matrix perturbations5 . As a result, the scores can change drastically
when A is modified slightly or when only partial information about the matrix is known. This
largely limits the possibility of quickly approximating the scores with sampling algorithms, and
motivates our adoption of a new leverage score for low-rank approximation.
Rather than use scores based on Ak , we employ regularized scores called ridge leverage scores,
which have been used for approximate kernel ridge regression [AM15] and in work on iteratively
computing standard leverage scores [LMP13, KLM+ 14]. We extend their applicability to low-rank
approximation. For a given regularization parameter λ, define the λ-ridge leverage score as:
def

τiλ (A) = aTi (AAT + λI)+ ai .

(6)

We will always set λ = kA − Ak k2F /k and thus, for simplicity, use “ith ridge leverage score” to refer
+

kA−Ak k2F
I
ai .
to τ̄i (A) = aTi AAT +
k
For prior low-rank leverage scores, Ak truncates the spectrum of A, removing all but its top k
singular values. Regularization offers a smooth alternative: adding λI to AAT ‘washes out’ small
singular directions, causing them to be sampled with proportionately lower probability.
This paper proves that regularization can not only replace truncation, but is more natural and
stable. In particular, while τ̄i depends on the value of kA − Ak k2F , it does not depend on a specific
low-rank approximation. This is sufficient for stability since kA − Ak k2F changes predictably under
matrix perturbations even when Ak itself does not.
Before showing our sampling guarantees for ridge leverage scores, we prove that the sum of
these scores is not too large. Thus, when we use them for sampling, we will achieve column subsets
and projection-cost preserving samples of small size. Specifically we have:
Pn
Lemma 4.
i=1 τ̄i (A) ≤ 2k.
Proof. We rewrite (6) using A’s singular value decomposition:
τ̄i (A) =

aTi

= aTi
5

!+
kA\k k2F
T
UU
ai
UΣ U +
k
+

UΣ̄2 UT ai = aTi UΣ̄−2 UT ai ,
2

T

It is often fine to use a near-optimal low-rank approximation in place of Ak , but similar instability issues remain.

6

where Σ̄2i,i = σi2 (A) +
n
X
i=1

(Σ2 Σ̄−2 )i,i =

2

tr(Σ Σ̄

kA−Ak k2F
k



τ̄i (A) = tr AT UΣ̄−2 UT A = tr VΣΣ̄−2 ΣVT = tr(Σ2 Σ̄−2 )

σi2 (A)
σi2 (A)+

−2

kA−Ak k2
F
k

)=k+

n
X

i=k+1

3

. We then have:

. For i ≤ k we simply upper bound this by 1. So:
σi2
σi2 +

kA−Ak k2F
k

≤k+

n
X

σi2

i=k+1

kA−Ak k2F
k

=k+

Pn

2
i=k+1 σi
kA−Ak k2F
k

≤ k + k.

Core Sampling Results

Before considering how to efficiently compute ridge leverage scores, we prove that they can be used
to construct sketches satisfying the guarantees of Definitions 2 and 3. To do so, we introduce a
natural intermediate guarantee (Theorem 5), from which our results on column subset selection and
projection-cost preservation follow. This approach is the first to treat these guarantees in a unified
way and we hope it will be useful in future work on sketching methods for low-rank approximation.
Specifically, we will show that our selected columns spectrally approximate A up to additive error
depending on the ridge parameter λ = kA − Ak kF /k. This approximation is akin to the ubiquitous
subspace embedding guarantee [Sar06] which is used as a primitive for full rank problems like linear
regression and generally requires sampling Θ(d) columns.
Intuitively,
√ sampling by ridge leverage scores is equivalent to sampling by the standard leverage
scores of [A, λIn×n ]. A matrix Chernoff bound can be used to show that sampling by these scores
will yield C satisfying the subspace embedding property: (1 − ǫ)CCT  AAT + λI  (1 + ǫ)CCT .
(Recall that M  N indicates that sT Ms ≤ sT Ns for every vector s.)
However, we do not actually sample columns of the identify, only columns of A. Subtracting
off the identity yields the mixed additive-multiplicative bound of Theorem 5.
Theorem 5 (Additive-Multiplicative Spectral Approximation). For i ∈ {1, . . . , d}, let τ̃i ≥ τ̄i (A)
P
be an overestimate for the ith ridge leverage score. Let pi = Pτ̃i τ̃i . Let t = c log(k/δ)
i τ̃i for some
ǫ2
i

sufficiently large constant c. Construct C by sampling t columns of A, each set to
probability pi . With probability 1 − δ, C satisfies:

√1 ai
tpi

with

ǫ
ǫ
(1 − ǫ)CCT − kA − Ak k2F In×n  AAT  (1 + ǫ)CCT + kA − Ak k2F In×n
(7)
k
k


columns.
By Lemma 4, if each τ̃i is within a constant factor of τ̄i (A) then C has O k log(k/δ)
2
ǫ
Note that Theorem 5 and our other sampling results hold for independent sampling without replacement. A proof is included in Appendix B.

Proof. Following Lemma 4, we have τ̄i (A) = aTi UΣ̄−2 UT ai , where Σ̄2i,i = σi2 (A) +

Let Y = Σ̄−1 UT CCT − AAT UΣ̄−1 . We can write
Y=

t 
X
j=1

Σ̄

−1

T

U



cj cTj

1
− AAT
t
7



UΣ̄

−1



def

=

t
X
j=1

[Xj ] .

kA\k k2F
k

.

For each j ∈ 1, . . . , t, Xj is given by:


1
1
−1 T
T
T
Xj = · Σ̄ U
ai ai − AA
UΣ̄−1 with probability pi .
t
pi
i
h
E Y = 0 since E p1i ai aTi − AAT = 0. Furthermore, CCT = UΣ̄Y Σ̄U + AAT . Showing

kYk2 ≤ ǫ gives −ǫI  Y  ǫI, and since UΣ̄2 UT = AAT +
(1 − ǫ)AAT −

kA\k k2F
k

I would give:

ǫkA\k k2F
ǫkA\k k2F
I  CCT  (1 + ǫ)AAT +
I.
k
k

After rearranging and adjusting constants on ǫ, this statement is equivalent to (7).
To prove that kYk2 is small with high probability we use a stable rank (intrinsic dimension)
matrix Bernstein inequality from [Tro15] that was first proven in [Min13] following work in [HKZ12].
This inequality requires upper bounds on the spectral norm of each Xj and on variance of Y.
kA

k2

1
F
I. This is a well known property of
We use the fact that, for any i, τ̄i (A)
ai aTi  AAT + \k
k
leverage scores, shown for example in the proof of Lemma 11 in [CLM+ 15]. It lets us bound:
!
2
kA
k
1
\k F
· Σ̄−1 UT ai aTi UΣ̄−1  Σ̄−1 UT AAT +
I UΣ̄−1 = I.
τ̄i (A)
k

So we have:
1
1
ǫ2
P
Xj + Σ̄−1 UT AAT UΣ̄−1 
·
· τ̄i (A) · I 
t
tpi
c log(k/δ) i τ̃i

P

i τ̃i

τ̃i

· τ̄i (A) · I 

ǫ2
I.
c log(k/δ)

Additionally,

ǫ2
ǫ2
1 −1 T
P
Σ̄ U AAT UΣ̄−1 =
· Σ̄−2 Σ2 
I,
t
c log(k/δ) i τ̃i
c log(k/δ)

where the inequality follows from the fact that:
X
X



τ̄i (A) = tr AT UΣ̄−2 UT A = tr UΣ̄−2 Σ2 UT = tr Σ̄−2 Σ2 ≥ kΣ̄−2 Σ2 k2 .
τ̃i ≥
i

i

Overall this gives kXj k2 ≤
E(Y 2 ) = t · E(X2j ) =

ǫ2
c log(k/δ) .

1X
pi ·
t



Next we bound the variance of Y.
1 −1 T
Σ̄ U ai aTi UΣ̄−2 UT ai aTi UΣ̄−1
p2i


1 −1 T
T
−2 T
T
−1
−1 T
T
−2 T
T
−1
−2 Σ̄ U ai ai UΣ̄ U AA UΣ̄ + Σ̄ U AA UΣ̄ U AA UΣ̄
pi
P

τ̃i
1
1X
−1 T
T
−1
· τ̄i (A) · Σ̄ U ai ai UΣ̄
− Σ̄−1 UT AAT UΣ̄−2 UT AAT UΣ̄−1

t
τ̃i
t
2
ǫ

Σ̄−1 UT AAT UΣ̄−1
c log(k/δ)
ǫ2
ǫ2
Σ2 · Σ̄−2 
D,
(8)

c log(k/δ)
c log(k/δ)
8

where we set Di,i = 1 for i ∈ 1, . . . , k and Di,i =

σi2

σi2 +

kA\k k2
F
k

for all i ∈ k + 1, ..., n. By the stable

rank matrix Bernstein inequality given in Theorem 7.3.1 of [Tro15], for ǫ < 1,
P [kYk ≥ ǫ] ≤

4 tr(D)
·e
kDk2



−ǫ2 /2

ǫ2
(kDk2 +ǫ/3)
c log(k/δ)

.

(9)

Clearly kDk2 = 1. Furthermore, following Lemma 4, tr(D) ≤ 2k. Plugging into (9), we see that
P [kYk ≥ ǫ] ≤ 8ke−

c log(k/δ)
2

) ≤ δ/2,

if we choose the constant c large enough. So we have established (7).

3.1

Projection-Cost Preserving Sampling

We now use Theorem 5 to prove that sampling by ridge leverage scores is sufficient for constructing projection-cost preserving samples. The following theorem is a basic building block in our
O(nnz(A)) time low-rank approximation algorithm.
Theorem 6 (Projection-Cost Preservation). For i ∈ {1, . . . , d}, let τ̃i ≥ τ̄i (A) be an overestimate
P
for the ith ridge leverage score. Let pi = Pτ̃i τ̃i . Let t = c log(k/δ)
i τ̃i for any ǫ < 1 and some
ǫ2
i

sufficiently large constant c. Construct C by sampling t columns of A, each set to
probability pi . With probability 1 − δ, for any rank k orthogonal projection X,

√1 ai
tpi

with

(1 − ǫ)kA − XAk2F ≤ kC − XCk2F ≤ (1 + ǫ)kA − XAk2F .
Note that the theorem also holds for independent sampling without replacement, as shown in
Appendix B. By Lemma 4, when each approximation τ̃i is within a constant factor of the true ridge
leverage score τ̄i (A), we obtain a projection-cost preserving sample with t = O(k log(k/δ)/ǫ2 ).
To simplify bookkeeping, we only worry about proving a version of Theorem 6 with (1 ± aǫ)
error for some constant a, and assume ǫ ≤ 1/2. By simply adjusting our constant oversampling
parameter, c, we can recover the result as stated.
The challenge in proving Theorem 6 comes from the mixed additive-multiplicative error of
Theorem 5. Pure multiplicative error, e.g. from a subspace embedding, or pure additive error, e.g.
from a “Frequent Directions” sketch [GLPW15], are easily converted to projection-cost preservation
results [Mus15], but merging the analysis is intricate. To do so, we split AAT and CCT into their
projections onto the top “head” singular vectors of A and onto the remaining “tail” singular vectors.
Restricted to the span of A’s top singular vectors, Theorem 5 gives a purely multiplicative bound.
Restricted to vectors spanned by A’s lower singular vectors, the bound is purely additive.
Proof. For notational convenience, let Y denote I − X, so kA − XAk2F = tr(YAAT Y) and kC −
XCk2F = tr(YCCT Y).
3.1.1

Head/Tail Split

2 ≥ kA − A k2 /k. Let P
Let m be the index of the smallest singular value of A such that σm
m
k F
T
T
denote Um Um and P\m denote U\m U\m = I − Pm . We split:

tr(YAAT Y) = tr(YPm AAT Pm Y) + tr(YP\m AAT P\m Y) + 2 tr(YPm AAT P\m Y)
= tr(YAm ATm Y) + tr(YA\m AT\m Y).
9

(10)

The “cross terms” involving Pm A and P\m A equal 0 since the two matrices have mutually orthogT and VT , respectively). Additionally, we split:
onal rows (spanned by Vm
\m
tr(YCCT Y) = tr(YPm CCT Pm Y) + tr(YP\m CCT P\m Y) + 2 tr(YPm CCT P\m Y)

(11)

In (11) cross terms do not cancel because, in general, Pm C and P\m C will not have orthogonal
rows, even though they have orthogonal columns. Regardless, while these terms make our analysis
more difficult, we proceed with comparing corresponding parts of (10) and (11).
3.1.2

Head Terms

We first bound the terms involving Pm , beginning by showing that:
1+ǫ
1−ǫ
Pm CCT Pm  Am ATm 
Pm CCT Pm .
1+ǫ
1−ǫ

(12)

For any vector x, let y = Pm x. Note that xT Am ATm x = yT AAT y since Am ATm = Pm AAT Pm
and since Pm Pm = Pm . So, using (7) we can bound:
(1 − ǫ)yT CCT y − ǫ

kA\k k2F T
kA\k k2F T
y y ≤ xT Am ATm x ≤ (1 + ǫ)yT CCT y + ǫ
y y.
k
k

(13)

By our definition of m, y is orthogonal to all singular directions of A except those with squared
singular value greater than or equal to kA\k k2F /k. It follows that
xT Am ATm x = yT AAT y ≥

kA\k k2F T
y y,
k

and accordingly, from the left side of (13), that (1 − ǫ)yT CCT y ≤ (1 + ǫ)xT Am ATm x. Additionally,
from the right side of (13), we have that (1 + ǫ)yT CCT y ≥ (1 − ǫ)xT Am ATm x. Since yT CCT y =
xT Pm CCT Pm x, these inequalities combine to prove (12). From (12) we can bound the diagonal
entries of YAm ATm Y in terms of the corresponding diagonal entries of YPm CCT Pm Y, which are
all positive, and conclude that:
1−ǫ
1+ǫ
tr(YPm CCT Pm Y) ≤ tr(YAm ATm Y) ≤
tr(YPm CCT Pm Y).
1+ǫ
1−ǫ
Assuming ǫ < 1/2, this is equivalent to:
(1 − 4ǫ) tr(YAm ATm Y) ≤ tr(YPm CCT Pm Y) ≤ (1 + 4ǫ) tr(YAm ATm ).
3.1.3

(14)

Tail Terms

For the lower singular directions of A, Theorem 5 does not give a multiplicative spectral approximation, so we do things a bit differently. Specifically, we start by noting that:
tr(YA\m AT\m Y) = tr(A\m AT\m ) − tr(XA\m AT\m X) and

tr(YP\m CCT P\m Y) = tr(P\m CCT P\m ) − tr(XP\m CCT P\m X).
We handle tr(A\m AT\m ) = kA\m k2F and tr(P\m CCT P\m ) = kP\m Ck2F first. Since C is constructed


via an unbiased sampling of A’s columns, E kP\m Ck2F = kA\m k2F and a scalar Chernoff bound
10

is sufficient for showing that this value concentrates around its expectation. Our proof is included
as Lemma 20 in Appendix A and implies the following bound:
−ǫkA\k k2F ≤ tr(A\m AT\m ) − tr(P\m CCT P\m ) ≤ ǫkA\k k2F .

(15)

Next, we compare tr(XA\m AT\m X) to tr(XP\m CCT P\m X). We first claim that:
P\m CCT P\m −

4ǫ
4ǫ
kA\k k2F I  A\m AT\m  P\m CCT P\m + kA\k k2F I.
k
k

(16)

The argument is similar to the one for (12). For a vector x, let y = P\m x. xT A\m AT\m x =
yT AAT y since A\m AT\m = P\m AAT P\m and since P\m P\m = P\m . Applying (7) gives:
kA\k k2F T
kA\k k2F T
T
T
T
T
y y ≤ x A\m A\m x ≤ (1 + ǫ)y CC y + ǫ
y y.
(1 − ǫ)y CC y − ǫ
k
k
T

T

Noting that yT y ≤ xT x and assuming ǫ ≤ 1/2 gives the following two inequalities:
kA\k k2F T
x x ≤ (1 + 2ǫ)xT A\m AT\m x,
k
kA\k k2F T
x x.
(1 − 2ǫ)xT A\m AT\m x ≤ yT CCT y + 2ǫ
k

yT CCT y − 2ǫ

kA

(17)
(18)

k2

\k F
By our choice of m, xT A\m AT\m x ≤
xT x. So, substituting y with P\m x and rearranging
k
(17) and (18) gives (16).
Now, since X is a rank k projection matrix, it can be written as X = ZZT where Z ∈ Rn×k is
a matrix with k orthonormal columns, z1 , . . . , zk . By cyclic property of the trace,

tr(XA\m AT\m X) = tr(ZT A\m AT\m Z) =
Similarly, tr(XP\m CCT P\m X) =

Pk

T
T
i=1 zi P\m CC P\m zi

k
X
i=1

zTi A\m AT\m zi .

and we conclude from (16) that:

tr(XP\m CCT P\m X) − 4ǫkA\k k2F ≤ tr(XA\m AT\m X) ≤ tr(XP\m CCT P\m X) + 4ǫkA\k k2F ,
which combines with (15) to give the final bound:
tr(YA\m AT\m Y) − 5ǫkA\k k2F ≤ tr(YP\m CCT P\m Y) ≤ tr(YA\m AT\m Y) + 5ǫkA\k k2F .
3.1.4

(19)

Cross Terms

Finally, we handle the cross term 2 tr(YPm CCT P\m Y). We do not have anything to compare this
term to, so we just need to show that it is small. To do so, we rewrite:
tr(YPm CCT P\m Y) = tr(YAAT (AAT )+ Pm CCT P\m ),

(20)

which is an equality since the columns of Pm CCT P\m fall in the span of A’s columns. We eliminate
the trailing Y using the cyclic property of the trace. hM, Ni = tr(M(AAT )+ NT ) is a semi-inner

11

product since AAT is positive semidefinite. Thus, by the Cauchy-Schwarz inequality,
| tr(YAAT (AAT )+ Pm CCT P\m )| ≤
q
tr(YAAT (AAT )+ AAT Y) · tr(P\m CCT Pm (AAT )+ Pm CCT P\m )
q
T
T
= tr(YAAT Y) · tr(P\m CCT Um Σ−2
m Um CC P\m )
q
q
2
= tr(YAAT Y) · kP\m CCT Um Σ−1
(21)
m kF .

To bound the second term, we separate:
kP\m CC

T

2
Um Σ−1
m kF

=

m
X
i=1

kP\m CCT ui k22 σi−2 .

(22)

We next show that the summand is small for every i. Take pi to be a unit vector in the direction
of CCT ui ’s projection onto P\m . I.e. pi = P\m CCT ui /kP\m CCT ui k2 . Then:
kP\m CCT ui k22 = (pTi CCui )2 .


√
Now, suppose we construct the vector m = σi−1 ui + kA kkF pi . From (7) we know that:

(23)

\k

(1 − ǫ)mT CCT m −

ǫkA\k k2F
k

mT m ≤ mT AAT m,

which expands to give:
(1 −

σi−2 uTi AAT ui +

√

2 k
pT CCT pi + (1 − ǫ)
+ (1 − ǫ)
pT CCT ui ≤
σi kA\k kF i
kA\k k2F i
ǫkA\k k2F T
ǫkA\k k2F T
k
k
T
T
m
m
=
1
+
m m.
+
p
AA
p
+
i
k
k
kA\k k2F
kA\k k2F i
k

ǫ)σi−2 uTi CCT ui

(24)

There are no cross terms on the right side because pi lies in the span of U\m and is thus orthogonal
to ui over AAT . Now, from (12) we know that uTi CCT ui ≥ (1 − 2ǫ)uTi AAT ui ≥ (1 − 2ǫ)σi2 . From
(16) we also know that pTi CCT pi ≥ pTi AAT pi − 4ǫ

kA\k k2F
k

. Plugging into (28) gives:
√
2 k
k
−2 T
T
T
T
p AA pi − 4ǫ + (1 − ǫ)
(1 − 3ǫ)σi ui AA ui + (1 − ǫ)
pT CCT ui
σi kA\k kF i
kA\k k2F i
≤1+

k

kA\k k2F

pTi AAT pi +

Noting that pTi AAT pi ≤

ǫkA\k k2F T
m m.
k

(25)

kA\k k2F
k

since pi lies in the column span of U\m , rearranging (25) gives:
√
ǫkA\k k2F T
2 k
m m ≤ 12ǫ.
pTi CCT ui ≤ 8ǫ +
(1 − ǫ)
σi kA\k kF
k
 √ 2
√
The second inequality follows from the fact that σi−1 ≤ kA kkF so kmk22 ≤ kA2 kkF . Assuming
\k
\k
again that ǫ ≤ 1/2 gives our final bound:
√
k
pi CCT uTi ≤ 12ǫ
σi kA\k kF
(pi CCT uTi )2 ≤ 144ǫ2
12

σi2 kA\k k2F
.
k

(26)

Plugging into (22) gives:
2
kP\m CCT Um Σ−1
m kF ≤

m
X

144ǫ2

i=1

σi2 kA\k k2F −2
σi ≤ 288ǫ2 kA\k k2F .
k

(27)

Note that we get an extra factor of 2 because m ≤ 2k. Returning to (21), we conclude that:
q
q
| tr(YAAT (AAT )+ Pm CCT P\m )| ≤ tr(YAAT Y) · 288ǫ2 kA\k k2F ≤ 17ǫ tr(YAAT Y). (28)

The last inequality follows from the fact that kA\k k2F ≤ tr(YAAT Y) since A\k is the best rank k
approximation to A. tr(YAAT Y) = kA−XAk2F is the error of a suboptimal rank k approximation.
3.1.5

Final Bound

Ultimately, from(11), (14), (19), and (28), we conclude:
(1 − 4ǫ) tr(YAm ATm Y) + tr(YA\m AT\m Y) − 5ǫkA\k k2F − 34ǫ tr(YAAT Y) ≤ tr(YCCT Y)
≤ (1 + 4ǫ) tr(YAm ATm ) + tr(YA\m AT\m Y) + 5ǫkA\k k2F + 34 tr(YAAT Y).

Applying the fact that kA\k k2F ≤ tr(YAAT Y) proves Theorem 6 for a constant factor of ǫ.

3.2

Column Subset Selection

Although not required for our main low-rank approximation algorithm, we also prove that ridge
leverage score sampling can be used to obtain (1 + ǫ) error column subsets (Definition 2). The
column subset selection problem is of independent interest and the following result allows ridge
leverage scores to be used in our single-pass streaming algorithm for this problem (Section 6).
Theorem 7. For i ∈ {1, .. . , d}, let τ̃i ≥ 
τ̄i (A) be an overestimate for the ith ridge leverage score.
P
Let pi = Pτ̃i τ̃i . Let t = c log k + log(1/δ)
i τ̃i for ǫ < 1 and some sufficiently large constant c.
ǫ
i
Construct C by sampling t columns of A, each set to ai with probability pi . With probability 1 − δ:

kA − CC+ A k k2F ≤ (1 + ǫ)kA − Ak k2F .
P
Furthermore, C contains a subset of O ( i τ̃i /ǫ) columns that satisfies Definition 2 and can be
identified in polynomial time.

Note that (CC+ A)k is a rank k matrix in the column span of C, so Theorem 7 implies that
C is a (1 + ǫ) error column subset according to Definition 2. By Lemma 4, if each τ̃i is within a
constant factor of τ̄i (A), the approximate ridge leverage scores sum to O(k) so Theorem 7 gives a
column subset of size O (k log k + k/ǫ), which contains a near optimally sized column subset with
O (k/ǫ) columns. Again, the theorem also holds for sampling without replacement (see Appx. B).
Our proof relies on establishing a connection between ridge leverage sampling and well known
adaptive sampling techniques for column subset selection [DRVW06, DV06]. We start with the
following lemma on adaptive sampling for column subset selection:

Lemma 8 (Theorem 2.1 of [DRVW06]). Let C be any subset of A’s columns and let Z be
an orthonormal matrix whose
columns span those of C. If we sample an additional set S of

k log(1/δ) kA−ZZT Ak2F
columns from A with probability proportional to k(A − ZZT A)i k22 , then
O
· kA k2
ǫ
\k F

[S ∪ C] is a (1 + ǫ) error column subset for A with probability (1 − δ).

6

6
Theorem 2.1 was originally stated as an expected error result, but it can be seen to hold with constant probability
via Markov’s inequality and accordingly with (1 − δ) probability when oversampling by a factor of log(1/δ)

13

When C is a constant error column subset, then kA − ZZT Ak2F ≤ kA − (ZZT A)k k2F =
O(kA\k k2F ) and accordingly we only need O(k log(1/δ)/ǫ) additional adaptive samples. So one
potential algorithm for column subset selection is as follows: apply Theorems 5 and 6, sampling
O(k log(k/δ)) columns by ridge leverage score to obtain a constant error projection-cost preserving sample, will also be a constant error column subset. Then sample O(k log(1/δ)/ǫ) additional
columns adaptively against C.
However, it turns out that ridge leverage scores well approximate adaptive sampling probabilities computed with respect to any constant error additive-multiplicative spectral approximation
satisfying Theorem 5! That is, surprisingly, they achieve the performance of adaptive sampling
without being adaptive at all. Simply sampling O(k log(1/δ)/ǫ) more columns by ridge leverage
score and invoking Lemma 8 suffices to achieve (1 + ǫ) error.
Proof of Theorem 7. We formally prove that C is itself a good column subset before showing our
stronger guarantee, that it also contains a column subset of optimal size, up to constants.
3.2.1

Primary Column Subset Selection Guarantee

P
We split our sample C, into C
,
which
contains
the
first
c
log(k/δ)
1
i τ̃i columns and C2 , which
P
contains the next c log(1/δ)/ǫ i τ̃i columns. Note that in our final sample complexity the log(1/δ)
factor in the size of C1 is not shown as it is absorbed into the larger size of C2 when log(1/δ) > log(k)
and into the log(k) otherwise. By Theorem 6, we know that, appropriately reweighted, C1 is a
constant error projection-cost preserving sample of A. This means that C1 is also a constant error
column subset. Let Z be an orthonormal matrix whose columns span the columns of C1 .
To invoke Lemma 8 to boost C1 to a (1 + ǫ) column subset, we need to sample columns with
probabilities proportional to k(A − ZZT A)i k22 . This is equivalent to sampling proportional to:

(aTi − aTi ZZT )(ai − ZZT ai ) = aTi ai − 2aTi ZZT ai + aTi ZZT ZZT ai = aTi I − ZZT ai .
We can assume kA\k k2F > 0 or else C1 must fully span A’s columns and we’re done. Scaling τ̄i (A):
kA\k k2F
τ̄i (A) = aTi
k

k
kA\k k2F

T

AA + I

!+

ai .

Since C1 satisfies Theorem 5 with constant error, for large enough constant c1 ,
!
!
Tk
k
kkC
C
k
1 1 2
AAT + I  c1
C1 CT1 + I  c1 I +
ZZT .
kA\k k2F
kA\k k2F
kA\k k2F
Furthermore, I − ZZT  I + cZZT
c1

+

for any positive c so,
!+
!+
Tk
k
kkC
C
2
1
1
AAT + I
ZZT
 I − ZZT .
 I+
kA\k k2F
kA\k k2F

c1 kA\k k2F
c1 kA\k k2F
T
τ̄i (A) ≥ k(A − ZZP
A)i k22 for all i and hence
τ̃i ≥ k(A − ZZT A)i k22 .
k
k
C2 is a set of c log(1/δ)/ǫ · i τ̃i columns sampled with probability proportional to approximate
ridge leverage scores. Consider forming C′2 by setting (C2 )i = 0 with probability:

So

k(A − ZZT A)j(i) k22
c1 kA\k k2F
τ̃j(i)
k

14

,

where j(i) is just the index of the column of A that (C2 )i is equal to. Clearly, if not equal to 0,
each column of C′2 is equal to ai with probability proportional to the adaptive sampling probability
k(A − ZZT A)i k22 . Additionally, in expectation, the number of nonzero columns will be:


!
T
2
X
X
k(A
−
ZZ
A)
k
τ̃
ck log(1/δ) kA − ZZT Ak2F
j 2
Pj
·
=
.
c log(1/δ)/ǫ ·
τ̃i ·
2
c1 kA\k k2F
c1 ǫ
τ̃i
kA
k
\k
i
F
τ̃
j
i
k

j

By a Chernoff bound, with probability 1 − δ/2 at least half this number of columns will be
nonzero, and by Lemma 8, for large enough c, conditioning on the above column count bound
holding, [C1 ∪ C′2 ] is a (1 + ǫ) error column subset for A with probability 1 − δ/2. Just noting that
span([C1 ∪ C′2 ]) ⊆ span([C1 ∪ C2 ]) and union bounding over the two possible fail conditions, gives
that [C1 ∪ C2 ] = C is a (1 + ǫ) column subset with probability at least 1 − δ.
3.2.2

Stronger Containment Guarantee

P
It now remains to show the second condition of Theorem 7: C contains a subset of O( i τ̃i /ǫ)
columns that also satisfies Definition 2. This follows from noting that we can apply, for example,
the polynomial time deterministic column selection algorithm of [CEM+ 15] to produce a matrix
C′1 with O(k) columns that is both a constant error additive-multiplicative spectral approximation
and a constant error projection-cost preserving sample for C1 . If C′1 has constant error for C1 , it
does for A as well and so is a constantPerror column subset.
O(log(1/δ))
. By our argument
C2 contains O(log(1/δ)) sets of O( i τ̃i /ǫ) columns, C12 , C22 , . . . , C2
probability.
So with
above, for each Ci2 , [C′1 , Ci2 ] is a (1+ ǫ) error column subset of A with constant
P
P
i
′
probability 1 − δ, at least one [C1 , C2 ] is good. This set contains just O(k + i τ̃i /ǫ) = O( i τ̃i /ǫ)
columns, giving the theorem.

4

Monotonicity of Ridge Leverage Scores

With our main sampling results in place, we focus on the algorithmic problem of how to efficiently
approximate the ridge leverage scores of a matrix A. In the offline setting, we will show that these
scores can be approximated in O(nnz(A)) time using a recursive sampling algorithm. We will also
show how to compute and sample by the scores in a single-pass column stream.
Both of these applications will require a unique stability property of the ridge leverage scores:
Lemma 9 (Ridge Leverage Score Monotonicity). For any A ∈ Rn×d and vector x ∈ Rn , for every
i ∈ 1, . . . , d we have:
τ̄i (A) ≤ τ̄i (A ∪ x),
where A ∪ x is simply A with x appended as its final column.
This statement is extremely natural, given that leverage scores are meant to be a measure of
importance. It ensures that the importance of a column can only decrease when additional columns
are added to A. While it holds for standard leverage scores, surprisingly no prior low-rank leverage
scores satisfy this property.
We begin by defining the generalized ridge leverage score as the ridge leverage score of a column
estimated using a matrix other than A itself.

15

′

Definition 10 (Generalized Ridge Leverage Score). For any A ∈ Rn×d and M ∈ Rn×d , the ith
generalized ridge leverage score of A with respect to M is defined as:
 



aT MMT + kM−Mk k2F I + a for a ∈ span MMT + kM−Mk k2F I
i
i
i
k
k
τ̄iM (A) =
∞
otherwise.

This definition is the intuitive one. Since our goal is typically to compute over-estimates of
kM−Mk k2F
τ̄i (A) using M, if ai does not fall in the span of MMT +
I we conservatively set its
k
generalized leverage score to ∞ instead of 0. Note that this case only applies when M is rank k
kM−Mk k2F
and thus
I is 0.
k
We now prove a general monotonicity theorem, from which Lemma 9 follows immediately by
setting M = A and A = A ∪ x.
′

Theorem 11 (Generalized Monotonicity Bound). For any A ∈ Rn×d and M ∈ Rn×d with MMT 
AAT we have:
τ̄i (A) ≤ τ̄iM (A).
Proof. We first note that kM − Mk k2F ≤ kA − Ak k2F since, letting Pk be the projection onto the
top k column singular vectors of A, by the optimality of Mk we have:
kM − Mk k2F ≤ k(I − Pk )Mk2F ≤ k(I − Pk )Ak2F = kA − Ak k2F .
Accordingly,
MMT +

kA − Ak k2F
kM − Mk k2F
I  AAT +
I.
k
k
kM−M k2

k F
Let R be a projection matrix onto the column span of MMT +
I. Since for any PSD
k
+
+
matrices B and C with the same column span, B  C implies B  C (see [MA77]) we have:


+
+
kA − Ak k2F
kM − Mk k2F
T
R  R AA +
R.
I
I
R MM +
k
k


kM−Mk k2F
For any ai not lying in span MMT +
I
, τ̄iM (A) = ∞ and the theorem holds trivially.
k
Otherwise, we have Rai = ai and so:


τ̄i (A) =

aTi R



T

kA − Ak k2F
AA +
I
k
T

+

Rai ≤

aTi R



kM − Mk k2F
MM +
I
k
T

+

Rai = τ̄iM (A).

This gives the theorem.

5

Recursive Ridge Leverage Score Approximation

With Theorem 11 in place, we are ready to prove that ridge leverage scores can be approximated in
O(nnz(A)) time. Our work closely follows [CLM+ 15], which shows how to approximate traditional
leverage scores via recursive sampling.

16

5.1

Intuition and Preliminaries

The central idea behind recursive sampling is as follows: if we uniformly sample, for example, 1/2
of A’s columns to form C and compute ridge leverage score estimates with respect to just these
columns, by monotonicity, the estimates will upper bound A’s true ridge leverage scores. While
some of these upper bounds will be crude, we can show that their overall sum is small.
Accordingly, we can use the estimates to sample O(k log k) columns from A to obtain a constant
factor additive-multiplicative spectral approximation by Theorem 5, as well as a constant factor
projection-cost preserving sample by Theorem 6. This approximation is enough to obtain constant
factor estimates of the ridge leverage scores of A.
C may still be relatively large (e.g. half the size of A), but it can be recursively approximated
via the same sampling scheme, eventually giving our input sparsity time algorithm.
We first give a foundational lemma showing that an approximation of the form given by Theorems 5 and 6 is enough to give constant factor approximations to ridge leverage scores.
Lemma 12. Assume that, for an ǫ ≤ 1/2, we have C satisfying equation (7) from Theorem 5:
ǫ
ǫ
(1 − ǫ)CCT − kA − Ak k2F I  AAT  (1 + ǫ)CCT + kA − Ak k2F I,
k
k
along with equation (3) from Definition 3:
(1 − ǫ)kA − XAk2F ≤ kC − XCk2F ≤ (1 + ǫ)kA − XAk2F , ∀ rank k X.
Then for all i,
(1 − 4ǫ)τ̄i (A) ≤ τ̄iC (A) ≤ (1 + 4ǫ)τ̄i (A).
Proof. Let Pk be the projection onto A’s top k column singular vectors. By the optimality of Ck
in approximating C and the projection-cost preservation condition, we know that kC − Ck k2F ≤
kC−Pk Ck2F ≤ (1+ǫ)kA−Ak k2F . Also, letting P̃k be the projection onto C’s top k column singular
vectors, we have (1 − ǫ)kA − Ak k2F ≤ (1 − ǫ)kA − P̃k Ak2F ≤ kC − Ck k2F . So overall:
(1 − ǫ)kA − Ak k2F ≤ kC − Ck k2F ≤ (1 + ǫ)kA − Ak k2F .

(29)

Using the guarantee from Theorem 5 we have:
(1 − ǫ)CCT +

kA − Ak k2F
(1 + ǫ)kA − Ak k2F
(1 − ǫ)kA − Ak k2F
I  AAT +
I  (1 + ǫ)CCT +
I.
k
k
k

Combining with our bound on kC − Ck k2F gives:
T

(1 − ǫ)CC +

(1−ǫ)
(1+ǫ) kC

− Ck k2F

k

kA − Ak k2F
I  AA +
I  (1 + ǫ)CCT +
k
T

(1+ǫ)
(1−ǫ) kC

− Ck k2F

k

I,

and when ǫ ≤ 1/2, we can simplify to:




kA − Ak k2F
kC − Ck k2F
kC − Ck k2F
T
T
T
I  AA +
I  (1 + 4ǫ) CC +
I .
(1 − 4ǫ) CC +
k
k
k
If kA − Ak k2F = 0, and thus by (29) kC − Ck k2F = 0, then A and C must have the same column
span or else it could not hold that (1 − 4ǫ)CCT  AAT  (1 + 4ǫ)CCT . On the other hand, if

17

kA−A k2

kC−C k2F

k
k F
I and CCT +
kA − Ak k2F > 0, and thus by (29) kC − Ck k2F > 0, both AAT +
k
k
n
span all of R . Either way, the two matrices have the same span and so by [MA77] we have:

kA\k k2F
(1 − 4ǫ) AAT +
I
k

!+



kC\k k2F
CCT +
I
k

!+

kA\k k2F
 (1 + 4ǫ) AAT +
I
k

!+

I

,

which gives the lemma.
Our next lemma, which is analogous to Theorem 2 of [CLM+ 15], shows that by reweighting a
small number of columns in A, we can obtain a matrix with all ridge leverage scores bounded by
a small constant, which ensures that it can be well approximated by uniform sampling.
Lemma 13 (Ridge Leverage Score Bounding Column Reweighting). For any A ∈ Rn×d and any
score upper bound u > 0, there exists a diagonal matrix W ∈ Rd×d with 0  W  I such that:
∀i, τ̄i (AW) ≤ u,

(30)

and
|{i : Wii 6= 1}| ≤

3k
.
u

(31)

Proof. This result follows from Theorem 2 of [CLM+ 15], to which we refer the reader for details. To
show the existence of a reweighting W satisfying (30) and (31), we will argue that a simple iterative
process (which we never actually need to implement) converges on the necessary reweighting.
Specifically, if a column has too high of a leverage score, we simply decrease its weight until
τ̄i (AW) ≤ u. We want to argue that, given AW0 with τ̄i (AW0 ) > u, we can decrease the weight
on ai to produce W1 with τ̄i (AW1 ) ≤ u. By Lemma 5 of [CLM+ 15] we can always decrease
the weight on ai to ensure τi (AW1 ) ≤ u, where τi (·) is the traditional leverage score. And since
+


k(AW 1 )\k k2F
2 AT + , τ̄ (AW ) ≤ τ (AW ), so an equivalent or smaller
I
AW12 AT +

AW
i
1
i
1
1
k
weight decrease suffices to decrease τ̄i (AW1 ) below u.
Furthermore, we can see that τ̄i (AW) is continuous with respect to W. This is due to the fact
that both the traditional leverage scores of AW (shown in Lemma 6 of [CLM+ 15]) and k(AW)\k k2F
are continuous in W. From Theorem 2 of [CLM+ 15], continuity implies that iteratively reweighting
individual columns converges, and thus there is always exists a reweighting satisfying (30).
It remains to show that this reweighting satisfies (31). By continuity, we can always decrease
τ̄i (AW0 ) to exactly u unless τ̄i (AW) = 1, in which case the only option is to set the weight on the
column to 0 and hence set τ̄i (AW) = 0. However, if kA\k k2F > 0, then every ridge leverage score
is strictly less than 1. If kA\k k2F = 0, then A has rank k, the ridge leverage scores are the same as
the true leverage scores, and the number of columns with leverage
score 1 is at most k. Therefore,
P
+
by Theorem 2 of [CLM 15], monotonicity, and the fact that i τ̄i (AW) ≤ 2k for any W, we have
the lemma.

5.2

Uniform Sampling for Ridge Leverage Score Approximation

Using Lemmas 12 and 13 we can prove the key step of our recursive sampling method: if we
uniformly sample columns from A and use them to estimate ridge leverage scores, these scores can
be used to resample a set of columns that give constant factor ridge leverage scores approximations.

18

Theorem 14 (Ridge Leverage Score Approximation via Uniform Sampling). Given A ∈ Rn×d ,
construct Cu by independently sampling each column of A with probability 21 . Let
n
o
τ̃i = min 1, τ̄iCu (A) .

If we form C by sampling each column of A independently with probability pi = min {1, τ̃i c1 log(k/δ)}
√
and reweighting by 1/ pi if selected, then for large enough constant c1 , with probability 1 − δ, C
will have just O(k log(k/δ)) columns and will satisfy the conditions of Lemma 12 for some constant
error. Accordingly, we have:
1
τ̄i (A) ≤ τ̄iC (A) ≤ 2τ̄i (A).
2
Cu
Proof. Clearly Cu CTu  AAT , so by the monotonicity shown in
n Theorem o11 we have τ̄i (A) ≥
τ̄i (A). Since τ̄i (A) is always ≤ 1, it follows that τ̃i = min 1, τ̄iCu (A) ≥ τ̄i (A). Then we
can just use the τ̃i ’s obtained from Cu in independent sampling versions of Theorems 5 and 6,
which can be proven from Lemmas 21 and 22 in Appendix B. Accordingly, with probability 1 −
δ/3, C gives a constant factor additive-multiplicative spectral approximation and projection-cost
preserving sample of A. Hence by Lemma 12, τ̄iC (A) is a constant factor approximation to τ̄i (A).
To prove the theorem, we still have to show that C does not have too many columns. Its
expected number of columns is:
X
X
pi =
min {1, τ̃i c1 log(k/δ)} .
i

i

1
By Lemma 13 instantiated with u = 2c2 log(k/δ)
, we know that there is some reweighting matrix W
1
for all i. We have:
with only 3k · 2c2 log(k/δ) entries not equal to 1 such that τ̄i (AW) ≤ 2c2 log(k/δ)

X
i

pi =

X

i:Wii 6=1

pi +

X

pi

i:Wii =1

≤ 6kc2 log(k/δ) +

X

i:Wii =1

c log(k/δ) · τ̄iCu (A)

= 6kc2 log(k/δ) + c1 log(k/δ) ·
≤ 6kc2 log(k/δ) + c1 log(k/δ) ·
≤ 6kc2 log(k/δ) + c1 log(k/δ) ·

X

τ̄iCu (AW)

i:Wii =1

X

τ̄iCu W (AW)

i:Wii =1

X

τ̄iCu W (AW).

(32)

i

1
Now, since every ridge leverage score of AW is bounded by 2c2 log(k/δ)
, if c2 is set large enough,
the uniformly sampled Cu W is a proper ridge leverage score oversampling of AW, except that its
columns were not reweighted by a factor of 2 (they were each sampled with probability 1/2).
Accordingly, with probability 1 − δ/3, 2Cu W satisfies the approximation conditions of Lemma
12 for AW with ǫ = 1/2. Thus, for all i, 21 τ̄iCu W (AW) = τ̄i2Cu W (AW) ≤ 3τ̄i (AW). By Lemma
P
P
4, i τ̄i (AW) ≤ 2k so overall i τ̄iCu W (AW) ≤ 12k. Plugging back in to (32), we conclude that
C has O(k log(k/δ)) columns in expectation, and actually with probability 1 − δ/3 by a Chernoff
bound. Union bounding over our failure probabilities gives the theorem.

19

5.3

Basic Recursive Algorithm

Theorem 14 immediately proves correct Algorithm 1 for ridge leverage score approximation:
Algorithm 1 Repeated Halving
input: A ∈ Rn×d
output: A reweighted column sample C ∈ Rn×O(k log(k/δ)) satisfying the guarantees of Theorems
5 and 6 with constant error.
1: Uniformly sample d2 columns of A to form Cu
2: If Cu has > O(k log k) columns, recursively apply Repeated Halving to compute a constant
factor approximation C̃u for Cu with O(k log k) columns.
3: Compute generalized ridge leverage scores of A with respect to C̃u
4: Use these estimates to sample columns of A to form C
5: return C
Note that, by Lemma 12, generalized ridge leverage scores computed with respect to C̃u are
constant factor approximations to generalized ridge leverage scores computed with respect to Cu .
Accordingly, by Theorem 14, we conclude that C is a valid ridge leverage score sampling of A.
Before giving our full input sparsity time result, we warm up with a simpler theorem that
obtains a slightly suboptimal runtime.
Lemma 15. A simple implementation of Algorithm 1 that succeeds with probability 1 − δ runs in
O (nnz(A) log(d/δ)) + Õ(nk2 ) time.
For clarity of exposition, we use Õ(·) to hide log factors in k, d, and 1/δ on the lower order term.
Proof. The algorithm has log(d/k) levels of recursion and, since we sample our matrix uniformly,
nnz(A) is cut approximately in half at each level, with high probability. It thus suffices to show
that the work done at the top level is O (nnz(A) log(d/δ)) + Õ(nk 2 ).
To compute the generalized ridge leverage scores of A with respect to C̃u we must (approximately) compute, for each ai ,
aTi

kC̃u − (C̃u )k k2F
I
C̃u C̃Tu +
k

!+

ai .

(33)



kC̃ −(C̃ ) k2
We are going to ignore that C̃u C̃Tu + u k u k F I could be sparse and well conditioned (and
thus ideal for iterative solvers) and use direct methods for simplicity.
kC̃ −(C̃ ) k2

u
u k F
and let R ∈ Rn×Õ(k) be an orthonormal basis containing the left
Let λ denote
k
singular vectors of C̃u . We can rewrite:



C̃u C̃Tu + λI = C̃u C̃Tu + λRRT + λ I − RRT ,

and accordingly, using the fact that RRT and (I − RRT ) are orthogonal,


+ 
+ 1

I − RRT .
C̃u C̃Tu + λI
= C̃u C̃Tu + λRRT
+
λ

20

Now, using an SVD of C̃u , which can be computed in Õ(nk 2 ) time, we compute λ and then
+

write C̃u C̃Tu + λRRT
as RΣ−2 RT for some diagonal matrix Σ ∈ RÕ(k)×Õ(k) . Accordingly, to
evaluate (33), we need just need to compute:






1
1
T
T −1 T
T
−2 T
T
I − RR
ai = k R Σ R + √ I − RR
ai RΣ R +
ai k22 .
λ
λ

Since R has Õ(k) columns, naively evaluating this norm for all of A’s columns would require a
total of Õ(nnz(A)k) time. However, we can accelerate the computation via a Johnson-Lindenstrauss
embedding technique thathas become standard for computing
regular leverage scores [SS11].

1
T
−1
T
T
√
Specifically, denoting R Σ R + λ I − RR
as M, we can embed M’s columns into

O(log(d/δ)) × n dimensions by multiplying on the left by a matrix Π ∈ RO(log(d/δ)×n with scaled
random Gaussian or random sign entries. Even though M is n×n, we can perform the multiplication
in Õ(nk log(d/δ)) by working with our factored form of the matrix.
By standard Johnson-Lindenstrauss results, kΠMai k22 will be within a constant factor of
kMai k22 for all i with probability 1 − δ. Furthermore, we can evaluate kΠMai k22 for all ai in
O (nnz(A) log(d/δ)) total time.
 Our final cost for approximating all ridge leverage scores is thus
O (nnz(A) log(d/δ)) + Õ nk2 time, which gives the lemma.

5.4

True Input-Sparsity Time

Sharpening Lemma 15 to eliminate log factors on the nnz(A) runtime term requires standard
optimizations for approximating leverage scores with respect to a subsample [LMP13, CLM+ 15].
In particular, we can actually apply a Johnson-Lindenstrauss embedding matrix to M with just
θ −1 rows for some small constant θ. Doing so will approximate each ridge leverage score to within
a factor of dθ with high probability (see Lemma 4.5 of [LMP13] for example).

This level of approximation is sufficient to resample O kdθ log(k/δ) columns from A to form
an approximation C′ that satisfies the guarantees of Theorems 5 and 6. To form C, we further
sample C′ down to O(k log(k/δ)) columns using its ridge leverage scores, which takes Õ(nk 2 d2θ )
time. Finally, under the reasonable assumption that ǫ and δ are poly(n), we can also assume
d = poly(n). Otherwise, nnz(A) ≥ d dominates the Õ(nk2 d2θ ) term. This yields the following:
Lemma 16. An optimized
implementation of Algorithm 1, succeeding with probability 1 − δ, runs

−1
in time O θ nnz(A) + Õ(n1+θ k2 ) time, for any θ ∈ (0, 1].

Once we have used Algorithm 1 to obtain C satisfying the guarantees of Theorems 5 and 6 with
constant error, we can approximate A’s ridge leverage scores and resample one final time to obtain
an ǫ error projection-cost preserving sketch. This immediately yields our main algorithmic result:

Theorem 1. For any
 1], there exists an iterative column sampling algorithm that, in time
 θ1+θ∈ 2(0,

n
k
−1
, returns Z ∈ Rn×k satisfying:
O θ nnz(A) + Õ
ǫ4
kA − ZZT Ak2F ≤ (1 + ǫ)kA − Ak k2F .

(34)

All significant linear algebraic operations of the algorithm involve matrices whose columns are subsets of those of A, and thus inherit any structure from the original matrix, including sparsity.
Proof.
 θ/2We use the same technique as Lemma 16, but in the last round of sampling we select
columns to obtain an O(ǫ) factor projection-cost preserving sample, C. Setting Z
O kn log(k/δ)
ǫ2
to the top k column singular vectors of C, which takes Õ(n1+θ k2 /ǫ4 ) time, gives (34) [CEM+ 15].
21

6

Streaming Ridge Leverage Score Sampling

We conclude with an application of our results to novel low-rank sampling algorithms for singlepass column streams. While random projection algorithms work naturally in the streaming setting, the study of single-pass streaming column sampling has been limited to the “full-rank” case
[KL13, CMP15, KLM+ 14]. Column subset selection algorithms based on simple norm sampling are
adaptable to streams, but do not give relative error approximation guarantees [DKM06a, FKV04].
Relative error algorithms are obtainable by combining our projection-cost preserving sampling
procedures with the “merge-and-reduce” framework for coresets [BS80, AHPV04, HPM04]. This
approach relies on the composability of projection-cost preserving samples: a (1 + ǫ) error sample
for A unioned with a (1 + ǫ) error sample for B gives a (1 + ǫ) error sample for [A, B] [FSS13].
However, merge-and-reduce requires storage of O(log4 dk/ǫ2 ) scaled columns from A, where d is
the length of our stream (and its value is known ahead of time).
Our algorithms eliminate the logc d stream length dependence, storing a fixed number of columns
that only depends on ǫ and k. We note that our space bounds are given in terms of the number of
real numbers stored. We do not bound the required precision of these numbers, which would include
at least a single logarithmic dependence on d. In particular, we employ a Frequent Directions sketch
that requires words with at least Θ(log(nd)) bits of precision. Rigorously bounding maximum wordsize required for Frequent Directions and our algorithms could be an interesting direction for future
work.

6.1

General Approach

The basic idea behind our algorithms is quite simple and follows intuition from prior work on
standard leverage score sampling [KL13]. Suppose we have some space budget t for storing a
column sample C. As soon as we have streamed in t columns, we can downsample by ridge
leverage scores to say t/2 columns. As more columns are received, we will eventually reach our
storage limit and need to downsample columns again. Doing so naively would compound error: if
we resampled r times, our final sample would have error (1 + ǫ)r .
However, we can avoid compounding error by exploiting Lemma 9, which ensures that, as new
columns are added, the ridge leverage scores of columns already seen only decrease. Whenever we
add a column to C, we can record the probability it was kept with. In the next round of sampling,
we only discard that column with probability equal to the proportion that its ridge leverage score
decreased by (or keep the column with probability 1 if the score remained constant). New columns
are simply sampled by ridge leverage score. This process ensures that, at any point in the stream,
we have a set of columns sampled by true ridge scores with respect to the matrix seen so far.
Accordingly, we will have a (1 + ǫ) error column subset or projection-cost preserving sample at the
end of the stream.
This overview hides a number of details, the most important of which is how to compute ridge
leverage scores at any given point in the stream with respect to the columns of A observed so far. We
do not have direct access to these columns since we have only stored a subset of them. We could use
the fact that our current sample is projection-cost preserving and can be used to approximate ridge
leverage scores (see Lemma 12). However, this approach would introduce sampling dependencies
between columns and would require a logarithmic dependence on stream length to ensure that our
approximation does not fail at any round of sampling.

22

6.2

Frequent Directions for Approximating Ridge Leverage Scores

Instead, we use a constant error deterministic “Frequent Directions” sketch to estimate ridge leverage scores. Introduced in [Lib13] and further analyzed in a series of papers culminating with
[GLPW15], Frequent Directions sketches are easily maintained in a single-pass column stream of
A. The sketch always provides an approximation B ∈ Rn×(ℓ+1)k guaranteeing:
BBT  AAT  BBT +

1 kA\k k2F
.
ℓ
k

(35)

B does not contain columns from A, so it could be dense even for a sparse input matrix. However,
we will only be setting ℓ to a small constant. Precise information about A will be stored in our
column sample C, which maintains sparsity.
We first show that B can be used to compute constant factor approximations to the ridge
leverage scores of A.
Lemma 17. For every column ai ∈ A, define
+

kAk2F − kBk k2F
def T
T
ai .
τ̃i = ai BB +
I
k
If B ∈ Rn×3k is a Frequent Directions sketch for A with accuracy parameter ℓ = 2, then
1
τ̄i (A) ≤ τ̃i ≤ 2τ̄i (A).
2

kAk2F is obviously computable in a single-pass column stream, so τ̃i can be evaluated in the
streaming setting as long as we have access to ai .
Proof. By the Frequent Directions guarantee, either BBT = AAT giving the lemma trivially, or
kA k2F
I and
kA\k k2F ≥ 0. In this case, since BBT  AAT , kAk2F − kBk k2F > 0. So both AAT + \k
k

2 +
2 −kB k2
kA
k
kAk
\k
k
T
F
F
F
BBT +
I span all of Rn . Recalling that τ̄i (A) = aTi AA +
I ai , to prove
k
k
the lemma it suffices to show:




kA\k k2F
kAk2F − kBk k2F
kAk2F − kBk k2F
1
T
T
T
I  AA +
I  2 BB +
I .
(36)
BB +
2
k
k
k
Recall that the squared Frobenius norm of a matrix is equal to the sum of its squared singular
values. Additionally, a standard property of the relation M  N is that, for all i, the ith singular
value σi (M) ≤ σi (N). From the right hand side of (35) it follows that, when ℓ = 2, σi2 (B) ≥
σi2 (A) −

kA\k k2F
2k

. Accordingly, since kBk k2F is the sum of the top k singular values of B,

kA\k k2F
≤ 1.5kA\k k2F .
−
≤
+k
2k


kA k2F
kAk2F −kBk k2F
Since BBT  AAT , it follows that that BBT +
I
 AAT + 1.5 \k
I, which is
k
k
more than tight enough to give the left hand side of (36).
kAk2F

kBk k2F

kA\k k2F

Furthermore , kAk2F − kBk k2F ≥ kA\k k2F , and since ℓ = 2, BBT  AAT −


kA\k k2F
kAk2F − kBk k2F
T
I  AAT +
,
BB +
k
2k
which is more than tight enough to give the right hand side of (36).
23

kA\k k2F
2k

. Overall,

6.3

Streaming Column Subset Selection

Lemma 17 gives rise to a number of natural algorithms for rejection sampling by ridge leverage score.
The simplest approach is to emulate sampling columns from A independently without replacement
(see Lemmas 21 and 22). However, since sampling without replacement produces a variable number
of samples, this method would require a log d dependence to ensure that our space remains bounded
throughout the algorithm’s execution with high probability.
Instead, we apply our “with replacement” bounds, which sample a fixed number of columns, t.
We start by describing Algorithm 2 for column subset selection. The constant c used below is the
necessary oversampling parameter from Theorem 7. C ∈ Rn×t stores our actual column subset and
D ∈ Rn×t stores a queue of new columns. B is a Frequent Directions sketch with parameter ℓ = 2.
Algorithm 2 Streaming Column Subset
input: A ∈ Rn×d , accuracy ǫ, success probability (1 − δ)
output: C ∈ Rn×t such that t =h32c(k log k + k log(1/δ)/ǫ) and each column
i ci is equal to
1 τ̃j c(k log k+k log(1/δ)/ǫ) τ̃j c(k log k+k log(1/δ)/ǫ)
,
column aj with probability pj ∈ 2
and 0 otherwise,
t
t
Pn
where τ̃j ≥ 2τ̄j (A) for all j and j=1 τ̃j ≤ 16k.
1: count := 1, C := 0n×t , D := 0n×t , f robA := 0
⊲ Initialize storage
old
old
⊲ Initialize sampling probabilities
2: [τ̃1 , ..., τ̃t ] := 1
3: for i := 1, . . . , d do
⊲ Process column stream
4:
B := FreqDirUpdate(B, ai )
5:
if count ≤ t then
⊲ Collect t new columns
6:
dcount := ai .
⊲ Update kAk2F
7:
f robA := f robA + kai k22
8:
count := count + 1
9:
else
⊲ Prune columns

10:
[τ̃1 , ..., τ̃t ] := min [τ̃1old , ..., τ̃told ], ApproximateRidgeScores(B, C, f robA)
11:
[τ̃1D , ..., τ̃tD ] := ApproximateRidgeScores(B, D, f robA)
12:
for j := 1, . . . , t do
⊲ Rejection sample
13:
if cj 6= 0 then


14:

15:

16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
1:
2:
3:
4:
5:

With probability 1 − τ̃j /τ̃jold set cj := 0 and set τ̃jold := 1.

Otherwise set τ̃jold := τ̃j .
end if
if cj = 0 then
⊲ Sample from new columns in D
for ℓ := 1, . . . , t do
With probability τ̃ℓ c(k log k+kt log(1/δ)/ǫ) set cj := dℓ and set τ̃jold := τ̃ℓ
end for
end if
end for
count := 0
end if
end for
function ApproximateRidgeScores(B, M ∈ Rn×t , f robA)
for i := t + 1, . . . , d do
+

f robA−kBk k2F
τ̃k := 4mTi BBT +
I
mi
k
end for
return [τ̃1 , ..., τ̃t ]
24

6:

end function

To prove the correctness of Algorithm 2, we first note that, if our output C has columns
belonging to the claimed distribution, then with probability (1 − δ), C is a (1 + ǫ) error column
subset for A satisfying the guarantees of Theorem 7. Our procedure is not quite equivalent to
the sampling procedure from
a
Pn Theorem 7 because we have some positive probability of choosing
1
0 column (in fact, since j=1 τ̃j ≤ 16k, by our choice of t that probability is greater than 2 for
each column). However, Algorithm 2 samples from a distribution that is equivalent to sampling
from A with an all zeros column 0 tacked on and assigned a high ridge leverage score overestimate.
Furthermore, by inspecting Algorithm 2, we can see that each column is sampled independently, as
all ridge leverage score estimates are computed using the deterministic sketch B. Thus, we obtain
a column subset for [A ∪ 0], which is clearly also a column subset for A.
So, we just need to argue that we obtain an output according to the claimed distribution.
Consider the state of the algorithm after each set of t “Process column stream” iterations, or
equivalently, after each time the “Prune columns” else statement is entered. Denote A’s first t
columns as A(1) , its first 2t columns as A(2) , and in general, its first m · t columns as A(m) . These
submatrices represent the columns of A processed by the end of each epoch of t “Process column
stream” iterations. Let’s take as an inductive assumption that after every prior set of t steps, each
column in C equals:

c=

(

aj ∈ [A(m) ] with probability pj ∈
P
0 with probability (1 − j pj ),

h

1 τ̃j c(k log k+k log(1/δ)/ǫ) τ̃j c(k log k+k log(1/δ)/ǫ)
,
2
t
t

i

,

(37)

P
where τ̃j ≥ 2τ̄j (A(m) ) for all j and j τ̃j ≤ 16k. This is simply equivalent to our claimed output
property of C once all columns have been processed.
(37) holds for A(1) because all of its columns are initially stored in the buffer D and each
τ̃ c(k log k+k log(1/δ)/ǫ)
(see line 19). From Lemma 17 and our
c is set to dj with probability pj = j
t
chosen scaling by 4 (line 3 of ApproximateRidgeScores),
we know that τ̃j ≥ 2τ̄j (A(1) ). Additionally,
P
(1)
τ̃j ≤ 8τ̄j (A ), so it follows from Lemma 4 that j τ̃j ≤ 16k.
For future iterations, A(m) equals [A(m−1) , D]. Consider the columns in A(m−1) first. By our
(m−1) with probability
inductive
 assumption each column in C has already
 been set to aj ∈ A
pj ∈

old c(k log k+k log(1/δ)/ǫ)

1 τ̃j
2

t

,

τ̃jold c(k log k+k log(1/δ)/ǫ)
t

. Our “Rejection sample” step additionally

filters out any column sampled with probability τ̃j /τ̃jold , meaning that in total aj is sampled with
the desired probability from (37). We note that τ̃j /τ̃jold is trivially ≤ 1 since τ̃j was set to the
minimum of τ̃jold and the ridge leverage score of aj computed with respect to our updated Frequent
Directions sketch (see line 10).
If it was set based on the updated Frequent Directions sketch, then the argument that τ̃j ≥
2τ̄j (A(m) ) is the same as for A(1) . On the other hand, if τ̃j was set to equal τ̃jold , then we can apply
Lemma 9: from the inductive assumption, τ̃j = τ̃jold ≥ 2τ̄j (A(m−1) ) and τ̄j (A(m−1) ) ≥ τ̄j (A(m) )
from the monotonicity property so τ̃j ≥ 2τ̄j (A(m) ).
Next consider any aj ∈ D. Each column c is set to aj with the correct probability for (37),
but only conditioned on the fact that c = 0 before the “Sample from new columns in D” if
statement is reached. This conditioning should mean that we effectively sample each aj ∈ D
with lower probability. However,
P the probability cannot be much lower: by our choice of t and
the inductive assumption on j τ̃j , every column is set to 0 with at minimum 1/2 probability.
25

Accordingly, c is available at least half the time, meaning that we at least sample aj with probability
τ̃ c(k log k+k log(1/δ)/ǫ)
, which satisfies (37).
pj = 12 j
t
P
All that is left to argue is that j τ̃j ≤ 16k for A(m) . The argument is the same as for A(1) , the
only difference being that for some values of j, we could have set τ̃j = τ̃jold , which only decreases
the total sum. We conclude by induction that (37) holds for A itself, and thus C is a (1 + ǫ) error
column subset (Theorem 7). Algorithm 2 requires O(nk) space to store B and maintains at most
t = O(k log k + k log(1/δ)/ǫ) sampled columns. It thus proves Theorem 18:
Theorem 18 (Streaming Column Subset Selection). There exists a streaming algorithm that uses
just a single-pass over A’s columns to compute a (1 + ǫ) error column subset C with O(k log k +
k log(1/δ)/ǫ) columns. The algorithm uses O(nk) space in addition to the space required to store
C and succeeds with probability 1 − δ.
We note that, by using the stronger containment condition of Theorem 7 and the streaming
projection-cost preserving sampling algorithm described below we can easily modify the above
algorithm to output an optimally sized column subset with O(k/ǫ) columns. In order to select this
subset, we require a Frequent Directions sketch with ǫ error, so that we can evaluate each O(k/ǫ)
sized subset in our set of O(k log(1/δ)/ǫ) ‘adaptively sampled’ columns and return one giving ǫ
error. The higher accuracy Frequent Directions sketch incurs space overhead O(nk/ǫ).

6.4

Streaming Projection-Cost Preserving Samples

Our single-pass streaming procedure for projection-cost preserving samples is similar to Algorithm
2, although with one important difference. When constructing column subsets, we sampled new
columns in the buffer D while ignoring the fact that “available slots” in C (i.e. columns currently
set to 0) had already been consumed with some probability. This decision was deliberate, rather
than a convenience for analysis. We could not account for the probability of slots being unavailable
because calculating that probability precisely would require knowing the ridge leverage scores of
already discarded columns.
Fortunately, the probability of a column not being set to 0 was bounded by 1/2 and our
procedure hits its sampling target up to this factor. However, while a constant factor approximation
to sampling probabilities is also sufficient for our Theorem 6 projection-cost preservation result,
the fact that columns need to be reweighted by the inverse of their sampling probability adds a
complication: we do not know the true probability with which we sampled each column!
Unfortunately, approximating the reweighting
up to a constant factor is insufficient. We need
p
√
to reweight columns by a factor within (1 ± ǫ) of 1/ tpi for Theorem 5 and Lemma 20 to hold
(which are both required for Theorem 6). This is easily checked by noting that such a reweighting
is equivalent to replacing CCT with CWCT where (1 − ǫ)Id×d  W  (1 + ǫ)Id×d .
We achieve this accuracy by modifying our algorithm so that it maintains an even higher “open
rate” within C. Specifically, we choose t so that each column c has at least a (1 − ǫ) probability of
equaling 0 at any given point in our stream. The procedure is given as Algorithm 3. The constant
c is the required oversampling parameter from Theorem 7.
Algorithm 3 Streaming Projection-Cost Preserving Samples
input: A ∈ Rn×d , accuracy ǫ, success probability (1 − δ)
1
2 and each column c is equal to column
ck log(k/δ)/ǫ
output: C ∈ Rn×t such that t = 16ǫ
i
h
i
τ̃
ck
log(k/δ)/ǫ2 τ̃j ck log(k/δ)/ǫ2
j
1
√
a
with
probability
p
∈
(1
−
ǫ)
,
and 0 otherwise,
j
j
t
t
τ̃j ck log(k/δ)/ǫ2
Pn
where τ̃j ≥ 2τ̄j (A) for all j and j=1 τ̃j ≤ 16k.
26

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
1:
2:
3:
4:
5:
6:

count := 1, C := 0n×t , D := 0n×t , f robA := 0
⊲ Initialize storage
old
old
⊲ Initialize sampling probabilities
[τ̃1 , ..., τ̃t ] := 1
for i := 1, . . . , d do
⊲ Process column stream
B := FreqDirUpdate(B, ai )
if count ≤ t then
⊲ Collect t new columns
dcount := ai .
⊲ Update kAk2F
f robA := f robA + kai k22
count := count + 1
else
⊲ Prune columns
 old
old
[τ̃1 , ..., τ̃t ] := min [τ̃1 , ..., τ̃t ], ApproximateRidgeScores(B, C, f robA)
[τ̃1D , ..., τ̃tD ] := ApproximateRidgeScores(B, D, f robA)
for j := 1, . . . , t do
if cj 6= 0 then
⊲ Rejection sample


With probability 1 − τ̃j /τ̃jold set cj := 0 and set τ̃jold := 1.
q
Otherwise set τ̃jold := τ̃j and multiply cj by τ̃jold /τ̃j .
end if
if cj = 0 then
⊲ Sample from new columns in D
for ℓ := 1, . . . , t do
2
1
set cj := √
d and set τ̃jold := τ̃ℓ
With probability τ̃ℓ ck log(k/δ)/ǫ
t
2 ℓ
τ̃ℓ ck log(k/δ)/ǫ

end for
end if
end for
count := 0
end if
end for
function ApproximateRidgeScores(B, M ∈ Rn×t , f robA)
for i := t + 1,. . . , d do
+
f robA−kBk k2F
I
mi
τ̃k := 4mTi BBT +
k
end for
return [τ̃1 , ..., τ̃t ]
end function

The analysis of Algorithm 3 is equivalent to that of Algorithm 2, along with the additional
observation that our true sampling probability, pj , is within an ǫ factor of the sampling probability
τ̃ ck log(k/δ)/ǫ2

used for reweighting, j
. Note that while C contains just O(k log(k/δ)/ǫ2 ) non-zero
t
columns in expectation, during the course of a the column stream it could contain as many as
O(k log(k/δ)/ǫ3 ) columns. Regardless, it is always possible to resample from C after running
Algorithm 3 to construct an optimally sized sample for A with error (1 + 2ǫ). Overall we have:
Theorem 19 (Streaming Projection-Cost Preserving Sampling). There exists a streaming algorithm that uses just a single-pass over A’s columns to compute a (1 + ǫ) error projection-cost
preserving sample C with O(k log(k/δ)/ǫ2 ) columns. The algorithm requires a fixed O(nk) space
overhead along with space to store O(k log(k/δ)/ǫ3 ) columns of A. It succeeds with probability 1− δ.

27

References
[AHPV04] Pankaj K. Agarwal, Sariel Har-Peled, and Kasturi R. Varadarajan. Approximating
extent measures of points. J. ACM, 51(4):606–635, 2004.
[AM15]

Ahmed El Alaoui and Michael W. Mahoney. Fast randomized kernel methods with statistical guarantees. In Advances in Neural Information Processing Systems 28 (NIPS),
2015. Full version at arXiv:1411.0306v1.

[AMS01]

Dimitris Achlioptas, Frank Mcsherry, and Bernhard Schölkopf. Sampling techniques
for kernel methods. In Advances in Neural Information Processing Systems 14 (NIPS),
2001.

[Bac13]

Francis Bach. Sharp analysis of low-rank kernel matrix approximations. In Proceedings
of the 26th Annual Conference on Computational Learning Theory (COLT), 2013.

[BDN15]

Jean Bourgain, Sjoerd Dirksen, and Jelani Nelson. Toward a unified theory of sparse dimensionality reduction in euclidean space. Geometric and Functional Analysis (GAFA),
25(4):1009–1088, 2015. Preliminary version in the 47th Annual ACM Symposium on
Theory of Computing (STOC).

[BJS15]

Srinadh Bhojanapalli, Prateek Jain, and Sujay Sanghavi. Tighter low-rank approximation via sampling the leveraged element. In Proceedings of the 26th Annual ACM-SIAM
Symposium on Discrete Algorithms (SODA), 2015.

[BK13]

Wei Bi and James Kwok. Efficient multi-label classification with many labels. In
Proceedings of the 30th International Conference on Machine Learning (ICML), pages
405–413, 2013.

[BS80]

Jon Louis Bentley and James B. Saxe. Decomposable searching problems I: Static-todynamic transformation. Journal of Algorithms, 1(4):301–358, 1980.

[BW09a]

Mohamed-Ali Belabbas and Patrick J. Wolfe. On landmark selection and sampling
in high-dimensional data analysis. Philosophical Transactions of the Royal Society A,
367:4295–4312, 2009.

[BW09b]

Mohamed-Ali Belabbas and Patrick J. Wolfe. Spectral methods in machine learning:
New strategies for very large datasets. Proceedings of the National Academy of Sciences
of the USA, 106:369–374, 2009.

[BW14]

Christos Boutsidis and David P. Woodruff. Optimal CUR matrix decompositions. In
Proceedings of the 46th Annual ACM Symposium on Theory of Computing (STOC),
pages 353–362, 2014.

[BZMD15] Christos Boutsidis, Anastasios Zouzias, Michael W. Mahoney, and Petros Drineas.
Randomized dimensionality reduction for k-means clustering. IEEE Transactions on
Information Theory, 61(2):1045–1062, Feb 2015.
[CBSW15] Yudong Chen, Srinadh Bhojanapalli, Sujay Sanghavi, and Rachel Ward. Coherent
matrix completion. Journal of Machine Learning Research, 2015. Preliminary version
in the 31st International Conference on Machine Learning (ICML).

28

[CEM+ 15] Michael B. Cohen, Sam Elder, Cameron Musco, Christopher Musco, and Madalina
Persu. Dimensionality reduction for k-means clustering and low rank approximation.
In Proceedings of the 47th Annual ACM Symposium on Theory of Computing (STOC),
pages 163–172, 2015.
[CLM+ 15] Michael B. Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng,
and Aaron Sidford. Uniform sampling for matrix approximation. In Proceedings of the
6th Conference on Innovations in Theoretical Computer Science (ITCS), pages 181–
190, 2015.
[CMP15]

Michael B. Cohen, Cameron Musco, and Jakub Pachocki. Online row sampling. In Proceedings of the 19th International Workshop on Approximation Algorithms for Combinatorial Optimization Problems (APPROX), 2015.

[Coh16]

Michael B. Cohen. Nearly tight oblivious subspace embeddings by trace inequalities.
In Proceedings of the 27th Annual ACM-SIAM Symposium on Discrete Algorithms
(SODA), pages 278–287, 2016.

[CW09]

Kenneth L. Clarkson and David P. Woodruff. Numerical linear algebra in the streaming
model. In Proceedings of the 41st Annual ACM Symposium on Theory of Computing
(STOC), pages 205–214, 2009.

[CW13]

Kenneth L. Clarkson and David P. Woodruff. Low rank approximation and regression
in input sparsity time. In Proceedings of the 45th Annual ACM Symposium on Theory
of Computing (STOC), pages 81–90, 2013.

[DFK+ 04] Petros Drineas, Alan Frieze, Ravi Kannan, Santosh Vempala, and V Vinay. Clustering
large graphs via the singular value decomposition. Machine Learning, 56(1-3):9–33,
2004. Preliminary version in the 10th Annual ACM-SIAM Symposium on Discrete
Algorithms (SODA).
[DKM06a] Petros Drineas, Ravi Kannan, and Michael W. Mahoney. Fast Monte Carlo algorithms
for matrices II: Computing a low-rank approximation to a matrix. SIAM J. Comput.,
36(1):158–183, 2006.
[DKM06b] Petros Drineas, Ravi Kannan, and Michael W. Mahoney. Fast Monte Carlo algorithms
for matrices III: Computing a compressed approximate matrix decomposition. SIAM
J. Comput., 36(1):184–206, 2006. Preliminary version in the 14th Annual ACM-SIAM
Symposium on Discrete Algorithms (SODA).
[DMM06]

Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan. Sampling algorithms
for ℓ2 regression and applications. In Proceedings of the 17th Annual ACM-SIAM
Symposium on Discrete Algorithms (SODA), pages 1127–1136, 2006.

[DMM08]

Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan. Relative-error CUR
matrix decompositions. SIAM Journal on Matrix Analysis and Applications, 30(2):844–
881, 2008.

[DRVW06] Amit Deshpande, Luis Rademacher, Santosh Vempala, and Grant Wang. Matrix approximation and projective clustering via volume sampling. Theory of Computing,
2(1):225–247, 2006. Preliminary version in the 17th Annual ACM-SIAM Symposium
on Discrete Algorithms (SODA).
29

[DV06]

Amit Deshpande and Santosh Vempala. Adaptive sampling and fast low-rank matrix
approximation. In Proceedings of the 10th International Workshop on Randomization
and Computation (RANDOM), pages 292–303, 2006.

[FKV04]

Alan Frieze, Ravi Kannan, and Santosh Vempala. Fast Monte-Carlo algorithms for
finding low-rank approximations. J. ACM, 51(6):1025–1041, November 2004. Preliminary version in the 39th Annual IEEE Symposium on Foundations of Computer Science
(FOCS).

[FS02]

Shai Fine and Katya Scheinberg. Efficient SVM training using low-rank kernel representations. The Journal of Machine Learning Research, 2:243–264, 2002.

[FSS13]

Dan Feldman, Melanie Schmidt, and Christian Sohler. Turning big data into tiny data:
Constant-size coresets for k-means, PCA, and projective clustering. In Proceedings
of the 24th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages
1434–1453, 2013.

[GLPW15] Mina Ghashami, Edo Liberty, Jeff M. Phillips, and David P. Woodruff. Frequent
Directions: Simple and deterministic matrix sketching. arXiv:1501.01711, 2015.
[GM13]

Alex Gittens and Michael Mahoney. Revisiting the Nyström method for improved
large-scale machine learning. In Proceedings of the 30th International Conference on
Machine Learning (ICML), pages 567–575, 2013.

[HI15]

John T. Holodnak and Ilse C. F. Ipsen. Randomized approximation of the Gram matrix:
Exact computation and probabilistic bounds. SIAM Journal on Matrix Analysis and
Applications, 36(1):110–137, 2015.

[HKZ12]

Daniel Hsu, Sham Kakade, and Tong Zhang. Tail inequalities for sums of random
matrices that depend on the intrinsic dimension. Electron. Commun. Probab., 17:1–13,
2012.

[HPM04]

Sariel Har-Peled and Soham Mazumdar. On coresets for k-means and k-median clustering. In Proceedings of the 36th Annual ACM Symposium on Theory of Computing
(STOC), pages 291–300, 2004.

[HRZ+ 09] David Hall, Daniel Ramage, Jason Zaugg, Alexander Lehmann, Jonathan Merritt,
Keith Stevens, Jason Baldridge, Timothy Hunter, Dave DeCaprio, Daniel Duckworth,
Eric Christiansen, Marc Millstone, Mérő László, Alexey Noskov, Devon Bryant, Kentaroh Takagaki, Sam Halliday, Chris Stucchio, and Xiangrui Meng. ScalaNLP: Breeze.
http://www.scalanlp.org/, 2009.
[IBM14]

IBM Reseach Division, Skylark Team. libskylark: Sketching-based Distributed Matrix
Computations for Machine Learning. IBM Corporation, Armonk, NY, 2014.

[JK16]

Gorav Jindal and Pavel Kolev. An efficient parallel algorithm for spectral sparsification
of laplacian and SDDM matrix polynomials. arXiv:1507.07497, 2016.

[KL13]

Jonathan A. Kelner and Alex Levin. Spectral sparsification in the semi-streaming
setting. Theory of Computing Systems, 53(2):243–262, 2013. Preliminary version in the
28th International Symposium on Theoretical Aspects of Computer Science (STACS).

30

[KLM+ 14] Michael Kapralov, Yin Tat Lee, Cameron Musco, Christopher Musco, and Aaron Sidford. Single pass spectral sparsification in dynamic streams. In Proceedings of the
55th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages
561–570, 2014.
[KLP+ 16] Rasmus Kyng, Yin Tat Lee, Richard Peng, Sushant Sachdeva, and Daniel A. Spielman.
Sparsified cholesky and multigrid solvers for connection laplacians. In Proceedings of
the 48th Annual ACM Symposium on Theory of Computing (STOC), pages 842–850,
2016.
[Lib13]

Edo Liberty. Simple and deterministic matrix sketching. In Proceedings of the 19th
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining
(KDD), pages 581–588, 2013.

[Liu14]

Antoine Liutkus. Randomized SVD. http://www.mathworks.com/matlabcentral/fileexchange/478
2014. MATLAB Central File Exchange.

[LJS16]

Chengtao Li, Stefanie Jegelka, and Suvrit Sra. Fast DPP sampling for nyström with
application to kernel methods. In Proceedings of the 33rd International Conference on
Machine Learning (ICML), 2016.

[LMP13]

Mu Li, Gary L. Miller, and Richard Peng. Iterative row sampling. In Proceedings of
the 54th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages
127–136, 2013. Preliminary version at arXiv:1211.2713v1.

[LPS15]

Yin Tat Lee, Richard Peng, and Daniel A. Spielman. Sparsified cholesky solvers for
SDD linear systems. arXiv:1506.08204, 2015.

[LS14]

Yin Tat Lee and Aaron√Sidford. Path finding methods for linear programming: Solving
linear programs in Õ( rank) iterations and faster algorithms for maximum flow. In
Proceedings of the 55th Annual IEEE Symposium on Foundations of Computer Science
(FOCS), pages 424–433, 2014.

[LS15]

Yin Tat Lee and Aaron Sidford. Efficient inverse maintenance and faster algorithms
for linear programming. In Proceedings of the 56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), 2015.

[LSW15]

Yin Tat Lee, Aaron Sidford, and Sam Chiu-wai Wong. A faster cutting plane method
and its implications for combinatorial and convex optimization. In Proceedings of the
56th Annual IEEE Symposium on Foundations of Computer Science (FOCS), 2015.

[MA77]

George A. Milliken and Fikri Akdeniz. A theorem on the difference of the generalized inverses of two nonnegative matrices. Communications in Statistics - Theory and
Methods, 6(1):73–79, 1977.

[MD05]

Michael W. Mahoney and Petros Drineas. On the Nyström method for approximating a gram matrix for improved kernel-based learning. Journal of Machine Learning
Research, 6:2153–2175, 2005. Preliminary version in the 18th Annual Conference on
Computational Learning Theory (COLT).

[MD09]

Michael W. Mahoney and Petros Drineas. CUR matrix decompositions for improved
data analysis. Proceedings of the National Academy of Sciences of the USA, 106(3):697–
702, 2009.
31

[Min13]

Stanislav Minsker. On some extensions of Bernstein’s inequality for self-adjoint operators. arXiv:1112.5448, 2013.

[MM13]

Michael W. Mahoney and Xiangrui Meng. Low-distortion subspace embeddings in
input-sparsity time and applications to robust linear regression. In Proceedings of the
45th Annual ACM Symposium on Theory of Computing (STOC), pages 91–100, 2013.

[MM15]

Cameron Musco and Christopher Musco. Randomized block krylov methods for
stronger and faster approximate singular value decomposition. In Advances in Neural Information Processing Systems 28 (NIPS), 2015. Full version at arXiv:1504.05477.

[MM16]

Cameron Musco and Christopher Musco. Provably useful kernel matrix approximation
in linear time. arXiv:1605.07583, 2016.

[Mus15]

Cameron Musco. Dimensionality reduction for k-means clustering. Master’s thesis,
Massachusetts Institute of Technology, 2015.

[NN13]

Jelani Nelson and Huy L. Nguyen. OSNAP: Faster numerical linear algebra algorithms
via sparser subspace embeddings. In Proceedings of the 54th Annual IEEE Symposium
on Foundations of Computer Science (FOCS), pages 117–126, 2013.

[Oka10]

Daisuke
Okanohara.
redsvd:
https://code.google.com/p/redsvd/, 2010.

RandomizED

SVD.

[PVG+ 11] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel,
P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau,
M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python.
Journal of Machine Learning Research, 12:2825–2830, 2011.
[PZB+ 07]

Peristera Paschou, Elad Ziv, Esteban G. Burchard, Shweta Choudhry, William
Rodriguez-Cintron, Michael W. Mahoney, and Petros Drineas. PCA-correlated SNPs
for structure identification in worldwide human populations. PLoS Genet, 3(9):1672–
1686, 2007.

[RR07]

Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In
Advances in Neural Information Processing Systems 20 (NIPS), pages 1177–1184, 2007.

[Sar06]

Támas Sarlós. Improved approximation algorithms for large matrices via random projections. In Proceedings of the 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS), pages 143–152, 2006.

[SS11]

Daniel A. Spielman and Nikhil Srivastava. Graph sparsification by effective resistances.
SIAM Journal on Computing, 40(6):1913–1926, 2011. Preliminary version in the 40th
Annual ACM Symposium on Theory of Computing (STOC).

[Str14]

Martin Tobias Strauch. Column subset selection with applications to neuroimaging
data. PhD thesis, Universität Konstanz, 2014.

[Tro15]

Joel A. Tropp. An introduction to matrix concentration inequalities. Foundations and Trends in Machine Learning, 8(1-2):1–230, 2015. Preliminary version at
arXiv:1501.01571.

32

[VM15]

Sergey Voronin and Per-Gunnar Martinsson. RSVDPACK: Subroutines for computing
partial singular value decompositions via randomized sampling on single core, multi
core, and GPU architectures. arXiv:1502.05366, 2015.

[WS01]

Christopher K. I. Williams and Matthias Seeger. Using the Nyström method to speed
up kernel machines. In Advances in Neural Information Processing Systems 13 (NIPS),
pages 682–688, 2001.

A

Trace Bound for Ridge Leverage Score Sampling

Lemma 20. For i ∈ {1, . . . , d}, let τ̃i ≥ τ̄i (A) be an overestimate for the ith ridge leverage score.
P
· i τ̃i , for some sufficiently large constant c. Construct C by
Let pi = Pτ̃i τ̃i . Let t = c log(k/δ)
ǫ2
i

sampling t columns of A, each set to

2 ≥
singular value with σm

kA\k k2F
k

√1 ai
tpi

with probability pi . Let m be the index of the smallest

. With probability 1 − δ, C satisfies:

| tr(A\m AT\m ) − tr(U\m UT\m CCT U\m UT\m )| ≤ ǫkA\k k2F .

(38)

Proof. Letting P\m = U\m UT\m , we can rewrite (38) as:
|kP\m Ck2F − kA\m k2F | ≤ ǫkA\k k2F .
We can write kP\m Ck2F as a sum over column norms:
kP\m Ck2F =

t
X
j=1

kP\m cj k22 .

Now, for some i ∈ {1, . . . , d} and recalling our definition Σ̄2i,i = σi2 (A) +
kP\m ci k22 =

kA\k k2F
k

we have:

kP\m ai k22
1
ǫ2
·
kP\m ai k22 ≤
tpi
c log(k/δ)
τ̄i (A)
2
kP\m ai k22
ǫ

· T
=
c log(k/δ) ai UΣ̄−2 UT ai

=

kP\m ai k22
ǫ2

· T
c log(k/δ) ai P\m UΣ̄−2 UT P\m ai

≤

ǫ2
·
c log(k/δ)

≤

2ǫ2
· kA\k k2F ,
c log(k/δ)

≤

kP\m ai k22
ǫ2


·
c log(k/δ) aT P
−2 T
\m U\m Σ̄ U\m P\m ai
i
kP\m ai k22

k
2kA\k k2F

kP\m ai k22

where the second to last inequality follows from the fact that Σ̄2i,i = σi2 (A) +
i ≥ m. Therefore, U\m Σ̄−2 UT\m  2kAk k2 P\m .
\k F

33

kA\k k2F
k

≤

2kA\k k2F
k

for

So,
bound:

c log(k/δ)
2ǫ2 kA\k k2F

· kP\m ci k22 ∈ [0, 1]. We have E

hP

t
2
j=1 kP\m ci k2

i

= kA\m k2F so by a Chernoff



P kP\m Ck2F ≥ kA\m k2F + ǫkA\k k2F


!
t
2
2
X
c
log(k/δ)kA
k
ǫkA
k
c log(k/δ)
\m F 
\k F
kP\m ci k22 ≥ 1 +
= P 2
2ǫ kA\k k2F j=1
kA\m k2F
2ǫ2 kA\k k2F
≤ e−c log(k/δ)/4 ≤ δ/2,

kA

k2

if we set c sufficiently large. In the second to last step we use the fact that kA \k kF2 ≥ 21 by the
\m F


definition of m. We can similarly prove that P kP\m Ck2F ≤ kA\m k2F − ǫkA\k k2F ≤ δ/2. Union
bounding gives the result.

B

Independent Sampling Bounds

In this section we give analogies to Theorem 5 and Lemma 20 when columns are sampled independently using their ridge leverage scores rather than sampled with replacement.
o
n
,
1
for
Lemma 21. For i ∈ {1, . . . , d}, given τ̃i ≥ τ̄i (A) for all i, let pi = min τ̃i · c log(k/δ)
2
ǫ
some sufficiently large constant c. Construct C by independently sampling each column ai from
√
A with probability
pi and scaling selected columns by 1/ pi . With probability 1 − δ, C has
P
O log(k/δ)/ǫ2 · i τ̃i columns and satisfies:
ǫ
ǫ
(1 − ǫ)CCT − kA − Ak k2F In×n  AAT  (1 + ǫ)CCT + kA − Ak k2F In×n .
k
k

(7)

Proof. Again we rewrite the ridge leverage score definition using A’s singular value decomposition:
!+
2
kA
k
\k
F
τ̄i (A) = aTi UΣ2 UT +
UUT
ai
k
+

= aTi UΣ̄2 UT ai = aTi UΣ̄−2 UT ai ,
where Σ̄2i,i = σi2 (A) +

kA\k k2F
k

Xi =
Let Y =

P

. For each i ∈ 1, . . . , d define the matrix valued random variable:

(

1
pi


− 1 Σ̄−1 UT ai aTi UΣ̄−1 with probability pi

−Σ̄−1 UT ai aTi UΣ̄−1 with probability (1 − pi )

Xi . We have E Y = 0. Furthermore, CCT = UΣ̄Y Σ̄U + AAT . Showing kYk2 ≤ ǫ

gives −ǫI  Y  ǫI, and since UΣ̄2 UT = AAT +
(1 − ǫ)AAT −

kA\k k2F
k

I would give:

ǫkA\k k2F
ǫkA\k k2F
I  CCT  (1 + ǫ)AAT +
I.
k
k

After rearranging and adjusting constants on ǫ, this statement is equivalent to (7).

34

To prove that kYk2 is we use the same stable rank matrix Bernstein inequality used for our
with replacement results [Tro15]. If pi = 1 (i.e. τ̃i · c log(k/δ)/ǫ2 ≥ 1) then Xi = 0 so kXi k2 = 0.
1
T
τ̄i (A) ai ai

Otherwise, we use the fact that

 AAT +

kA\k k2F
k

I, which lets us bound:

!
2
kA
k
\k F
AAT +
I UΣ̄−1 = I.
k

1
· Σ̄−1 UT ai aTi UΣ̄−1  Σ̄−1 UT
τ̃i
2

2

ǫ
ǫ
I and hence kXi k2 ≤ c log(k/δ)
.
So we have Xi  p1i Σ̄−1 UT ai aTi UΣ̄−1  c log(k/δ)
Next we bound the variance of Y.
#
" 
2
X
X
1
2
2
E(Y ) =
E(Xi ) 
pi
− 1 + (1 − pi ) · Σ̄−1 UT ai aTi UΣ̄−2 UT ai aTi UΣ̄−1
pi
X 1
ǫ2

Σ̄−1 UT AAT UΣ̄−1
· τ̄i (A) · Σ̄−1 UT ai aTi UΣ̄−1 
pi
c log(k/δ)
ǫ2
ǫ2
Σ2 Σ̄−2 
D.
(39)

c log(k/δ)
c log(k/δ)

where again we set Di,i = 1 for i ∈ 1, . . . , k and Di,i =

σi2
σi2 +kA\k k2F /k

for all i ∈ k + 1, . . . , n. By the

stable rank matrix Bernstein inequality given in Theorem 7.3.1 of [Tro15], for ǫ < 1,
4 tr(D)
e
P [kYk ≥ ǫ] ≤
kDk2



−ǫ2 /2

ǫ2
(kDk2 +ǫ/3)
c log(k/δ)

.

(40)

Clearly kDk2 = 1. Furthermore,
tr(D) = k +

d
X

i=k+1

σi2
σi2 +

kA\k k2F
k

≤k+

d
X

i=k+1

σi2
kA\k k2F

=k+

k

Pd

2
i=k+1 σi
2
kA\k kF
k

≤ k + k.

Plugging into (9), we see that
P [kYk ≥ ǫ] ≤ 8ke−

c log(k/δ)
2

) ≤ δ/2,

if we choose the constant c large enough. So we have established (7).
All that remains to note is that, the expected number of columns in C is at most c log(k/δ)
·
ǫ2


Pd
log(k/δ) P
· i τ̃i columns with probability > 1 − δ/2 by a
i=1 τ̃i . Accordingly, C has at most O
ǫ2
standard Chernoff bound. Union bounding over failure probabilities gives the lemma.
o
n
,
1
for
Lemma 22. For i ∈ {1, . . . , d}, given τ̃i ≥ τ̄i (A) for all i, let pi = min τ̄i (A) · c log(k/δ)
2
ǫ
some sufficiently large constant c. Construct C by independently sampling each column ai from
√
A with probability pi and scaling selected columns by 1/ pi . Let m be the index of the smallest
2 ≥
singular value with σm

kA\k k2F
k

. With probability 1 − δ, C satisfies:

| tr(A\m AT\m ) − tr(U\m UT\m CCT U\m UT\m )| ≤ ǫkA\k k2F .

35

(41)

Proof. We need to show tr(A\m AT\m ) − tr(U\m UT\m BBT U\m UT\m ) ≥ −ǫkA\m k2F . Letting P\m =
U\m UT\m , we can rewrite this as:
kP\m Bk2F − kA\m k2F ≤ ǫkA\m k2F .
We can write kP\m Bk2F as a sum over column norms:
kP\m Bk2F =

d
X
i=1

Ii

1
kP ai k2 ,
pi \m 2

where Ii is an indicator random variable equal to 1 with probability pi and 0 otherwise.
We have:
kP\m ai k22
kP\m ai k22
ǫ2
ǫ2
1
2

kP ai k =
·
≤
·
pi \m 2 c log(k/δ)
τ̃i
c log(k/δ) aTi UΣ̄−2 UT ai

=

kP\m ai k22
ǫ2

· T
c log(k/δ) ai P\m UΣ̄−2 UT P\m ai

≤

ǫ2
·
c log(k/δ)

≤

2ǫ2
· kA\k k2F ,
c log(k/δ)

≤

kP\m ai k22
ǫ2


·
c log(k/δ) aT P
−2 UT
U
Σ̄
\m
\m
i
\m P\m ai
kP\m ai k22

k
2kA\k k2F

kP\m ai k22

where the second to last inequality follows from the fact that Σ̄2i,i = σi2 (A) +
i ≥ m. Therefore, U\m Σ̄−2 UT\m  2kAk k2 P\m .
So

P

c log(k/δ)
2ǫ2 kA\m k2F



kA\k k2F
k

≤

2kA\k k2F
k

\k F

·

kP\m Bk2F

1
2
pi kP\m ai k2

≥ (1 +

∈ [0, 1] and by a Chernoff bound we have:

ǫ)kA\m k2F



d
c log(k/δ)
c log(k/δ) X 1
=P
Ii kP\m ai k22 ≥ (1 + ǫ)
2
2
2ǫ2
2ǫ kA\m kF i=1 pi

"

≤ e−c log(k/δ)/4 ≤ δ/2,

if we set c sufficiently large.

36

#

for

