arXiv:1702.08415v1 [cs.DS] 27 Feb 2017

An SDP-Based Algorithm for
Linear-Sized Spectral Sparsification
Yin Tat Lee
Microsoft Research
Redmond, USA
yile@microsoft.com

He Sun
The University of Bristol
Bristol, UK
h.sun@bristol.ac.uk

Abstract
For any undirected and weighted graph G = (V, E, w) with n vertices and m edges, we call
a sparse subgraph H of G, with proper reweighting of the edges, a (1 + ε)-spectral sparsifier if
(1 − ε)x⊺ LG x 6 x⊺ LH x 6 (1 + ε)x⊺ LG x
holds for any x ∈ Rn , where LG and LH are the respective Laplacian matrices of G and H.
Noticing that Ω(m) time is needed for any algorithm to construct a spectral sparsifier and
a spectral sparsifier of G requires Ω(n) edges, a natural question is to investigate, for any
e
constant ε, if a (1 + ε)-spectral sparsifier of G with O(n) edges can be constructed in O(m)
time,
e
where the O notation suppresses polylogarithmic factors. All previous constructions on spectral
sparsification [ST11, SS11, BSS12, Zou12, ALO15, LS15] require either super-linear number of
edges or m1+Ω(1) time.
In this work we answer this question affirmatively by presenting an algorithm that, for any
undirected graph G and ε > 0, outputs a (1 + ε)-spectral sparsifier of G with O(n/ε2 ) edges
O(1)
e
in O(m/ε
) time. Our algorithm is based on three novel techniques: (1) a new potential
function which is much easier to compute yet has similar guarantees as the potential functions
used in previous references; (2) an efficient reduction from a two-sided spectral sparsifier to a
one-sided spectral sparsifier; (3) constructing a one-sided spectral sparsifier by a semi-definite
program.
Keywords: spectral graph theory, spectral sparsification

1

Introduction

A sparse graph is one whose number of edges is reasonably viewed as being proportional to the
number of vertices. Since most algorithms run faster on sparse instances of graphs and it is more
space-efficient to store sparse graphs, it is useful to obtain a sparse representation H of G so that
certain properties between G and H are preserved, see Figure 1 for an illustration. Over the past
three decades, different notions of graph sparsification have been proposed and widely used to design
approximation algorithms. For instance, a spanner H of a graph G is a subgraph of G so that the
shortest path distance between any pair of vertices is approximately preserved [Che89]. Benczúr
and Karger [BK96] defined a cut sparsifier of a graph G to be a sparse subgraph H such that
the value of any cut between G and H are approximately the same. In particular, Spielman and
Teng [ST11] introduced a spectral sparsifer, which is a sparse subgraph H of an undirected graph G
such that many spectral properties of the Laplacian matrices between G and H are approximately
preserved. Formally, for any undirected graph G with n vertices and m edges, we call a subgraph
H of G, with proper reweighting of the edges, a (1 + ε)-spectral sparsifier if
(1 − ε)x⊺ LG x 6 x⊺ LH x 6 (1 + ε)x⊺ LG x
holds for any x ∈ Rn , where LG and LH are the respective Laplacian matrices of G and H. Spectral
sparsification has been proven to be a remarkably useful tool in algorithm design, linear algebra,
combinatorial optimisation, machine learning, and network analysis.
2

3

2

1

3

1
4
6

5

4

7

8

14

13

14

13

15

11

16

16

19

19

18

18

20
21

17
22

10

12

15

11

9
8

10

12

7

6

5

9

20
21

17

23

22

24

23

24
25

26

25
26

28
27
29

G

28
27
29

H

Figure 1: The graph sparsification is a reweighted subgraph H of an original graph G such that certain properties are preserved.
These subgraphs are sparse, and are more space-efficient to be stored than the original graphs. The picture above uses the
thickness of edges in H to represent their weights.

In the seminal work on spectral sparsification, Spielman and Teng [ST11] showed that, for any
undirected graph G of n vertices, a spectral sparsifier of G with only O(n logc n/ε2 ) edges exists
and can be constructed in nearly-linear time1 , where c > 2 is some constant. Both the runtime of
their algorithm and the number of edges in the output graph involve large poly-logarithmic factors,
and this motivates a sequence of simpler and faster constructions of spectral sparsifiers with fewer
edges [BSS12, ALO15, LS15]. In particular, since any constant-degree expander graph of O(n)
edges is a spectral sparsifier of an n-vertex complete graph, a natural question is to study, for any
n-vertex undirected graph G and constant ε > 0, if a (1 + ε)-spectral sparsifier of G with O(n)
edges can be constructed in nearly-linear time. Being considered as one of the most important
open question about spectral sparsification by Batson et al. [BSST13], there has been many efforts
1
We say a graph algorithm runs in nearly-linear time if the algorithm runs in O(m · poly log n) time, where m and
n are the number of edges and vertices of the input graph.

1

for fast constructions of linear-sized spectral sparsifiers, e.g. [ALO15, LS15], however the original
problem posed in [BSST13] has remained open.
In this work we answer this question affirmatively by presenting the first nearly-linear time
algorithm for constructing a linear-sized spectral sparsifier. The formal description of our result is
as follows:
Theorem 1.1. Let G be any undirected graph
with n vertices

 and m edges. For any 0 < ε < 1, there
e m/εO(1) work, O
e 1/εO(1) depth, and produces a (1 + ε)-spectral
is an algorithm that runs inO
sparsifier of G with O n/ε2 edges2 .

Theorem 1.1 shows that a linear-sized spectral sparsifier can be constructed in nearly-linear time
in a single machine setting, and in polylogarithmic time in a parallel setting. The same algorithm
can be applied to the matrix setting, whose result is summarised as follows:
Pm
m
n×n . Let M =
Theorem
i=1 Mi and
Pm 1.2. Given a set of m PSD matrices {Mi }i=1 , where Mi ∈ R
Z = i=1 nnz(Mi ), where nnz(Mi ) is the number of non-zero entries in Mi . For any 1 > ε > 0,


e (Z + nω )/εO(1) work, O
e 1/εO(1) depth and produces a (1+ε)there is an algorithm that runs in O

2 components, i.e., there is an non-negative coefficients {c }m
spectral sparsifier of M with O n/ε
i i=1

2
such that |{ci |ci 6= 0}| = O n/ε , and
(1 − ε) · M 

m
X

ci Mi  (1 + ε) · M.

(1.1)

i=1

Here ω is the matrix multiplication constant.

1.1

Related work

