Online Row Sampling

arXiv:1604.05448v1 [] 19 Apr 2016

Michael B. Cohen
MIT∗
micohen@mit.edu

Cameron Musco
MIT†
cnmusco@mit.edu

Jakub Pachocki
Carnegie Mellon University‡
pachocki@cs.cmu.edu

April 20, 2016

Abstract
Finding a small spectral approximation for a tall n × d matrix A is a fundamental numerical
primitive. For a number of reasons, one often seeks an approximation whose rows are sampled
from those of A. Row sampling improves interpretability, saves space when A is sparse, and
preserves row structure, which is especially important, for example, when A represents a graph.
However, correctly sampling rows from A can be costly when the matrix is large and cannot
be stored and processed in memory. Hence, a number of recent publications focus on row
sampling in the streaming setting, using little more space than what is required to store the
outputted approximation [KL13, KLM+ 14].
Inspired by a growing body of work on online algorithms for machine learning and data analysis, we extend this work to a more restrictive online setting: we read rows of A one by one and
immediately decide whether each row should be kept in the spectral approximation or discarded,
without ever retracting these decisions. We present an extremely simple algorithm that approximates A up to multiplicative error ǫ and additive error δ using O(d log d log(ǫkAk22 /δ)/ǫ2 ) online
samples, with memory overhead proportional to the cost of storing the spectral approximation.
We also present an algorithm that uses O(d2 ) memory but only requires O(d log(ǫkAk22 /δ)/ǫ2 )
samples, which we show is optimal.
Our methods are clean and intuitive, allow for lower memory usage than prior work, and
expose new theoretical properties of leverage score based matrix approximation.

∗

Work supported by NSF grant CCF-1111109.
Work partially supported by NSF Graduate Research Fellowship No. 1122374, AFOSR grant FA9550-13-1-0042
and the NSF Center for Science of Information.
‡
Work supported by NSF grant CCF-1065106.
†

1
1.1

Introduction
Background

A spectral approximation to a tall n×d matrix A is a smaller, typically Õ(d)×d matrix Ã such that
kÃxk2 ≈ kAxk2 for all x. Typically one asks for a multiplicative approximation, which guarantees
that (1 − ǫ)kAxk22 ≤ kÃxk22 ≤ (1 + ǫ)kAxk22 . In other notation, (1 − ǫ)A  Ã  (1 + ǫ)A.
Such approximations have many applications, most notably for solving least squares regression
over A [CW13, CLM+ 15]. If A is the vertex edge incidence matrix of a graph, Ã is a spectral
sparsifier [ST04]. It can be used to approximate effective resistances, spectral clustering, mixing
time and random walk properties, and many other computations.
A number of recent papers focus on fast algorithms for spectral approximation. Using sparse
random subspace embeddings [CW13, NN13, MM13], it is possible to find Ã in input sparsity time,
essentially by randomly recombining the rows of A into a smaller number of rows. In some cases
these embeddings are not enough, as it is desirable for the rows of Ã to be a subset of rows sampled
from A. If A is sparse, this ensures that Ã is also sparse. If A represents a graph, it ensures that
Ã is also a graph, specifically a weighted subgraph of the original.
It is well known that sampling O(d log d/ǫ2 ) rows of A with probabilities proportional to their
leverage scores yields a (1 + ǫ) multiplicative factor spectral approximation to A. Further, this
sampling can be done in input sparsity time, either using subspace embeddings to approximate
leverage scores, or using iterative sampling techniques [LMP13], some that only work with subsampled versions of the original matrix [CLM+ 15].

1.2

Streaming and Online Row Sampling

When A is very large, input sparsity runtimes are not enough – memory restrictions also become
important. Hence, recent work has tackled row sampling in a streaming model of computation.
[KL13] presents a simple algorithm for sampling rows from an insertion only stream, using space
approximately proportional to the size of the final approximation. [KLM+ 14] gives a sparse-recovery
based algorithm that works in dynamic streams with row insertions and deletions, also using nearly
optimal space. Unfortunately, to handle dynamic streams, the algorithm in [KLM+ 14] is complex,
requires additional restrictions on the input matrix, and uses significantly suboptimal runtime to
recover a spectral approximation from its low memory representation of the input stream.
While the algorithm in [KL13] is simple and efficient, we believe that its proof is incomplete,
and do not see an obvious way to fix it. The main idea behind the algorithm is to sample rows by
their leverage scores with respect to the stream seen so far. These leverage scores may be coarse
overestimates of the true scores. However as more rows are streamed in, better estimates can be
obtained and the sampled rows pruned to a smaller set. Unfortunately, the probability of sampling
a row becomes dependent on which other rows are sampled. This seems to break the argument in
that paper, which essentially claims that their process has the same distribution as would a single
round of leverage score sampling.
In this paper we initiate the study of row sampling in an online setting. As in an insertion
stream, we read rows of A one by one. However, upon seeing a row, we immediately decide
whether it should be kept in the spectral approximation or discarded, without ever retracting these
decisions. We present a similar algorithm to [KL13], however, since we never prune previously
sampled rows, the probability of sampling a row only depends on whether previous rows in the

1

stream were sampled. This limited dependency structure allows us to rigorously argue that a
spectral approximation is obtained.
In addition to addressing gaps in the literature on streaming spectral approximation, our restricted model extends work on online algorithms for a variety of other machine learning and data
analysis problems, including principal component analysis [BGKL15], clustering [LSS16], classification [BB05, CDK+ 06], and regression [CDK+ 06]. In practice, online algorithms are beneficial
since they can be highly computationally and memory efficient. Further, they can be applied in
scenarios in which data is produced in a continuous stream and intermediate results must be output
as the stream is processed. Spectral approximation is a widely applicable primitive for approximate
learning and computation, so studying its implementation in an online setting is a natural direction.

1.3

Our Results