In the seminal paper on spectral sparsification, Spielman and Teng [ST11] showed that a spectral
sparsifier of any undirected graph G can be constructed by decomposing G into multiple nearly
expander graphs, and sparsifying each subgraph individually. This method leads to the first nearlylinear time algorithm for constructing a spectral sparsifier with O(n logc n/ε2 ) edges for some c >
2. However, both the runtime of their algorithm and the number of edges in the output graph
involve large poly-logarithmic factors. Spielman and Srivastava [SS11] showed that a (1 + ε)spectral sparsifier of G with O(n log n/ε2 ) edges can be constructed by sampling the edges of G
with probability proportional to their effective resistances, which is conceptually much simpler than
the algorithm presented in [ST11].
Noticing that any constant-degree expander graph of O(n) edges is a spectral sparsifier of an nvertex complete graph, Spielman and Srivastava [SS11] asked if any n-vertex graph has a spectral
sparsifier with O(n) edges. To answer this question, Batson, Spielman and Srivastava [BSS12]
presented a polynomial-time algorithm that, for any undirected graph G of n vertices, produces a
spectral sparsifier of G with O(n) edges. At a high level, their algorithm, a.k.a. the BSS algorithm,
proceeds for O(n) iterations, and in each iteration one edge is chosen deterministically to “optimise”
the change of some potential function. Allen-Zhu et al. [ALO15] noticed that a less “optimal” edge,
based on a different potential function, can be found in almost-linear time and this leads to an
almost-quadratic time algorithm. Generalising their techniques, Lee
 and Sun [LS15] showed that a
linear-sized spectral sparsifier can be constructed in time O m1+c for an arbitrary small constant c.
All of these algorithms proceed for Ω(nc ) iterations, and every iteration takes Ω(m1+c ) time for some
constant c > 0. Hence, to break the Ω(m1+c ) runtime barrier faced in all previous constructions
multiple new techniques are needed.
2

e hides a factor of logc n for some positive constant c.
Here, the notation O(·)

2

1.2

Organisation

The remaining part of the paper is organised as follows. We introduce necessary notions about
matrices and graphs in Section 2. In Section 3 we overview our algorithm and proof techniques.
For readability, more detailed discussions and technical proofs are deferred to Section 4.

2
2.1

Preliminaries
Matrices

For any n × n real and symmetric matrix A, let λmin (A) = λ1 (A) 6 · · · 6 λn (A) = λmax (A) be the
eigenvalues of A, where λmin (A) and λmax (A) represent the minimum and maximum eigenvalues of
A. We call a matrix A positive semi-definite (PSD) if x⊺ Ax > 0 holds for any x ∈ Rn , and a matrix
n
A positive definite if x⊺ Ax > 0 holds for any
 x ∈⊺ R−1. For any positive definite matrix A, we define
the corresponding ellipsoid by Ellip(A) , x : x A x 6 1 .

2.2

Graph Laplacian

Let G = (V, E, w) be a connected, undirected and weighted graph with n vertices, m edges, and
weight function w : E → R>0 . We fix an arbitrary orientation of the edges in G, and let B ∈ Rm×n
be the signed edge-vertex incidence matrix defined by

if v is e’s head,

 1
BG (e, v) = −1
if v is e’s tail,


0
otherwise.

We define an m × m diagonal matrix WG by WG (e, e) = we for any edge e ∈ E[G].
The Laplacian matrix of G is an n × n matrix L defined by

if u ∼ v,

 −w(u, v)
deg(u)
if u = v,
LG (u, v) =


0
otherwise,
P
where deg(v) = u∼v w(u, v). It is easy to verify that
X
⊺
W G BG x =
wu,v (xu − xv )2 > 0
x⊺ L G x = x⊺ BG
u∼v

holds for any x ∈ Rn . Hence, the Laplacian matrix of any undirected graph is a PSD matrix. Notice
that, by setting xu = 1 if u ∈ S and xu = 0 otherwise, x⊺ LG x equals to the value of the cut between
S and V \ S. Hence, a spectral sparsifier is a stronger notion than a cut sparsifer.

2.3

Other notations

m
For any sequence {αi }m
i=1 , we use nnz(α) to denote the number of non-zeros in {αi }i=1 . For any
two matrices A and B, we write A  B to represent B − A is PSD, and A ≺ B to represent B − A
is positive definite. For any two matrices A and B of the same dimension, let A • B , tr (A⊺ B),
and


A 0
A⊕B =
.
0 B

3

3

Overview of our algorithm

Without loss of generality we study the problem of sparsifying the sum of PSD matrices. The oneto-one correspondence between the construction of a graph sparsifier and the following Problem 1
was presented in [BSS12].
P
Problem 1. Given a set S of m PSD matrices M1 , · · · , Mm with m
i=1 Mi = I and 0 < ε < 1,
m
2
find non-negative coefficients {ci }i=1 such that |{ci |ci 6= 0}| = O n/ε , and
(1 − ε) · I 

m
X

ci Mi  (1 + ε) · I.

(3.1)

i=1

For intuition, one can think all Mi are rank-1 matrices, i.e., Mi = vi vi⊺ for some vi ∈ Rn . Given
the correspondence between PSD matrices and ellipsoids, Problem 1 essentially asks to use O(n/ε2 )
vectors from S to construct an ellipsoid, whose shape is close to be a sphere. To construct such an
ellipsoid with desired shape, all previous algorithms [BSS12, ALO15, LS15] proceed by iterations:
in each iteration j the algorithm chooses one or more vectors, denoted by vj1 , · · · , vjk , and adds
P
∆j , kt=1 vjt vj⊺t to the currently constructed matrix by setting Aj = Aj−1 + ∆j . To control the
shape of the constructed ellipsoid, two barrier values, the upper barrier uj and the lower barrier
ℓj , are maintained such that the constructed ellipsoid Ellip(Aj ) is sandwiched between the outer
sphere uj · I and the inner sphere ℓj · I for any iteration j. That is, the following invariant always
maintains:
ℓj · I ≺ Aj ≺ uj · I.
(3.2)
To ensure (3.2) holds, two barrier values ℓj and uj are increased properly after each iteration, i.e.,
uj+1 = uj + δu,j ,

ℓj+1 = ℓj + δℓ,j

for some positive values δu,j and δℓ,j . The algorithm continues this process, until after T iterations
Ellip(AT ) is close to be a sphere. This implies that AT is a solution of Problem 1, see Figure 2 for
an illustration.

Iteration j

Iteration j + 1

Final iteration T

Figure 2: Illustration of the algorithms for constructing a linear-sized spectral sparsifier. Here, the light grey and green circles
in iteration j represent the spheres ℓj · I and uj · I, and the dark grey ellipse sandwiched between the two circles corresponds to
the constructed ellipsoid in iteration j. After each iteration j, the algorithm increases the value of ℓj and uj by some δℓ,j and
δu,j so that the invariant (3.2) holds in iteration j + 1. This process is repeated for T iterations, so that the final constructed
ellipsoid is close to be a sphere.

However, turning the scheme described above into an efficient algorithm we need to consider
the following issues:

4

• Which vectors should we pick in each iteration?
• How many vectors can be added in each iteration?
• How should we update uj and ℓj properly so that the invariant (3.2) always holds?
These three questions closely relate to each other: on one hand, one can always pick a single
“optimal” vector in each iteration based on some metric, and such conservative approach requires
a linear number of iterations T = Ω(n/ε2 ) and super-quadric time for each iteration. On the
other hand, one can choose multiple less “optimal” vectors to construct ∆j in iteration j, but this
makes the update of barrier values more challenging to ensure the invariant (3.2) holds. Indeed,
the previous constructions [ALO15, LS15] speed up their algorithms at the cost of increasing the
sparsity, i.e., the number of edges in a sparsifier, by more than a multiplicative constant.
To address these, we introduce three novel techniques for constructing a spectral sparsifier: First
of all, we define a new potential function which is much easier to compute yet has similar guarantee
as the potential function introduced in [BSS12]. Secondly we show that solving Problem 1 with
two-sided constraints in (3.1) can be reduced to a similar problem with only one-sided constraints.
Thirdly we prove that the problem with one-sided constraints can be solved by a semi-definite
program.

3.1

A new potential function

To ensure that the constructed ellipsoid A is always inside the outer sphere u · I, we introduce a
potential function Φu (A) defined by
−1

Φu (A) , tr exp (uI − A)



=

n
X

exp

i=1



1
u − λi (A)



.

It is easy to see that, when Ellip(A) gets closer to the outer sphere, λi (u · I − A) becomes smaller
and the value of Φu (A) increases. Hence, a bounded value of Φu (A) ensures that Ellip(A) is inside
the sphere u · I. For the same reason, we introduce a potential function Φℓ (A) defined by
−1

Φℓ (A) , tr exp (A − ℓI)



=

n
X
i=1

exp



1
λi (A) − ℓ



.

to ensure that the inner sphere is always inside Ellip(A). We also define
Φu,ℓ (A) , Φu (A) + Φℓ (A),

(3.3)

as a bounded value of Φu,ℓ (A) implies that the two events occur simultaneously. Our goal is to
design a proper update rule to construct {Aj } inductively, so that Φuj ,ℓj (Aj ) is monotone nonincreasing after each iteration. Assuming this, a bounded value of the initial potential function
guarantees that the invariant (3.2) always holds.
To analyse the change of the potential function, we first notice that




−1
−1
Φu,ℓ (A + ∆) > Φu,ℓ (A) + tr e(uI−A) (uI − A)−2 ∆ − tr e(A−ℓI) (A − ℓI)−2 ∆

by the convexity of the function Φu,ℓ . We prove that, as long as the matrix ∆ satisfies 0  ∆ 
δ(uI − A)2 and 0  ∆  δ(A − ℓI)2 for some small δ, the first-order approximation gives a good
approximation.
5

Lemma 3.1. Let A be a symmetric matrix. Let u, ℓ be the barrier values such that u − ℓ 6 1 and
ℓI ≺ A ≺ uI. Assume that ∆ ≻ 0, ∆  δ(uI − A)2 and ∆  δ(A − ℓI)2 for δ 6 1/10. Then, it
holds that


−1
Φu,ℓ (A + ∆) 6Φu,ℓ (A) + (1 + 2δ)tr e(uI−A) (uI − A)−2 ∆


−1
− (1 − 2δ)tr e(A−ℓI) (A − ℓI)−2 ∆ .

We remark that this is not the first paper to use a potential function to guide the growth of the
ellipsoid. In [BSS12], the potential function


Λu,ℓ,p (A) = tr (uI − A)−p + tr (A − ℓI)−p
(3.4)

is used with p = 1. The main drawback is that Λu,ℓ,1 does not differentiate the following two cases:
• Multiple eigenvalues of A are close to the boundary (both u and ℓ).
• One of the eigenvalues of A is very close to the boundary (either u or ℓ).

It is known that, when one of the eigenvalues of A is very close to the boundary, it is more difficult to
find an “optimal” vector. [ALO15] shows this problem can be alleviated by using p ≫ 1. However,
this choice of p makes the function Λu,ℓ,p less smooth and hence one has to take a smaller step size
δ, as shown in the following lemma by [ALO15].
Lemma 3.2 ([ALO15]). Let A be a symmetric matrix and ∆ be a rank-1 matrix. Let u, ℓ be the
barrier values such that ℓI ≺ A ≺ uI. Assume that ∆ ≻ 0, ∆  δ(uI − A) and ∆  δ(A − ℓI) for
δ 6 1/(10p) and p > 10. Then, it holds that




Λu,ℓ,p (A + ∆) 6Λu,ℓ,p (A) + p(1 + pδ)tr (uI − A)−(p+1) ∆ − p(1 − pδ)tr (A − ℓI)−(p+1) ∆ .

Notice that, comparing with the potential function (3.4), our new potential function (3.3) blows
up much faster when the eigenvalues of A are closer to the boundaries ℓ and u. This allows the
problem of finding an “optimal” vector much easier than using Λu,ℓ,p . At the same time, we avoid
the problem of taking a small step ∆  1/p · (uI − A) by taking a “non-linear” step ∆  (uI − A)2 .
As there cannot be too many eigenvalues close to the boundaries, this “non-linear” step allows us
to take a large step except on a few directions.

3.2

A simple construction based on Oracle

The second technique we introduce is the reduction from a spectral sparsifier with two-sided constraints to the one with one-sided constraints. Geometrically, it is equivalent to require the constructed ellipsoid inside another ellipsoid, instead of being sandwiched between two spheres as
problem to the following problem:
depicted in Figure 2. Ideally, we want to reduce the two-sided
Pm
m
for a P
set of PSD matrices M = {Mi }i=1 such that
i=1 Mi  I, find a sparse representation
∆= m
α
M
with
small
nnz(α)
such
that
∆

I.
However,
in the reduction we use such mai=1Pi i
m
trix ∆ = i=1 αi Mi to update Aj and we need the length of Ellip(∆) is large on the direction
that Ellip(Aj ) is small. To encode this information, we define the generalised one-sided problem as
follows:
Definition 3.3 (One-sided Oracle). Let 0  B PI, C+  0, C−  0 be symmetric mam
trices, M = {Mi }m
i=1 be a set of matrices such that
i=1 Mi = I. We call a randomised algorithm Oracle (M, B, C+ , C− ) a one-sided oracle
with
speed S ∈ (0, 1] and error ε > 0, if
Pm
Oracle (M, B, C+ , C− ) outputs a matrix ∆ = i=1 αi Mi such that
6


1. nnz(α) 6 λmin (B) · tr B −1 .

2. ∆  B and αi > 0 for all i.

3. E [C • ∆] > S ·λmin (B)·tr(C)−εS ·λmin (B)·tr(C|·| ), where C = C+ −C− , and C|·| = C+ +C− .
We show in Section 3.3 the existence of a one-sided oracle with speed S = Ω(1) and error ε = 0,
in which case the oracle only requires C as input, instead of C+ and C− . However, to construct
such an oracle efficiently an additional error is introduced, which depends on C+ + C− .
For the main algorithm Sparsify (M, ε), we maintain the matrix Aj inductively as we discussed
at the beginning of Section 3. By employing the ⊕ operator, we reduce the problem of constructing
∆j with two-sided constraints to the problem of constructing ∆j ⊕ ∆j with one-sided constraints.
We also use C to ensure that the length of Ellip(∆) is large on the direction where the length of
Ellip(Aj ) is small. See Algorithm 1 for formal description.
Algorithm 1 Sparsify (M, ε)
j = 0, A0 = 0;
ℓ0 = − 14 , u0 = 14 ;
while uj − ℓj < 1 do
Bj = (uj I − Aj )2 ⊕ (Aj − ℓj I)2 ;
C+ = (1 − 2ε)(Aj − ℓj I)−2 exp (Aj − ℓj I)−1 ;
C− = (1 + 2ε)(uj I − Aj )−2 exp(uj I − Aj )−1 ;