Our primary contribution is a very simple algorithm for leverage score sampling in an online manner.
The main difficultly with row sampling using leverage scores is that leverage scores themselves are
not easy to compute. They are given by li = aTi (AT A)−1 ai , and so require solving systems in AT A
if computed naively. This is not only expensive, but also impossible in an online setting, where we
do not have access to all of A.
A critical observation is that it always suffices to sample rows by overestimates of their true
leverage scores. The number of rows that must be sampled is proportional to the sum of these
overestimates. Since the leverage score of a row can only go up when we remove rows from the
matrix, a simple way to obtain an overestimate is to compute leverage score using just a subset of
the other rows of A. That is, letting Aj contain just j of A’s n rows, we can overestimate li by
˜li = aT (AT Aj )−1 ai
i
j
[CLM+ 15] shows that if Aj is a subset of rows sampled uniformly at random, then the expected
leverage score of ai is d/j. This simple fact immediately gives a result for online sampling from a
randomly ordered stream. If we compute the leverage score of the current row ai against all previously seen rows (or some approximation to these rows), then the expected sum of our overestimates
will be bounded by d + d/2 + ... + ... + d/n = O(d log n). So, sampling O(d log d log n/ǫ2 ) rows will
be enough obtain a (1 + ǫ) multiplicative factor spectral approximation.
What if we cannot guarantee a randomly ordered input stream? Is there any hope of being
able to compute good leverage score estimates in an online manner? Surprisingly the answer to
this is yes - we can in fact run nearly the exact same algorithm and be guaranteed that the sum of
estimated leverage scores is low, regardless of stream order. Roughly, each time we receive a row
which has high leverage score with respect to the previous rows, it must compose a significant part
of A’s spectrum. If A does not continue to grow unboundedly, there simply cannot be too many
of these significant rows.
Specifically, we show that if we sample by the ridge leverage scores [AM] over all previously
seen rows, which are the leverage scores computed over ATi Ai + λI for some small regularizing
factor λ, then with just O(d log d log(ǫkAk22 /δ)/ǫ2 ) samples we obtain a (1 + ǫ) multiplicative, δ
additive error spectral approximation. That is, with high probability we sample a matrix Ã with
(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI.
To gain intuition behind this bound, note that we can convert it into a multiplicative one by
setting δ = ǫσmin (A)2 (as long as we have some estimate of σmin (A)). This setting of δ will require
taking O(d log d log(κ(A))/ǫ2 ) samples. If we have a polynomial bound on the condition number
2

of A, as we do, for instance, for graphs with polynomially bounded edges weights, this becomes
O(d log2 d/ǫ2 ) – nearly matching the O(d log d/ǫ2 ) achievable if sampling by true leverage scores.
Our online sampling algorithm is extremely simple. When each row comes in, we compute
the online ridge leverage score, or an estimate of it, and then irrevocably either add the row to
our approximation or remove it. As mentioned, it is similar in form to the streaming algorithm
of [KL13], except that it does not require pruning previously sampled rows. This allows us to
avoid difficult dependency issues. Additionally, without pruning, we do not even need to store all
previously sampled rows. As long as we store a constant factor spectral approximation our previous
samples, we can compute good approximations to the online ridge leverage scores. In this way, we
can store just O(d log d log(ǫkAk22 /δ)) rows in working memory (O(d log 2 d) if we want a spectral
graph sparsifier), filtering our input stream into an O(d log d log(κ(A))/ǫ2 ) sized output stream.
Note that this memory bound in fact improves as ǫ decreases, and regardless, can be significantly
smaller than the output size of the algorithm.
In additional to our main sampling result, we use our bounds on online ridge leverage score
approximations to show that an algorithm in the style of [BSS12] allows us to remove a log d
factor and sample just O(d log(ǫkAk22 /δ)/ǫ2 ) (Theorem 4.1). This algorithm is more complex and
can require O(d2 ) working memory. However, in Theorem 5.1 we show that it is asymptotically
optimal. The log(ǫkAk22 /δ) factor is not an artifact of our analysis, but is truly the cost of restricting
ourselves to online sampling. No algorithm can obtain a multiplicative (1 + ǫ) additive δ spectral
approximation taking fewer than Ω(d log(ǫkAk22 /δ)/ǫ2 ) rows in an online manner.

2

Overview

Let A be an n × d matrix with rows a1 , . . . , an . A natural approach to row sampling from A is
picking an a priori probability with which each row is kept, and then deciding whether to keep
each row independently. A common choice is for the sampling probabilities to be proportional to
the leverage scores of the rows. The leverage score of the i-th row of A is defined to be
aTi (AT A)† ai ,
where the dagger symbol denotes the pseudoinverse. In this work, we will be interested in approximating AT A with some (very) small multiple of the identity added. Hence, we will be interested
in the λ-ridge leverage scores [AM]:
aTi (AT A + λI)−1 ai ,
for a parameter λ > 0.
In many applications, obtaining the (nearly) exact values of aTi (AT A + λI)−1 ai for sampling
is difficult or outright impossible. A key idea is that as long as we have a sequence l1 , . . . , ln of
overestimates of the λ-ridge leverage scores, that is for i = 1, . . . , n
li ≥ aTi (AT A + λI)−1 ai ,
we can sample by these overestimates and obtain rigorous guarantees on the quality of the obtained
spectral approximation. This notion is formalized in Theorem 2.1.

3

Theorem 2.1 Let A be an n × d matrix with rows a1 , . . . , an . Let ǫ ∈ (0, 1), δ > 0, λ := δ/ǫ, c :=
8 log d/ǫ2 . Assume we are given l1 , . . . , ln such that for all i = 1, . . . , n,
li ≥ aTi (AT A + λI)−1 ai .
For i = 1, . . . , n, let pi := min(cli , 1). Construct Ã by independently sampling each row ai of A with
√
probability pi , and rescaling it by 1/ pi if it is included in the sample. Then, with high probability,
(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI,

P
and the number of rows in Ã is O ( ni=1 li ) log d/ǫ2 .

Proof: This sort of guarantee for leverage score sampling is well known. See
√ for example Lemma
4 of [CLM+ 15]. If we sampled both the rows of A and the rows of λI with the leverage
scores over (AT A + λI), we would have (1 − ǫ)(AT A + λI)  ÃT Ã  (1 + ǫ)(AT A + λI).
However, we do not sample the rows of the identity. Since we could have sampled them each
with probability 1, we can simply subtract λI = (δ/ǫ)I from the multiplicative bound and have:

(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI.
The idea of using overestimates of leverage scores to perform row sampling has been applied
successfully to various problems (see e.g. [KMP10, CLM+ 15]). However, in these applications,
access to the entire matrix is required beforehand. In the streaming and online settings, we have
to rely on partial data to approximate the true leverage scores. The most natural idea is to just
use the portion of the matrix seen thus far as an approximation to A. This leads us to introduce
the online λ-ridge leverage scores:
li := min(aTi (ATi−1 Ai−1 + λI)−1 ai , 1),
where Ai (i = 0, . . . , n) is defined as the matrix consisting of the first i rows of A1 .
Since clearly ATi Ai  AT A for all i, it is not hard to see that li does overestimate the true
λ-ridge leverage
score for row ai . A more complex question, however, is establishing an upper
P
bound on ni=1 li so that we can bound the number of samples needed by Theorem 2.1.
A core result of this work, stated in Theorem 2.2, is establishing such an upper bound; in fact,
this bound is shown to be tight up to constants (Theorem 5.1) and is nearly-linear in most cases.
Theorem 2.2 Let A be an n × d matrix with rows a1 , . . . , an . Let Ai for i ∈ {0, . . . , n} be the
matrix consisting of the first i rows of A. For λ > 0, let
li := min(aTi (ATi−1 Ai−1 + λI)−1 ai , 1).
be the online λ-ridge leverage score of the ith row of A. Then
n
X
i=1

li = O(d log(kAk22 /λ)).

1
We use the proposed scores li for simplicity, however note that the following, perhaps more natural, definition of
online leverage scores would also be effective:

li′ := aTi (ATi Ai + λI)−1 ai .

4

Theorems 2.2 and 2.1 suggest a simple algorithm for online row sampling: simply use the
online λ-ridge leverage scores, for λ := δ/ǫ. This produces a spectral approximation with only
O(d log d log(ǫkAk22 /δ)/ǫ2 ) rows. Unfortunately, computing li exactly requires us to store all the
rows we have seen in memory (or alternatively to store the sum of their outer products, ATi Ai ).
In many cases, such a requirement would defeat the purpose of streaming row sampling.
A natural idea is to use the sample we have kept thus far as an approximation to Ai when
computing li . It turns out that the approximate online ridge leverage scores ˜li computed in this
way will not always be good approximations to li ; however, we can still prove that they satisfy the
requisite bounds and yield the same row sample size! We formalize these results in the algorithm
Online-Sample (Figure 1) and Theorem 2.3.
Ã = Online-Sample(A, ǫ, δ), where A is an n × d matrix with rows a1 , . . . , an , ǫ ∈ (0, 1),
δ > 0.
1. Set λ := δ/ǫ, c := 8 log d/ǫ2 .
2. Let Ã0 be a 0 × d matrix.
3. For i = 1, . . . , n:
(a) Let ˜li := min((1 + ǫ)aTi (ÃTi−1 Ãi−1 + λI)−1 ai , 1).
(b) Let pi := min(c˜li , 1).
#
"

Ã
i−1


with probability pi ,
√
a
/
pi
i
(c) Set Ãi :=



Ãi−1
otherwise.

4. Return Ã := Ãn .

Figure 1: The basic online sampling algorithm
Theorem 2.3 Let Ã be the matrix returned by Online-Sample(A, ǫ, δ). With high probability,
(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI,
and the number of rows in Ã is O(d log d log(ǫkAk22 /δ)/ǫ2 ).
To save computation, we note that, with a small modification to our analysis, we can run
Online-Sample with batch processing of rows. Specifically, say we start from the ith position in
the stream. we can store the next b = O(d) rows. We can then compute sampling probabilities
for these rows all at once using a system solver for (ÃTi+b Ãi+b + λI). Using a trick introduced in
[SS11], by applying a Johnson-Lindenstrauss random projection to the rows whose scores we are
computing, we need just O(log(1/δ)) system solves to compute constant factor approximations to
the ridge scores with probability 1 − δ. If we set δ = 1/poly(n) then we can union bound over our
whole stream, using this trick with each batch of O(d) input rows. The batch probabilities will
only be closer to the true ridge leverage scores than the non-batch probabilities and we will enjoy
the same guarantees as Online-Sample.
5

Additionally, it turns out that with a simple trick, it is possible to reduce the memory usage
of the algorithm by a factor of ǫ−2 , bringing it down to O(d log d log(ǫkAk22 /δ)) (assuming the row
sample is output to an output stream). Note that this expression gets smaller with ǫ; hence we
obtain a row sampling algorithm with memory complexity independent of desired multiplicative
precision. The basic idea is that, instead of keeping all previously sampled rows in memory, we
store a smaller set of rows that give a constant factor spectral approximation, still enough to give
good estimates of the online ridge leverage scores.
This result is presented in the algorithm Slim-Sample (Figure 2) and Lemma 3.5. A particularly interesting consequence for graphs with polynomially bounded edge weights is:
Corollary 2.4 Let G be a simple graph on d vertices, and ǫ ∈ (0, 1). We can construct a (1 + ǫ)sparsifier of G of size O(d log2 d/ǫ2 ), using only O(d log2 d) working memory in the online model.
2 (A) and noting that the
Proof: This follows simply from applying Theorem 2.3 with δ = ǫ/σmin
condition number of a graph on d vertices whose edge weights are within a multiplicative poly(d)
of each other is polynomial in d. So log(ǫkAk22 /δ) = log(κ2 (A)) = O(log d).


We remark that the algorithm of Corollary 2.4 can be made to run in nearly linear time in the
stream size. We combine Slim-Sample with the batch processing idea described above. Because
A is a graph, our matrix approximation is always a symmetric diagonally dominant matrix, with
O(d) nonzero entries. We can solve systems in it in time Õ(d). Using the Johnson-Lindenstrauss
random projection trick of [SS11], we can compute approximate ridge leverage scores for a batch
of O(d) rows with failure probability polynomially small in n in Õ(d log n) time. Union bounding
over the whole stream, we obtain nearly linear runtime.
To complement the row sampling results discussed above, we explore the limits of the proposed
online setting. In Section 4 we present the algorithm Online-BSS, which obtains spectral approximations with O(d log(ǫkAk22 /δ)/ǫ2 ) rows in the online setting (with larger memory requirements
than the simpler sampling algorithms). Its analysis is given in Theorem 4.1. In Section 5, we show
that this number of samples is in fact the best achievable, up to constant factors (Theorem 5.1).
The log(ǫkAk22 /δ) factor is truly the cost of requiring rows to be selected in an online manner.

3

Analysis of Sampling Schemes

We begin by bounding the sum of online λ-ridge leverage scores. The intuition behind the proof
of Theorem 2.2 is that whenever we add a row with a large online leverage score to a matrix, we
increase its determinant significantly, as follows from the matrix determinant lemma (Lemma 3.1).
Thus we can reduce upper bounding the online leverage scores to bounding the matrix determinant.
Lemma 3.1 (Matrix determinant lemma) Assume S is an invertible square matrix and u is
a vector. Then
det(S + uuT ) = (det S)(1 + uT S−1 u).

6

Proof of Theorem 2.2:

By Lemma 3.1, we have

det(ATi+1 Ai+1 + λI) = det(ATi Ai + λI) · 1 + aTi+1 (ATi Ai + λI)−1 ai+1
≥ det(ATi Ai + λI) · (1 + li+1 )



≥ det(ATi Ai + λI) · eli+1 /2 .
Hence,

det(AT A + λI) = det(ATn An + λI)
P

≥ det(λI) · e
P

= λd e
We have det(AT A + λI) ≤ (kAk22 + λ)d . Therefore

li /2

P

(kAk22 + λ)d ≥ λd e

li /2

.

li /2

.

Taking logarithms of both sides, we obtain
X
li /2,
d log(kAk22 + λ) ≥ d log λ +
X
li ≤ 2d log(1 + kAk22 /λ).



We now turn to analyzing the algorithm Online-Sample. Because the samples taken by the
algorithm are not independent, we are not able to use a standard matrix Chernoff bound like the
one in Theorem 2.1. However, we do know that whether we take row i does not depend on later
rows; thus, we are able to analyze the process as a martingale. We will use a matrix version of the
Freedman inequality given by Tropp.
Theorem 3.2 (Matrix Freedman inequality [Tro11]) Let Y0 , Y1 , . . . , Yn be a matrix martingale whose values are self-adjoint matrices with dimension d, and let X1 , . . . , Xn be the difference
sequence. Assume that the difference sequence is uniformly bounded in the sense that
kXk k2 ≤ R almost surely, for k = 1, . . . , n.
Define the predictable quadratic variation process of the martingale:
Wk :=

k
X
j=1

 
Ej−1 X2j , for k = 1, . . . , n.

Then, for all ǫ > 0 and σ 2 > 0,


P kYn k2 ≥ ǫ and kWn k2 ≤ σ

2



−ǫ2 /2
≤ d · exp − 2
σ + Rǫ/3




We begin by showing that the output of Online-Sample is in fact an approximation of A, and
that the approximate online leverage scores are lower bounded by the actual online leverage scores.
7

Lemma 3.3 After running Online-Sample, it holds with high probability that
(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI,
and also
˜
li ≥ aTi (AT A + λI)−1 ai
for i = 1, . . . , n.
Proof:

Let
ui := (AT A + λI)−1/2 ai .

We construct a matrix martingale Y0 , Y1 , . . . , Yn ∈ Rd×d with the difference sequence X1 , . . . , Xn .
Set Y0 = 0. If kYi−1 k2 ≥ ǫ, we set Xi := 0. Otherwise, let
(
(1/pi − 1)ui uTi
if ai is sampled in Ã,
Xi :=
T
−ui ui
otherwise.
Note that in this case we have
Yi−1 = (AT A + λI)−1/2 (ÃTi−1 Ãi−1 − ATi−1 Ai−1 )(AT A + λI)−1/2 .
Hence, since kYi−1 k2 < ǫ, we have
˜li = min((1 + ǫ)aT (ÃT Ãi−1 + λI)−1 ai , 1)
i
i−1
≥ min((1 + ǫ)aTi (ATi−1 Ai−1 + λI + ǫ(AT A + λI))−1 ai , 1)

≥ min((1 + ǫ)aTi ((1 + ǫ)(AT A + λI))−1 ai , 1)
= aTi (AT A + λI)−1 ai

= uTi ui ,
and so pi ≥ min(cuTi ui , 1). If pi = 1, then Xi = 0. Otherwise, we have pi ≥ cuTi ui and so
kXi k2 ≤ 1/c
and
 
Ei−1 X2i  pi · (1/pi − 1)2 (ui uTi )2 + (1 − pi ) · (ui uTi )2
= (ui uTi )2 /pi

 ui uTi /c.
And so, for the predictable quadratic variation process of the martingale {Yi }:
Wi :=

i
X
k=1

 
Ek−1 X2k ,
8

we have
kWi k2 ≤

i
X

ui uTi /c

k=1

2

≤ 1/c.

Therefore by, Theorem 3.2, we have
P [kYn k2 ≥ ǫ] ≤ d · exp



−ǫ2 /2
1/c + ǫ/(3c)

≤ d · exp(−cǫ2 /4)



= 1/d.
This implies that with high probability

k(AT A + λI)−1/2 (ÃT Ã + λI)(AT A + λI)−1/2 − Ik2 ≤ ǫ
and so
(1 − ǫ)(AT A + λI)  ÃT Ã + λI  (1 + ǫ)(AT A + λI).
Subtracting λI = (δ/ǫ)I from all sides, we get
(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI.

If we set c in Online-Sample to be proportional to log n rather than log d, we would be able to
take a union bound over all the rows and guarantee that with high probability all the approximate
online leverage scores ˜li are close to true online leverage scores li . Thus Theorem 2.2 would imply
that Online-Sample only selects O(d log n log(kAk22 /λ)/ǫ2 ) rows with high probability.
In order to remove the dependency on n, we have to sacrifice achieving close approximations
to li at every step. Instead, we show that the sum of the computed approximate online leverage
scores is still small with high probability, using a custom Chernoff bound.
Lemma 3.4 After running Online-Sample, it holds with high probability that
n
X
i=1

Proof:

˜li = O(d log(kAk2 /λ)).
2

Define
δi := log det(ÃTi Ãi + λI) − log det(ÃTi−1 Ãi−1 + λI).

The proof closely follows the idea from the proof of Theorem 2.2. We will aim to show that large
values of ˜li correlate with large values of δi . However, the sum of δi can be bounded by
h the logarithmi
T
of the ratio of the determinants of Ã Ã + λI and λI. First, we will show that Ei−1 exp(˜li /8 − δi )
is always at most 1. We begin by an application of Lemma 3.1.
h
i
Ei−1 exp(˜li /8 − δi ) = pi · eli /8 (1 + aTi (ÃTi−1 Ãi−1 + λI)−1 ai /pi )−1 + (1 − pi )eli /8
≤ pi · (1 + li /4)(1 + aTi (ÃTi−1 Ãi−1 + λI)−1 ai /pi )−1 + (1 − pi )(1 + li /4).
9

If c˜li < 1, we have pi = c˜li and ˜li = (1 + ǫ)aTi (ÃTi−1 Ãi−1 + λI)−1 ai , and so:
h
i
Ei−1 exp(˜li /8 − δi ) ≤ c˜li · (1 + li /4)(1 + 1/((1 + ǫ)c))−1 + (1 − c˜li )(1 + li /4)
= (1 + li /4)(cli (1 + 1/((1 + ǫ)c))−1 + 1 − cli )

≤ (1 + li /4)(1 + cli (1 − 1/(4c) − 1))

= (1 + li /4)(1 − li /4)
≤ 1.

Otherwise, we have pi = 1 and so:
h
i
Ei−1 exp(˜li /8 − δi ) ≤ (1 + li /4)(1 + aTi (ÃTi−1 Ãi−1 + λI)−1 ai )−1
≤ (1 + li /4)(1 + li )−1

≤ 1.

We will now analyze the expected product of exp(˜li /8 − δi ) over the first k steps. We group the
expectation over the first k steps into one over the first k − 1 steps, aggregating the expectation for
the last step by using one-way independence. For k ≥ 1 we have
!
#
"
!#
"
k−1
k
h
i
X
X
˜li /8 − δi Ek−1 exp(˜lk /8 − δk )
exp
=
E
l˜i /8 − δi
E exp
first
k
−
1
steps
i=1
i=1
!#
"
k−1
X
˜li /8 − δi
,
≤ E exp
i=1

and so by induction on k
"

E exp

n
X
i=1

˜li /8 − δi

!#

≤ 1.

Hence by Markov’s inequality
P

"

n
X

˜li > 8d + 8

n
X
i=1

i=1

#

δi ≤ e−d .

By Lemma 3.3, with high probability we have
ÃT Ã + λI  (1 + ǫ)(AT A + λI).
We also have with high probability
det(ÃT Ã + λI) ≤ (1 + ǫ)d (kAk22 + λ)d ,

log det(ÃT Ã + λI) ≤ d(1 + log(kAk22 + λ)).

10

Hence, with high probability it holds that
n
X
i=1

δi = log det(ÃT Ã + λI) − d log(λ)
≤ d(1 + log(kAk22 + λ) − log(λ))
= d(1 + log(1 + kAk22 /λ)).

And so, with high probability,
n
X
i=1

˜
li ≤ 8d + 8

n
X

δi

i=1

≤ 9d + 8d log(1 + kAk22 /λ)
= O(d log(kAk22 /λ)).


The thesis follows immediately from Lemmas 3.3 and 3.4.

Proof of Theorem 2.3:



We now consider a simple modification of Online-Sample that removes dependency on ǫ from
the working memory usage with no additional cost.
Ã = Slim-Sample(A, ǫ, δ), where A is an n × d matrix with rows a1 , . . . , an , ǫ ∈ (0, 1),
δ > 0.
1. Set λ := δ/ǫ, c := 8 log d/ǫ2 .
2. Let Ã0 be a 0 × d matrix.
3. Let ˜l1 , . . . , ˜ln be the approximate online leverage scores computed by an independent
instance of Online-Sample(A, 1/2, δ/(2ǫ)).
4. For i = 1, . . . , n:
(a) Let pi := min(c˜li , 1).
#
"

Ãi−1


√
ai / p i
(b) Set Ãi :=



Ãi−1

with probability pi ,
otherwise.

5. Return Ã := Ãn .

Figure 2: The low-memory online sampling algorithm
Lemma 3.5 Let Ã be the matrix returned by Slim-Sample(A, ǫ, δ). Then, with high probability,
(1 − ǫ)AT A − δI  ÃT Ã  (1 + ǫ)AT A + δI,
11

and the number of rows in Ã is O(d log d log(ǫkAk22 /δ)/ǫ2 ).
Moreover, with high probability the algorithm Slim-Sample’s memory requirement is dominated
by storing O(d log d log(ǫkAk22 /δ)) rows of A.
Proof: As the samples are independent, the thesis follows from Theorem 2.1 and Lemmas 3.3
and 3.4.


4

Asymptotically Optimal Algorithm

In addition to sampling by online leverage scores, there is also a variant of the “BSS” method
[BSS12] that applies in our setting. Like the original [BSS12], this approach removes the log d
factor from the row count of the output spectral approximation, matching the lower bound for
online sampling given in Theorem 5.1.
Unlike [BSS12] itself, our algorithm is randomized – it is similar to, and inspired by, the randomized version of BSS from [LS15], especially the simpler “Algorithm 1” from that paper (the
main difference from that is considering each row separately). In fact, this algorithm is of the same
form as the basic sampling algorithm, in that when each row comes in, a probability pi is assigned
to it, and it is kept (and rescaled) with probability pi and rejected otherwise. The key difference is
the definition of the pi .
There are also some differences in the nature of the algorithm and its guarantees. Notably, the
pi cannot be computed solely based on the row sample output so far–it is necessary to “remember”
the entire matrix given so far. This means that the BSS method is not memory efficient, using O(d2 )
space. Additionally, online leverage score sampling gives bounds on both the size of the output
spectral approximation and its accuracy with high probability. In contrast, this method gives an
expected bound on the output size, while it never fails to output a correct spectral approximation.
Note that these guarantees are essentially the same as those in the appendix of [LS15].
One may, however, improve the memory dependence in some cases simply by running it on the
output stream of the online leverage score sampling method. This reduces the storage cost to the
size of that spectral approximation. The BSS method still does not produce an actual space savings
(in particular, there is a still a log d factor in space), but it does reduce the number of rows in the
output stream while only blowing up the space usage by O(1/ǫ2 ) (due to requiring the storage of
an ǫ-quality approximation rather than only O(1)).
L
The BSS method maintains two matrices, BU
i and Bi , acting as upper and lower “barriers”.
The current spectral approximation will always fall between them:
U
T T
BL
i ≺ Ãi Ãi ≺ Bi .

This guarantee, at the end of the algorithm, will ensure that Ã is a valid approximation.
Below, we give the actual BSS algorithm and its performance guarantees.
Theorem 4.1
1. The online BSS algorithm always outputs Ã such that
(1 − ǫ)AT A − δI ≺ ÃT ÃT ≺ (1 + ǫ)AT A + δI
12

2. The probability that a row ai is included in Ã is at most ǫ82 li , where li is the online 2δ
ǫ -ridge
−1
2δ
T
T
leverage score of ai . That is li = min(ai Ai Ai + ǫ I
ai , 1). The expected number of rows
8 Pn
2
2
in Ã is thus at most ǫ2 i=1 li = O(d log(ǫkAk2 /δ)/ǫ ).
Ã = Online-BSS(A, ǫ, δ), where A is an n × d matrix with rows a1 , . . . , an , ǫ ∈ (0, 1),
δ > 0.
1. Set cU =

2
ǫ

+ 1 and cL =

2
ǫ

− 1.

L
2. Let Ã0 be a 0 × d matrix, BU
0 = δI, B0 = −δI.

3. For i = 1, . . . , n:
L
L
T
U
T
(a) Let XU
i−1 = (Bi−1 − Ãi−1 Ãi−1 ), Xi−1 = (Ãi−1 Ãi−1 − Bi−1 ).

−1
T
L
−1
(b) Let pi := min(cU aTi (XU
i−1 ) ai + cL ai (Xi−1 ) ai , 1).
#
"

 Ãi−1

with probability pi ,
√
ai / p i
(c) Set Ãi :=



Ãi−1
otherwise.

U
T
L
L
T
(d) Set BU
i = Bi−1 + (1 + ǫ)ai ai , Bi = Bi−1 + (1 − ǫ)ai ai .

4. Return Ã := Ãn .
Figure 3: The Online BSS Algorithm
L
We first note the basic invariant that XU
Proof of Theorem 4.1 part 1:
i and Xi always
remain positive definite–or equivalently,
U
T T
BL
i ≺ Ãi Ãi ≺ Bi .

We may prove this by induction on i. The base case follows from the initialization of Ã0 , BU
0
.
For
each
successive
step,
we
consider
two
possibilities.
and BL
0
The first is that pi = 1. In that case, ÃT Ã always increases by exactly ai aTi , BU by (1 + ǫ)ai aTi
and BL by (1 − ǫ)ai aTi . Thus XU and XL increase by exactly ǫai aTi , which is positive semidefinite,
and so remain positive definite.
In the other case, pi < 1. Now, XU decreases by at most the increase in ÃTi ÃTi , or
Mi =

ai aTi
.
p

−1
T
U
U
Since cU > 1, p > aTi (XU
i−1 ) ai , so ai ai ≺ pXi−1 and Mi ≺ Xi−1 . Subtracting this then
U
L
must leave X positive definite. Similarly, X decreases by at most the increase in BL , which is
L
−1
T
(1 − ǫ)ai aTi ≺ ai aTi . Since cL > 1 and p < 1, aTi (XL
i−1 ) ai < 1, and ai ai ≺ Xi−1 . Subtracting
L
this similarly leaves X positive definite. Finally, we note that
T
BU
n = (1 + ǫ)A A + δI
T
BL
n = (1 − ǫ)A A − δI.

13

This gives the desired result.



To prove part 2, we will use quantities of the form vT X−1 v. We will need a lemma describing
how this behaves under a random rank-1 update:
Lemma 4.2 Given a positive definite matrix X, two vectors u and v, two multipliers a and b and
a probability p, define the random variable X′ to be X − auuT with probability p and X − buuT
otherwise. Then if uT X−1 u = 1,


T ′−1
T −1
T −1 2 pa + (1 − p)b − ab
E v X v − v X v] = (v X u)
.
(1 − a)(1 − b)
Proof: We apply the Sherman-Morrison formula to each of the two possibilities (subtracting
auuT and buuT respectively). These give X′ values of respectively
X−1 + a
and

a
X−1 uuT X−1
= X−1 +
X−1 uuT X−1
T
−1
1 − au X u
1−a

b
X−1 uuT X−1
= X−1 +
X−1 uuT X−1 .
T
−1
1 − bu X u
1−b
The values of vT X′−1 v − vT X−1 v are then respectively
a
a
vT X−1 uuT X−1 v = (vT X−1 u)2
1−a
1−a
and
b
b
vT X−1 uuT X−1 v = (vT X−1 u)2
.
1−b
1−b
Combining these gives the stated result.
X−1 + b



First, we introduce some new matrices to help in the analysis:

ǫ T
ǫ T
CU
=
δI
+
A Aj
A
A
+
1
+
i
i,j
2 i
2 j

ǫ T
ǫ
CL
ATj Aj .
i,j = −δI − Ai Ai + 1 −
2
2
U
U
L
L
= BL
i , and for j ≤ i, Ci,j  Bj and Ci,j  Bj . We can then define:

Proof of Theorem 4.1 part 2:

U
L
Note that CU
i,i = Bi , Ci,i

U
T
Yi,j
= CU
i,j − Ãj Ãj

L
Yi,j
= ÃTj Ãj − CL
i,j .

U = XU , Y L = XL , and for j ≤ i, Y U  XU and Y L  XL .
We then have, similarly, Yi,i
i
i,i
i
i,j
j
i,j
j
We will assume that li < 1, since otherwise the claim is immediate (as probabilities cannot
exceed 1). Now, note that
U −1
L −1
aTi (Yi,0
) ai = aTi (Yi,0
) ai
ǫ
−1
= aTi
ATi Ai + δI
ai
2
−1

2δ
2
ai
ATi Ai + I
=
ǫ
ǫ
2
= li .
ǫ

14

Next, we will aim to show that for j < i − 1,




U
U
ai
E aTi Yi−1,j+1
ai ≤ E aTi Yi−1,j




L
L
ai
E aTi Yi−1,j+1
ai ≤ E aTi Yi−1,j

In particular, we will simply show that conditioned on any choices for the first j rows, the
U
U
expected value of aTi Yi−1,j+1
ai is no larger than aTi Yi−1,j
ai , and analogously for Y L .
Similar to the proof of part 1, we separately consider the case where pj+1 = 1. In that case, the
positive semidefinite matrix 2ǫ aj aTj is simply added to Y U and Y L . Adding this can only decrease
the values of aTi Y U ai and aTi Y L ai .
a
The pj+1 < 1 case is more tricky. Here, we define the vector wj+1 = √pj+1
. Importantly
j+1
−1
T
U
−1
pj+1 ≥ cU aTj+1 (XU
j ) aj+1 ≥ cU aj+1 (Yi−1,j ) aj+1

−1
T
L
−1
pj+1 ≥ cL aTj+1 (XL
j ) aj+1 ≥ cL aj+1 (Yi−1,j ) aj+1 .

This means that
1
cU
1
.
≤
cL

T
U
T
wj+1
(Yi−1,j
)−1 wj+1
≤
T
L
T
wj+1
(Yi−1,j
)−1 wj+1

Now, we additionally define
T
U
−1 T
sU
j+1 = wj+1 (Yi−1,j ) wj+1
T
L
−1 T
sL
j+1 = wj+1 (Yi−1,j ) wj+1
wj+1
uU
j+1 = q
sU
j+1
w
j+1
uL
.
j+1 = q
sL
j+1

We then deploy Lemma 4.2 to compute the expectations. For the contribution from the upper
U
T
U
U
barrier, we use X = Yi−1,j
, u = uU
j+1 , v = ai , a = −sj+1 (1 − pj+1 (1 + ǫ/2)), b = sj+1 pj+1 (1 + ǫ/2),
L
T
L
p = pj+1 . For the lower barrier, we use X = Yi−1,j
, u = uL
j+1 , v = ai , a = sj+1 (1 − pj+1 (1 − ǫ/2)),
b = −sL
j+1 pj+1 (1 − ǫ/2), p = pj+1 . In both cases we can see that the numerator of the expected
change is nonpositive. Finally, this implies that the probability that row i is sampled is

 T L −1 

−1
E [pi ] = cU E aTi (XU
i−1 ) ai + cL E ai (Xi−1 ) ai




U
L
= cU E aTi (Yi−1,i−1
)−1 ai + cL E aTi (Yi−1,i−1
)−1 ai




L
U
)−1 ai
≤ cU E aTi (Yi−1,0
)−1 ai + cL E aTi (Yi−1,0
2
(cU + cL )li
ǫ
8
= 2 li
ǫ
=

as desired.



15

5

Matching Lower Bound

Here we show that the row count obtained by Theorem 4.1 is in fact optimal. While it is possible
to obtain a spectral approximation with O(d/ǫ2 ) rowsin the offline setting, online sampling always

d log(ǫkAk22 /δ)
incurs a loss of Ω log(ǫkAk22 /δ) and must sample Ω
rows.
ǫ2
√
Theorem 5.1 Assume that ǫkAk22 ≥ c1 δ and ǫ ≥ c2 / d, for fixed constants c1 and c2 . Then
any algorithm that selects rows in an online manner and outputs a spectral approximation to AT A
with
error and δ additive error with probability at least 1/2 must sample
 (1 + ǫ)2 multiplicative

d log(ǫkAk2 /δ)
Ω
rows of A in expectation.
ǫ2

Note that the lower bounds we assume on ǫkAk22 and ǫ are very minor. They just ensure that
log(ǫkAk22 /δ) ≥ 1 and that ǫ is not so small that we can essentially sample all rows of A.
Proof: We apply Yao’s minimax principle, constructing, for any large enough M , a distribution on
inputs A with kAk22 ≤ M for which any deterministic
online
row selection algorithm that succeeds



with probability at least 1/2 must output Ω

d log(ǫM/δ)
ǫ2

rows in expectation. The best randomized

algorithm that works
with probability 1/2 on any input matrix with kAk22 ≤ M therefore must

select at least Ω d log(ǫM/δ)
rows in expectation on the worst case input, giving us the theorem.
ǫ2
Our distribution is as follows. We select an integer N uniformly at random from [1, log(M ǫ/δ)].
We then stream in the vertex edge incidence matrices of N complete graphs on d vertices. We
double the weight of each successive graph. Intuitively, spectrally
approximating a complete graph
√
2
requires selecting Ω(d/ǫ ) edges [BSS12] (as long as ǫ ≥ c2 / d for some fixed constant c2 ). Each
time we stream in a new graph with double the weight, we force the algorithm to add Ω(d/ǫ2 )
more edges to its output, eventually forcing it to output Ω(d/ǫ2 · N ) edges – Ω(d log(M ǫ/δ)/ǫ2 ) in
expectation.

Specifically, let Kd be the d2 × d vertex edge incidence matrix of the complete graph on d
vertices. KTd Kd is the Laplacian matrix of the complete graph on d vertices. We weight the
first graph so that its Laplacian has all its nonzero eigenvalues equal to δ/ǫ. (That is, each edge
δ
). In this way, even if we select N = ⌊log(M ǫ/δ)⌋ we will have overall kAk22 ≤
has weight dǫ
δ/ǫ + 2δ/ǫ + ...2⌊log(M ǫ/δ)⌋−1 δ/ǫ ≤ M .
Even if N = 1, all nonzero eigenvalues of AT A are at least δ/ǫ, so achieving (1+ǫ) multiplicative
error and δI additive error is equivalent to achieving (1 + 2ǫ) multiplicative error. AT A is a graph
Laplacian so has a null space. However, as all rows are orthogonal to the null space, achieving
additive error δI is equivalent to achieving additive error δIr where Ir is the identity projected to
the span of AT A. δIr  ǫAT A which is why we must achieve (1 + 2ǫ) multiplicative error.
In order for a deterministic algorithm to be correct with probability 1/2 on our distribution, it
must be correct for at least 1/2 of our ⌊log(M ǫ/δ)⌋ possible choices of N .
Let i be the lowest choice of N for which the algorithm is correct. By the lower bound of
[BSS12], the algorithm must output Ω(d/ǫ2 ) rows of Ai to achieve a (1 + 2ǫ) multiplicative factor
spectral approximation. Here Ai is the input consisting of the vertex edge incidence matrices of i
increasingly weighted complete graphs. Call the output on this input Ãi . Now let j be the second
lowest choice of N on which the algorithm is correct. Since the algorithm was correct on Ai to
within a multiplicative (1 + 2ǫ), to be correct on Aj , it must output a set of edges Ãj such that
(ATj Aj − ATi Ai ) − 4ǫATj Aj  ÃTj Ãj − ÃTi Ãi  (ATj Aj − ATi Ai ) + 4ǫATj Aj .
16

Since we double each successive copy of the complete graph, ATj Aj  2(ATj Aj − ATi Ai ). So,
ÃTj Ãj − ÃTi Ãi must be a 1 + 8ǫ spectral approximation to the true difference ATj Aj − ATi Ai .
Noting that this difference is itself just a weighting of the complete graph, by the lower bound in
[BSS12] the algorithm must select Ω(d/ǫ2 ) additional edges between the ith and j th input graphs.
Iterating this argument over all ⌊log(M ǫ/δ)⌋/2 inputs on which the algorithm must be correct, it
must select a total of Ω(d log(M ǫ/δ)/ǫ2 ) edges in expectation over all inputs.


6

Future Work

An obvious open question arising from our work is if one can prove that the algorithm of [KL13]
works despite dependencies arising due to the row pruning step. By operating in the online setting,
our algorithm avoids row pruning, and hence is able to skirt these dependencies, as the probability
that a row is sampled only depends on earlier rows in the stream. However, because the streaming
setting offers the potential for sampling fewer rows than in the online case, obtaining a rigorous
proof of [KL13] would be very interesting.
While our work focuses on spectral approximation, variants on (ridge) leverage score sampling
and the BSS algorithm are also used to solve low-rank approximation problems, including column
subset selection [BW14, CMM15] and projection-cost-preserving sketching [CEM+ 15, CMM15].
Compared with spectral approximation, there is less work on streaming sampling for low-rank
approximation, and understanding how online algorithms may be used in this setting would an
interesting extension of our work.

7

Acknowledgments

The authors would like to thank Kenneth Clarkson, Jonathan Kelner, Gary Miller, Christopher
Musco and Richard Peng for helpful discussions and comments.
Cameron Musco and Jakub Pachocki both acknowledge the Gene Golub SIAM Summer School
program on Randomization in Numerical Linear Algebra, where work on this project was initiated.

References
[AM]

Ahmed Alaoui and Michael W Mahoney. Fast randomized kernel ridge regression
with statistical guarantees. In Advances in Neural Information Processing Systems
28 (NIPS), pages 775–783. 1.3, 2

[BB05]

Antoine Bordes and Léon Bottou. The huller: a simple and efficient online SVM. In
Machine Learning: ECML 2005, pages 505–512. Springer, 2005. 1.2

[BGKL15] Christos Boutsidis, Dan Garber, Zohar Karnin, and Edo Liberty. Online principal
components analysis. In Proceedings of the 26th Annual ACM-SIAM Symposium on
Discrete Algorithms (SODA), pages 887–901, 2015. 1.2
[BSS12]

Joshua Batson, Daniel A Spielman, and Nikhil Srivastava. Twice-ramanujan sparsifiers.
SIAM Journal on Computing, 41(6):1704–1721, 2012. 1.3, 4, 5

17

[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. 6

[CDK+ 06] Koby Crammer, Ofer Dekel, Joseph Keshet, Shai Shalev-Shwartz, and Yoram Singer.
Online passive-aggressive algorithms. The Journal of Machine Learning Research,
7:551–585, 2006. 1.2
[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. 6
[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. 1.1, 1.3, 2
[CMM15] Michael B Cohen, Cameron Musco, and Christopher Musco. Ridge leverage scores for
low-rank approximation. arXiv:1511.07263, 2015. 6
[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. 1.1

[KL13]

Jonathan A Kelner and Alex Levin. Spectral sparsification in the semi-streaming setting.
Theory of Computing Systems, 53(2):243–262, 2013. (document), 1.2, 1.3, 6

[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. (document), 1.2
[KMP10]

Ioannis Koutis, Gary L Miller, and Richard Peng. Approaching optimality for solving
SDD linear systems. In Proceedings of the 51st Annual IEEE Symposium on Foundations
of Computer Science (FOCS), pages 235–244, 2010. 2

[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. 1.1

[LS15]

Yin Tat Lee and He Sun. Constructing linear-sized spectral sparsification in almostlinear time. In Proceedings of the 56th Annual IEEE Symposium on Foundations of
Computer Science (FOCS), pages 250–269, 2015. 4

[LSS16]

Edo Liberty, Ram Sriharsha, and Maxim Sviridenko. An algorithm for online k-means
clustering. In Proceedings of the Eighteenth Workshop on Algorithm Engineering and
Experiments (ALENEX), pages 81–89, 2016. 1.2

[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.
1.1

18

[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. 1.1

[SS11]

Daniel A Spielman and Nikhil Srivastava. Graph sparsification by effective resistances.
SIAM Journal on Computing, 40(6):1913–1926, 2011. 2, 2

[ST04]

Daniel A Spielman and Shang-Hua Teng. Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems. In Proceedings of the 36th
Annual ACM Symposium on Theory of Computing (STOC), pages 81–90, 2004. 1.1

[Tro11]

Joel Tropp. Freedman’s inequality for matrix martingales. Electronic Communications
in Probability, 16:262–270, 2011. 3.2

19