1
1
∆j ⊕ ∆j = Oracle {Mi ⊕ Mi }m
i=1 , Bj , 2 (C+ ⊕ C+ ) , 2 (C− ⊕ C− ) ;
Aj+1 ← Aj + ε · ∆j ;
· S · λmin (Bj ) and δℓ,j = ε · (1−2ε)(1−ε)
· S · λmin (Bj );
δu,j = ε · (1+2ε)(1+ε)
1−4ε
1+4ε
uj+1 ← uj + δu,j , ℓj+1 ← ℓj + δℓ,j ;
j ← j + 1.
end while
Return Aj .

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:

To analyse Algorithm 1, we use the fact that the returned ∆j satisfies the preconditions of
Lemma 3.1 and prove in Lemma 4.4 that


E Φuj+1 ,ℓj+1 (Aj+1 ) 6 Φuj ,ℓj (Aj )

for any iteration j. Hence, with high probability the bounded ratio between uT and ℓT after the
final iteration T implies that the Ellip(AT ) is close to be a sphere. In particular,
 2 for any ε < 1/20,
n
times, which is
a (1 + O(ε))-spectral sparsfier can be constructed by calling Oracle O log
ε2 ·S
described in the lemma below.
Lemma 3.4. Let 0 < ε < 1/20. Suppose we have one-sided oracle Oracle with speed S and error
ε. Then, with constant probability the algorithm Sparsify
 2  (M, ε) outputs a (1 + O(ε))-spectral

n
n
sparsifier with O ε2 ·S vectors by calling Oracle O log
times.
ε2 ·S

3.3

Solving Oracle via SDP

Now we show that the required solution of Oracle(M, B, C) indeed exists3 , and can be further
solved in nearly-linear time by a semi-definite program. We first prove that the required matrix
3

As the goal here is to prove the existence of Oracle with error ε = 0, the input here is C instead of C+ and C− .

7

satisfying the conditions of Definition 3.3 exists for some absolute constant S = Ω(1) and ε = 0. To
make a parallel discussion between Algorithm 1 and the algorithm we will present later, we use A
to denote the output of the Oracle instead of ∆. We adopt the ideas between the ellipsoid and two
spheres discussed before, but only consider one sphere for the one-sided case. Hence, we introduce
a barrier value uj for each iteration j, where u0 = 0. We will use the potential function
Ψj = tr (uj B − Aj )−1
in our analysis, where uj is increased by δj , (Ψj · λmin (B))−1
Pafter iteration j. Moreover, since we
only need to prove the existence of the required matrix A = m
i=1 αi Mi , our process proceeds for T
iterations, where only one vector is chosen in each iteration. To find a desired vector, we perform
random sampling, where each matrix Mi is sampled with probability prob(Mi ) proportional to
Mi • C, i.e.,
(Mi • C)+
prob(Mi ) , Pm
(3.5)
+,
t=1 (Mt • C)

where x+ , max{x, 0}. Notice that, since our goal is to construct A such that E[C • A] is lower
bounded by some threshold as stated in Definition 3.3, we should not pick any matrix Mi with
Mi • C < 0. This random sampling procedure is described in Algorithm 2, and the properties of
the output matrix is summarised in Lemma 3.5.
Algorithm 2 SolutionExistence (M, B, C)


1: A0 = 0, u0 = 1 and T = λmin (B) · tr(B −1 ) ;
2: for j = 0, 1, . . . , T − 1 do
3:
repeat
4:
Sample a matrix Mt with probability prob(Mt );
5:
Let ∆j = (4Ψj · prob(Mt ))−1 · Mt ;
6:
until ∆j  12 (uj B − Aj );
7:
Aj+1 = Aj + ∆j ;
8:
δj = (Ψj · λmin (B))−1 ;
9:
uj+1 = uj + δj ;
10: end for
11: Return u1 AT .
T
Lemma 3.5. Let 0  BP I and C be symmetric matrices, and M = {Mi }m
i=1 be a set of
m
PSD P
matrices such that
i=1 Mi = I. Then SolutionExistence (M, B, C) outputs a matrix
A= m
α
M
such
that
the
following holds:
i
i
i=1


1. nnz(α) = λmin (B) · tr(B −1 ) .
2. A  B, and αi > 0 for all i.

3. E [C • A] >

1
32

· λmin (B) · tr (C).

Lemma 3.5 shows that the required matrix A defined in Definition 3.3 exists, and can be found
by random sampling described in Algorithm 2. Our key observation is that such matrix A can be
constructed by a semi-definite program.

8

TheoremP
3.6. Let 0  B  I, C be symmetric matrices, and M = {Mi }m
i=1 be a set of matrices
m
−1
such that i=1 Mi = I. Let S ⊆ [m] be a random set of λmin (B) · tr(B ) coordinates, where
every index i is picked with probability prob(Mi ). Let A⋆ be the solution of the following semidefinite
program
!
X
X
αi Mi  B
(3.6)
αi Mi subject to A =
max C •
αi >0

Then, we have E [C • A⋆ ] >

i∈S

i∈S

1
32

· λmin (B) · tr(C).

Taking the SDP formuation (3.6) and the specific constraints of the Oracle’s input into account,
the next lemma shows that the required matrix used in each iteration of Sparsify (M, ε) can be
computed efficiently by solving a semidefinite program.
Lemma 3.7. The Oracle used in Algorithm Sparsify (M, ε) can be implemented in




e (Z + nω ) · ε−O(1)
e ε−O(1) depth
O
work and O

P
Pm
where Z = m
i=1 nnz(Mi ) is the total number of non-zeros in Mi . When the matrix
i=1 Mi = I
comes from spectral sparsification of graphs, each iteration of Sparsify (M, ε) can be implemented
in




e ε−O(1) depth.
e mε−O(1) work and O
O
Furthermore, the speed of this one-sided oracle is Ω(1) and the error of this one-sided oracle is ε.
Combining Lemma 3.4 and Lemma 3.7 gives us the proof of the main result.
Proof of Theorem 1.1 and Theorem 1.2. Lemma 3.7 shows that we can construct Oracle with Ω(1)
speed and ε error that runs in




−O(1)
ω
−O(1)
e
e
depth
work and O ε
O (Z + n ) · ε

for the matrix setting and



e mε−O(1) work
O

and



e ε−O(1) depth.
O

for the graph
setting. Combining this with Lemma 3.4, which states that it suffices to call Oracle

e 1/ε2 times, the main statements hold.
O


3.4

Further discussion

Before presenting a more detailed analysis of our algorithm, we compare our new approach with
the previous ones for constructing a linear-sized spectral sparsifier, and see how we address the
bottlenecks faced in previous constructions. Notice that all previous algorithms require super
poly-logarithmic number of iterations, and super linear-time for each iteration. For instance, our
previous algorithm [LS15] for constructing a sparsifier with O(pn) edges requires Ω(n1/p ) iterations
and Ω(n1+1/p ) time per iteration for the following reasons:

• Ω n1+1/p time is needed per iteration: Each iteration takes n/gΩ(1) time to pick the vector(s)
when (ℓ + g)I  A  (u − g)I. To avoid eigenvalues of A getting too close to the boundary
u or ℓ, i.e., g being too small, we choose the potential function whose value dramatically
increases when the eigenvalues of A get close u or ℓ. As the cost, we need to scale down the
added vectors by a n1/p factor.
9



• Ω n1/p iterations are needed: By random sampling, we choose O n1−1/p vectors each

iteration and use the matrix Chernoff bound to show that the “quality” of added O n1−1/p

vectors is just p = Θ(1) times worse than adding a single vector. Hence, this requires Ω n1/p
iterations.
In contrast, our new approach breaks these two barriers through the following way:
• A “non-linear” step: Instead of rescaling down the vectors we add uniformly, we pick much
fewer vectors on the direction that blows up, i.e., we impose the condition ∆  (uI − A)2
instead of ∆ 
 1/p · (uI − A). This allows us to use the new potential function (3.3) with
−1
form exp x
to control the eigenvalues in a more aggressive way.

• SDP filtering: By matrix Chernoof bound, we know that the probability that we sample a
few “bad” vectors is small. Informally, weapply semi-definite
programming to filter out those

O(1)
bad vectors, and this allows us to add Ω n/ log
(n) vectors in each iteration.

4

Detailed analysis

In this section we give detailed analysis for the statements presented in Section 3.

4.1

Analysis of the potential function

Now we analyse the properties of the potential function (3.3), and prove Lemma 3.1. The following
two facts from matrix analysis will be used in our analysis.
Lemma 4.1 (Woodbury Matrix Identity). Let A ∈ Rn×n , U ∈ Rn×k , C ∈ Rk×k and V ∈ Rk×n be
matrices. Suppose that A, C and C −1 + V A−1 U are invertible, it holds that
−1
(A + U CV )−1 = A−1 − A−1 U C −1 + V A−1 U
V A−1 .
Lemma 4.2 (Golden-Thompson Inequality). It holds for any symmetric matrices A and B that


tr eA+B 6 tr eA · eB .

Proof of Lemma 3.1. We analyse the change of Φu (·) and Φℓ (·) individually. First of all, notice
that
−1

(uI − A)−1/2 .
(uI − A − ∆)−1 = (uI − A)−1/2 I − (uI − A)−1/2 ∆(uI − A)−1/2

We define Π = (uI − A)−1/2 ∆(uI − A)−1/2 . Since 0  ∆  δ(uI − A)2 and u − ℓ 6 1, it holds that
Π  δ(uI − A)  δ(uI − ℓI)  δI,
and therefore
(I − Π)−1  I +

1
· Π.
1−δ

Hence, it holds that
−1

(uI − A − ∆)


1
 (uI − A)
· Π (uI − A)−1/2
I+
1−δ
1
· (uI − A)−1 ∆(uI − A)−1 .
= (uI − A)−1 +
1−δ
−1/2



10

By the fact that tr exp is monotone and the Golden-Thompson Inequality (Lemma 4.2), we have
that

Φu (A + ∆) = tr exp (uI − A − ∆)−1


1
−1
−1
−1
· (uI − A) ∆(uI − A)
6 tr exp (uI − A) +
1−δ



1
−1
−1
−1
6 tr exp(uI − A) exp
· (uI − A) ∆(uI − A)
.
1−δ
Since 0  ∆  δ(uI − A)2 and δ 6 1/10 by assumption, we have that (uI − A)−1 ∆(uI − A)−1  δI,
and


1
−1
−1
exp
· (uI − A) ∆(uI − A)
 I + (1 + 2δ) · (uI − A)−1 ∆(uI − A)−1 .
1−δ
Hence, it holds that


−1
Φu (A + ∆) 6 tr e(uI−A) · I + (1 + 2δ)(uI − A)−1 ∆(uI − A)−1
= Φu (A) + (1 + 2δ) · tr(e(uI−A) (uI − A)−2 ∆).
−1

By the same analysis, we have that


−1
Φℓ (A + ∆) 6 tr e(A−ℓI) · I − (1 − 2δ)(A − ℓI)−1 ∆(A − ℓI)−1
= Φℓ (A) − (1 − 2δ) · tr(e(A−ℓI) (A − ℓI)−2 ∆).
−1

Combining the analysis on Φu (A + ∆) and Φℓ (A + ∆) finishes the proof.



Lemma 4.3. Let A be a symmetric matrix. Let u, ℓ be the barrier values such that u − ℓ 6 1 and
ℓI ≺ A ≺ uI. Assume that 0 6 δu 6 δ · λmin (uI − A)2 and 0 6 δℓ 6 δ · λmin (A − ℓI)2 for δ 6 1/10.
Then, it holds that


−1
Φu+δu ,ℓ+δℓ (A) 6Φu,ℓ (A) − (1 − 2δ)δu · tr e(uI−A) (uI − A)−2


−1
+ (1 + 2δ)δℓ · tr e(A−ℓI) (A − ℓI)−2 .
Proof. Since 0 6 δu 6 δ·λmin (uI −A)2 and 0 6 δℓ 6 δ·λmin (A−ℓI)2 , we have that δu ·I  δ·(uI −A)2
and δℓ · I  δ · (A − ℓI)2 . The statement follows by a similar analysis for proving Lemma 3.1. 

4.2

Analysis of the reduction

Now we present the detailed analysis for the reduction from a spectral sparsifier to a one-sided
oracle. We first analyse Algorithm 1, and prove that in expectation the value of the potential
function is not increasing. Based on this fact, we will give a proof of Lemma
 2 3.4,
 which shows that
log n
a (1 + O(ε))-spectral sparsifier can be constructed by calling Oracle O ε2 ·S times.

Lemma 4.4. Let Aj and Aj+1 be the matrices constructed by Algorithm 1 in iteration j and j + 1,
and assume that 0 6 ε 6 1/20. Then, it holds that


E Φuj+1 ,ℓj+1 (Aj+1 ) 6 Φuj ,ℓj (Aj ).
11

Proof. By the description of Algorithm 1 and Definition 3.3, it holds that
∆j ⊕ ∆j  (uj I − Aj )2 ⊕ (Aj − ℓj I)2 ,
which implies that ∆j  (uj I − Aj )2 and ∆j  (Aj − ℓj I)2 . Since uj − ℓj 6 1 by the algorithm
description and 0 6 ε 6 1/20, by setting ∆ = ε · ∆j in Lemma 3.1, we have


−1
Φuj ,ℓj (Aj + ε · ∆j ) 6 Φuj ,ℓj (Aj ) + ε(1 + 2ε) · tr e(uj I−Aj ) (uj I − Aj )−2 ∆j


−1
− ε(1 − 2ε) · tr e(Aj −ℓj I) (Aj − ℓj I)−2 ∆j
= Φuj ,ℓj (Aj ) − ε · C • ∆j .

Notice that the matrices {Mi ⊕ Mi }m
i=1 as the input of Oracle always satisfy
Using this and the definition of Oracle, we know that

Pm

i=1 Mi ⊕ Mi

= I ⊕ I.


E [C • ∆j ] > S · λmin (Bj ) · tr(C) − ε · S · λmin (Bj ) · tr C|·| ,

Let αj = ε · S · λmin (Bj ). Then we have that


E Φuj ,ℓj (Aj+1 ) 6 Φuj ,ℓj (Aj ) − ε · E [C • ∆j ]



−1
6 Φuj ,ℓj (Aj ) + (1 + 2ε)(1 + ε) · αj · tr e(uj I−Aj ) (uj I − Aj )−2


−1
− (1 − 2ε)(1 − ε) · αj · tr e(Aj −ℓj I) (Aj − ℓj I)−2 .

(4.1)

On the other hand, using that 0 6 ε 6 1/20, S 6 1 and ∆j  (uj I − Aj ), we have that
δu,j 6 ε ·

(1 + 2ε)(1 + ε)
· λmin (uj I − Aj )2 6 2ε · λmin (uj I − Aj+1 )2
1 − 4ε

and

(1 − 2ε)(1 − ε)
· λmin (Aj − ℓj I)2 6 2ε · λmin (Aj+1 − ℓj I)2 .
1 + 4ε
Hence, Lemma 4.3 shows that


−1
Φuj +δu,j ,ℓj +δℓ,j (Aj+1 ) 6 Φuj ,ℓj (Aj+1 ) − (1 − 4ε)δu,j · tr e(uj I−Aj+1 ) (uj I − Aj+1 )−2


−1
+ (1 + 4ε)δℓ,j · tr e(Aj+1 −ℓj I) (Aj+1 − ℓj I)−2


−1
6 Φuj ,ℓj (Aj+1 ) − (1 − 4ε)δu,j · tr e(uj I−Aj ) (uj I − Aj )−2


−1
(4.2)
+ (1 + 4ε)δℓ,j · tr e(Aj −ℓj I) (Aj − ℓj I)−2 .
δℓ,j 6 ε ·

By combining (4.1),
(4.2), and setting
(1−4ε)δu,j = (1+2ε)(1+ε)αj , (1+4ε)δℓ,j = (1−2ε)(1−ε)αj ,


we have that E Φuj+1 ,ℓj+1 (Aj+1 ) 6 Φuj ,ℓj (Aj ).

Proof of Lemma 3.4. We
bound the number of times the algorithm calls the oracle. Notice
 first


 
1 −1
that Φu0 ,ℓ0 = 2 · tr exp 4 I
= 2e4 · n. Hence, by Lemma 4.4 we have E Φuj ,ℓj (Aj ) = O(n) for

any iteration j. By Markov’s inequality, it holds that Φuj ,ℓj (Aj ) = nO(1) with high probability in
n. In the remainder of the proof, we assume that this event occurs.

12

Since Bj = (uj I − Aj )2 ⊕ (Aj − ℓj I)2 by definition, it holds that


exp (λmin (Bj ))−1/2 6 Φuj ,ℓj (Aj ) = nO(1) ,

which implies that


λmin (Bj ) = Ω log−2 n .

On the other hand, in iteration j the gap between uj and ℓj is increased by

δu,j − δℓ,j = Ω ε2 · S · λmin (Bj ) .

(4.3)

(4.4)

Combining this with (4.3) gives us that

δu,j − δℓ,j = Ω



ε2 · S
log2 n



for any j. Since u0 − ℓ0 = 1/2 and the algorithmterminates
once uj − ℓj > 1 for some j, with high

log2 n
probability in n, the algorithm terminates in O ε2 ·S iterations.

Next we prove that the number of Mi ’s involved in the output is at most O ε2n·S . By the
properties of Oracle, the number of matrices in iteration j is at most λmin (Bj ) · tr(Bj−1 ). Since

x−2 6 exp x−1 for all x > 0, it holds for any iteration j that






tr Bj−1 = tr (uj I − Aj )−2 + tr (Aj − ℓj I)−2 6 Φuj ,ℓj (Aj )

By (4.4), we know that for 
added matrix
 Mi in iteration j, the average increase of the gap uj − ℓj


2
for each added matrix is Ω Φu ε,ℓ ·S(Aj ) . Since E Φuj ,ℓj (Aj ) = O(n), for every new added matrix,
j j
 2 
in expectation the gap between uj and ℓj is increased by Ω ε n·S . By the ending condition of the
algorithm, i.e.,
 uj − ℓj > 1, and Markov’s inequality, the number of matrices picked in total is at
most O ε2n·S with constant probability.
Finally we prove that the output is a (1 + O(ε))-spectral sparsifier. Since the condition number
of the output matrix Aj is at most
uj
=
ℓj



uj − ℓj
1−
uj

−1

,

it suffices to prove that (uj − ℓj )/uj = O(ε) and this easily follows from the ending condition of
the algorithm and
δu,j − δℓ,j
= O(ε).
δu,j


4.3

Existence proof for Oracle

Proof of Lemma 3.5. The property on nnz (α) follows from the algorithm description. For the
second property, notice that every chosen matrix ∆j in iteration j satisfies ∆j  21 (uj B − Aj ),
which implies that Aj  uj B holds for any iteration j. Hence, A = u1T AT  B, and αi > 0 since
Ψj > 0.
13

Now we prove the third statement. Let
β=

m
X

(Mi • C)+ .

i=1

Then, for each matrix Mij picked in iteration j, C • Aj is increased by
β
1
· C • Mij =
.
4Ψj · prob(Mij )
4Ψj

C • ∆j =
On the other hand, it holds that

T
−1
X

uT = u0 +

δj = 1 +

T
−1
X

(Ψj · λmin (B))−1 .

j=0

j=0

Hence, we have that


1
C •A=
·C•
uT

T
−1
X

βλmin (B)
·
4

>

βλmin (B)
·
4

>

βλmin (B)
·
4

∆j  =

PT −1

−1
j=0 β · (4Ψj )
PT −1
−1
j=0 (Ψj · λmin (B))

1+
PT −1 −1
j=0 Ψj
P −1 −1
λmin (B) + Tj=0
Ψj
PT −1
−1
j=0 (Ψj + Ψ0 )
PT −1
λmin (B) + j=0 (Ψj + Ψ0 )−1
PT −1
T −1
−1
β X
j=0 (Ψj + Ψ0 )
> ·
(Ψj + Ψ0 )−1 ,
8
λmin (B) + T · Ψ−1
0

j=0

=



(4.5)

j=0

where the last inequality follows by the choice of T . Hence, it suffices to bound Ψj .
Since ∆j  12 (uj B − Aj )  12 (uj+1 B − Aj ), we have that

Ψj+1 6 tr (uj+1 B − Aj )−1 + 2 · ∆j • (uj+1 B − Aj )−2 .
Since tr(uB − Aj )−1 is convex in u, we have that



tr (uj B − Aj )−1 > tr (uj+1 B − Aj )−1 + δj tr (uj+1 B − Aj )−2 B

Combining (4.6) and (4.7), we have that




Ψj+1 6 tr (uj B − Aj )−1 − δj · tr (uj+1 B − Aj )−2 B + 2 · ∆j • (uj+1 B − Aj )−2

= Ψj − δj · λmin (B) · tr (uj+1 B − Aj )−2 + 2 · ∆j • (uj+1 B − Aj )−2

(4.6)

(4.7)

(4.8)

Let Ej be the event that ∆j  12 (uj B − Aj ). Notice that our picked ∆j in each iteration always
satisfies Ej by algorithm description. Since


E ∆j • (uj B − Aj )−1 =

X

prob(Mi ) ·

i:(Mi •C)+ >0

14

1
1
· Mi • (uj B − Aj )−1 6 ,
4Ψj · prob(Mi )
4

by Markov inequality it holds that

1
1
P [Ej ] = P ∆j  (uj B − Aj ) > ,
2
2


and therefore


 E ∆j • (uj+1 B − Aj )−2

−2
E ∆j • (uj+1 B − Aj ) | Ej 6
P (Ej )


6 2 · E ∆j • (uj+1 B − Aj )−2
X
1
=
·
Mi • (uj+1 B − Aj )−2
2Ψj
+
i:(Mi •C) >0

1
6
· tr(uj+1 B − Aj )−2 .
2Ψj

Combining the inequality above, (4.8), and the fact that every ∆j picked by the algorithm satisfies
E, we have that


1
− δj · λmin (B) · tr (uj+1 B − Aj )−2 .
E [Ψj+1 ] 6 Ψj +
Ψj
By our choice of δj , it holds for any iteration j that E [Ψj+1 ] 6 Ψj , and
i
h
E (Ψj+1 + Ψ0 )−1 > E (Ψj + Ψ0 )−1 >

1
.
2 · Ψ0

Combining this with (4.5), it holds that
T −1

βX 
T
tr(C)
T
β T
β
E [C • A] >
E (Ψj + Ψ0 )−1 >
·
·
>
·
.
=
−1
8
16 Ψ0
16 tr (B )
16
tr (B −1 )
j=0


The result follows from the fact that T > λmin (B)tr B −1 /2.



Using the lemma above, we can prove that such A can be solved by a semidefinite program.

Proof of Theorem 3.6. Note that the probability we used in the statement is the same as the probability we used in SolutionExistence
(M, B, C). Therefore, Lemma 3.5 shows that there is a
P
matrix A of the form m
α
M
such
that
i
i
i=1
E [C • A] >

1
· λmin (B) · tr (C) ,
32

The statement follows by the fact that A⋆ is the solution of the semidefinite program (3.6) that
maximises C • A.


4.4

Implementing the SDP in nearly-linear time

Now, we discuss how to solve the SDP (3.6) in nearly-linear time. Since this SDP is a packing SDP,
it is known how to solve it in nearly-constant depth [JY11, ALO16, PT12]. The following result
will be used in our analysis.

15

Theorem 4.5 ([ALO16]). Given a SDP
max c⊺ x subject to
x>0

m
X

xi Ai  B

i=1

with Ai  0, B  0 and c ∈ Rm . Suppose that we are given a direct access to the vector c ∈ Rm
and an indirect access to Ai and B via an oracle OL,δ which inputs a vector x ∈ Rm and outputs a
vector v ∈ Rm such that
#
!
!
"

X
δ
xi Ai − B B −1/2 B −1/2
Ai • B −1/2 exp L · B −1/2
vi ∈ 1 ±
2
i

in WL,δ work and DL,δ depth for any x such that xi > 0 and
x such that
⊺

E [c x] > (1 − O(δ))OPT

with

Pm

i=1 xi Ai

m
X

 2B. Then, we can output

xi Ai  B

i=1

in


O WL,δ log m · log (nm/δ) /δ3 work

and

O (DL,δ log m · log(nm/δ)/δ) depth,

where L = (4/δ) · log(nm/δ).

Since we are only interested in a fast implementation of the one-sided oracle used in Sparsify (M, ε),
it suffices to solve the SDP (3.6) for this particular situation.
Proof of Lemma 3.7. Our basic idea is to use Theorem 4.5 as the one-sided oracle. Notice that
each iteration of Sparsify (M, ε) uses the one-sided oracle with the input

1 − 2ε
(A − ℓI)−2 exp(A − ℓI)−1 ⊕ (A − ℓI)−2 exp(A − ℓI)−1 ,
2

1 + 2ε
C− =
(uI − A)−2 exp(uI − A)−1 ⊕ (uI − A)−2 exp(uI − A)−1 ,
2
B = (uI − A)2 ⊕ (A − ℓI)2 ,

C+ =

where we drop the subscript j indicating iterations here for simplicity. To apply Theorem 3.6, we
first sample a subset S ⊆ [m], then solve the SDP
!
m
m
X
X
βi Mi ⊕ Mi  B.
subject to
βi Mi ⊕ Mi
max
(C+ − C− ) •
βi >0,βi =0 on i∈S
/

i=1

i=1

By ignoring the matrices Mi with i 6∈ S, this SDP is equivalent to the SDP
max c⊤ β
βi >0

subject to

m
X

βi Mi ⊕ Mi  B

i=1

where ci = (C+ − C− ) • (Mi ⊕ Mi ). Now assume that (i)
can approximate ci with δ(C+ + C− ) •
Pwe
m
e
(Mi ⊕ Mi ) additive error, and (ii) for any x such that i=1 xi Mi ⊕ Mi  2B and L = O(1/δ),
we
can approximate
!
!
X
−1/2
−1/2
−1/2
B −1/2
xi Mi ⊕ Mi − B B
(Mi ⊕ Mi ) • B
exp L · B
i

16

with 1 ± δ multiplicative error. Then, by Theorem 4.5 we can find a vector β such that
⊺

E [c β] > (1 − O(δ))OPT − δ
> (1 − O(δ))OPT − δ

m
X

i=1
m
X

βi (C+ + C− ) • (Mi ⊕ Mi )
βi (C+ + C− ) • B

i=1

> OPT − O(δ)(C+ + C− ) • B.
where we used that Mi ⊕ Mi  B and OPT 6 (C+ + C− ) • B. Since u − ℓ 6 1, we have that
B  I ⊕ I and hence
E [c⊺ β] > OPT − O(δ) · tr(C+ + C− )

1
λmin (B) · tr(C) − O(δ log2 n) · tr(C+ + C− )
>
32

where we apply Theorem 3.6 and (4.3) at the last line. Therefore, this gives an oracle with speed
1/32 and ε error by setting δ = ε/ log2 n.
The problem of approximating sample probabilities, {ci }, as well as implementing the oracle
OL,δ is similar with approximating leverage scores [SS11], and relative leverage scores [ALO15,
LS15]. All these references use the Johnson-Lindenstrauss lemma to reduce the problem of approximating matrix dot product or trace to matrix vector multiplication. The only difference is
that, instead of computing (A − ℓI)−(q+1) x and (uI − A)−(q+1) x for a given vector x in other references, we compute (A − ℓI)−2 exp(A − ℓI)−1 x and (uI − A)−2 exp(uI − A)−1 x. These can be
approximated by Taylor expansion and the number of terms required for Taylor expansion depends
on how close the eigenvalues of A are to the boundary (u or ℓ). In particular, we show in Sec2 ) terms in Taylor expansion suffices, where the gap g is the largest number
e
tion 4.5 that O(1/g
e
such that (ℓ + g)I  A  (u− g)I. Since 1/g2 = O(1)
by (4.3), each iteration can be impleO(1)
e 1/ε
e 1/εO(1) matrix vector multiplication. For
mented via solving O
linear systems and O
the matrices coming from graph sparsification, this can be done in nearly-linear work and nearlyconstant depth [PS14, KLP+ 16]. For general matrices, this can be done in input sparsity time and
nearly-constant depth [NN13, LMP13, CLM+ 15].


4.5

Taylor Expansion of x−2 exp(x−1 )

Theorem 4.6 (Cauchy’s Estimates). Suppose f is holomorphic on a neighborhood of the ball
B , {z ∈ C : |z − s| 6 r}, then we have that
f (k) (s) 6

k!
sup |f (z)| .
r k z∈B

Lemma 4.7. Let f (x) = x−2 exp(x−1 ). For any 0 < x 6 1, we have that
f (x) −

d
X
5
1 (k)
f (1)(x − 1)k 6 8(d + 1)e x −xd .
k!
k=0

In particular, if d >

c
x2

1
) for some large enough universal constant c, we have that
log( xε

f (x) −

d
X
1 (k)
f (1)(x − 1)k 6 ε.
k!
k=0

17

Proof. By the formula of the remainder term in Taylor series, we have that
Z
d
X
1 x (d+1)
1 (k)
k
f (1)(x − 1) +
f
(s)(x − s)d ds.
f (x) =
k!
d! 1
k=0

For any s ∈ [x, 1], we define D(s) = {z ∈ C : |z − s| 6 s − x2 }. Since |f (z)| 6 (x/2)−2 exp(2/x) on
z ∈ D(s), Cauchy’s estimates (Theorem 4.6) shows that
 
2
(d + 1)! 4
(d + 1)!
(d+1)
.
f
(s) 6
x d+1 sup |f (z)| 6
x d+1 2 exp
(s − 2 )
(s − 2 )
x
x
z∈B(s)
Hence, we have that
f (x) −

 
Z x
d
X
1 (k)
(d + 1)! 4
2
1
f (t)(x − t)k 6
exp
(x − s)d ds
k!
d! 1 (s − x2 )d+1 x2
x
k=0
2 Z
4(d + 1)e x 1 (s − x)d
=
x d+1 ds
x2
x (s − 2 )
2 Z
8(d + 1)e x 1
(1 − x)d ds
6
x3
x
5

68(d + 1) · e x −xd .



Acknowledgement
The authors would like to thank Michael Cohen for helpful discussions and suggesting the ideas to
improve the sparsity from O(n/ε3 ) to O(n/ε2 ).

References
[ALO15] Zeyuan Allen-Zhu, Zhenyu Liao, and Lorenzo Orecchia. Spectral sparsification and
regret minimization beyond multiplicative updates. In 47th Annual ACM Symposium
on Theory of Computing (STOC’15), pages 237–245, 2015.
[ALO16] Zeyuan Allen-Zhu, Yin Tat Lee, and Lorenzo Orecchia. Using optimization to obtain
a width-independent, parallel, simpler, and faster positive SDP solver. In 27th Annual
ACM-SIAM Symposium on Discrete Algorithms (SODA’16), pages 1824–1831, 2016.
[AT11] Haim Avron and Sivan Toledo. Randomized algorithms for estimating the trace of an
implicit symmetric positive semi-definite matrix. Journal of the ACM (JACM), 58(2):8,
2011.
e 2 ) time. In
[BK96] András Benczúr and David Karger. Approximating s-t minimum cuts in O(n
28th Annual ACM Symposium on Theory of Computing (STOC’96), pages 47–55, 1996.

[BSS12] Joshua D. Batson, Daniel A. Spielman, and Nikhil Srivastava. Twice-Ramanujan sparsifiers. SIAM Journal on Computing, 41(6):1704–1721, 2012.
18

[BSST13] Joshua D. Batson, Daniel A. Spielman, Nikhil Srivastava, and Shang-Hua Teng. Spectral
sparsification of graphs: theory and algorithms. Communications of the ACM, 56(8):87–
94, 2013.
[Che89] L. Paul Chew. There are planar graphs almost as good as the complete graph. Journal
of Computer and System Sciences, 39(2):205–219, 1989.
[CLM+ 15] Michael B Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng,
and Aaron Sidford. Uniform sampling for matrix approximation. In 6th Innovations in
Theoretical Computer Science (ITCS’15), pages 181–190, 2015.
[JY11] Rahul Jain and Penghui Yao. A parallel approximation algorithm for positive semidefinite programming. In 52th Annual IEEE Symposium on Foundations of Computer
Science (FOCS’11), pages 463–471, 2011.
[KLP+ 16] Rasmus Kyng, Yin Tat Lee, Richard Peng, Sushant Sachdeva, and Daniel A. Spielman.
Sparsified Cholesky and multigrid solvers for connection Laplacians. In 48th Annual
ACM Symposium on Theory of Computing (STOC’16), pages 842–850, 2016.
[LMP13] Mu Li, Gary L Miller, and Richard Peng. Iterative row sampling. In 54th Annual IEEE
Symposium on Foundations of Computer Science (FOCS’13), pages 127–136, 2013.
[LS15] Yin Tat Lee and He Sun. Constructing linear-sized spectral sparsification in almostlinear time. In 56th Annual IEEE Symposium on Foundations of Computer Science
(FOCS’15), pages 250–269, 2015.
[NN13] Jelani Nelson and Huy L Nguyên. Osnap: Faster numerical linear algebra algorithms
via sparser subspace embeddings. In 54th Annual IEEE Symposium on Foundations of
Computer Science (FOCS’13), pages 117–126, 2013.
[PS14] Richard Peng and Daniel A. Spielman. An efficient parallel solver for SDD linear systems.
In 46th Annual ACM Symposium on Theory of Computing (STOC’14), pages 333–342,
2014.
[PT12] Richard Peng and Kanat Tangwongsan. Faster and simpler width-independent parallel algorithms for positive semidefinite programming. In 24th ACM Symposium on
Parallelism in Algorithms and Architectures (SPAA’12), pages 101–108, 2012.
[SS11] Daniel A. Spielman and Nikhil Srivastava. Graph sparsification by effective resistances.
SIAM Journal on Computing, 40(6):1913–1926, 2011.
[ST11] Daniel A. Spielman and Shang-Hua Teng. Spectral sparsification of graphs. SIAM
Journal on Computing, 40(4):981–1025, 2011.
[Zou12] Anastasios Zouzias. A matrix hyperbolic cosine algorithm and applications. In 39th International Colloquium on Automata, Languages, and Programming (ICALP’12), pages
846–858, 2012.

19

