arXiv:1507.07497v4 [cs.DM] 21 Apr 2016

An Efficient Parallel Algorithm for Spectral Sparsification of
Laplacian and SDDM Matrix Polynomials
Gorav Jindal
Pavel Kolev∗
Max-Planck-Institut für Informatik, Saarbrücken, Germany
{gjindal,pkolev}@mpi-inf.mpg.de

Abstract
A mixture of discrete Binomial distributions (MDBD), denoted by BN,T (α, p), is a set of pairs
{(B(pi , N ), αi )}Ti=1 , where B(·, ·) denotes the Binomial distribution, all pi ∈ (0, 1) are distinct,
PT
PT
i=1 αi 6 1 and all αi ∈ (0, 1). A vector γ is induced by MDBD if γi =
j=1 αj · BN,i (pj ) for
N i
N −i
all i ∈ [0 : N ], where BN,i (p) = i p (1 − p)
.
We prove for “large” class C of continuouspprobability density functions (p.d.f.), that for
every w ∈ C there exists MDBD with T > N φw /δ that δ-approximates a discretized p.d.f.
PN
w(i/N
b
) , w(i/N )/[ ℓ=0 w(ℓ/N )] for all i ∈ [3 : N −3], where φw > maxx∈[0,1] |w(x)|. Moreover,
we propose an efficient parallel algorithm that on input p.d.f. w ∈ C and parameter δ > 0,
outputs MDBD that induces a vector γ which δ-approximates w.
b Also, we give an efficient
parallel algorithm that on input a discretized p.d.f. w
b induced by MDBD with T = N + 1 and
the corresponding vector p, outputs exactly the coefficients α.
Cheng et al. [4] proposed the first sequential algorithm that on input a discretized p.d.f. β,
B = D − M that is either Laplacian or SDDM matrix and parameter ε ∈ (0, 1), outputs in time
cN ≈ε D − D PN βi (D−1 M )i .
b −2 mN 2 )1 a spectral sparsifier of a matrix-polynomial D − M
O(ε
i=0
However, given MDBD BN,T (α, p) that induces a discretized p.d.f. γ, to apply the algorithm
in [4] one has to explicitly precompute in O(N T ) time the vector γ. Instead, we give two
algorithms (sequential and parallel) that bypass this explicit precomputation.
We propose a faster sequential algorithm that on input MDBD BN,T (α, p) with N = 2k
b −2 m + ε−4 nT ) time the desired spectral sparsifier. Moreover, our
for k ∈ N+ outputs in O(ε
b −2 m + ε−4 nT ) work and O(log N · poly(log n) + log T )
algorithm is parallelizable and runs in O(ε
depth. Our main algorithmic contribution is to propose the first efficient parallel algorithm that
on input continuous p.d.f. w ∈ C, matrix B = D − M as above, outputs a spectral sparsifier of
matrix-polynomial whose coefficients approximate component-wise the discretized p.d.f. w.
b
Our results yield the first efficient and parallel algorithm that runs in nearly linear work
and poly-logarithmic depth and analyzes the long term behaviour of Markov chains in nontrivial settings. In addition, we strengthen the Spielman and Peng’s [21] parallel SDD solver by
introducing a simple parallel preprocessing step.

∗
This work has been funded by the Cluster of Excellence “Multimodal Computing and Interaction” within the
Excellence Initiative of the German Federal Government.
1b
O(·) notation hides poly(log n, log N ) factors.

Contents
1 Introduction

2

2 Our
2.1
2.2
2.3
2.4

2
3
4
6
6

Results
Representational Power of MDBD . . . . . . . . . . . . . . . . .
Spectral Sparsification of Matrix Polynomials induced by MDBD
Analyzing the Long Term Behaviour of Markov Chains . . . . . .
Faster SDDM Solver . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

3 Algorithmic Background on Spectral Sparsification
3.1 Prior Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Spectral Sparsification of T -Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
8

4 Core Iterative Algorithm
9
4.1 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.2 Iterative Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
5 Spectral Sparsification of Binomial T -Matrix Polynomials

11

6 Spectral Sparsification of T -Matrix Polynomials Induced by MDBD

11

7 Parallelization of Algorithm SS MDBD

12

8 Faster SDDM Solver

14

9 Representational Power of MDBD

14

10 Approximating Discretized PDF

16

11 Efficient Parallel Solver for Transpose Bernstein-Vandermonde Systems

18

A Generalized Escaping Probability

21

B Spectral Sparsification of T -Matrices
21
B.1 Structural Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
C Bernstein Basis Matrix

23

D Approximating Two Canonical PDFs
23
D.1 Uniform Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
D.2 Exponential Families . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
E Schur Complement

25

1

1

Introduction

In their seminal work Spielman and Teng [26] introduced the notion of spectral sparsifiers and
proposed the first nearly linear time algorithm for spectral sparsification. In consecutive work,
Spielman and Srivastava [24] proved that spectral sparsifiers with O(ε−2 n log n) edges exist and
e logc n · log(wmax /wmin ))2 time for any undirected graph G = (V, E, w).
can be computed in O(m
The computational bottleneck of their algorithm is to approximate the solutions of logarithmically
many SDD3 systems.
Recently, Koutis, Miller and Peng [15] developed an improved solver for SDD systems that works
e log n · log(1/ε)) time. In a survey result [13, Theorem 3] Kelner and Levin showed that in
in O(m
e log2 n) time all effective resistances can be approximated up to a constant factor. This yields
O(m
a (1 ± ε)-spectral sparsifier with only a constant factor blow-up of non-zero edges O(ε−2 n log n).
Although there are faster by a poly log-factor sparsification algorithms [16] they output spectral
sparsifiers with poly log-factor more edges.
Spielman and Peng [21] introduced the notion of sparse approximate inverse chain of SDDM4
matrices. They proposed the first parallel algorithm that finds such chains and runs in work
e log3 n·log2 κ) and depth O(logc n·log κ), where κ is the condition number of the SDDM matrix
O(m
e −2 m log3 n) time a
with m non-zero entries and dimension n. Furthermore, they showed that in O(ε
e A
e ≈ε D−AD −1 A can be computed with nnz(A)
e 6 O(ε−2 n log n). In a follow
spectral sparsifier D−
up work, Cheng et al. [3] designed an algorithm that computes a sparse approximate generalized
e such that C
eC
e T ≈ε M p for any SDDM matrix M and |p| 6 1. The chain C
e is constructed
chain C
fi D −1 M
fi
]
iteratively and it involves a normalization step that produces a sparsifier D−Mi+1 ≈ε D−M
that is expressed in terms of the original diagonal matrix D, for all iterations i.
Sinclair and Jerrum [23] analyzed Markov chains with transition matrices W = [I + D −1 A]/2,
corresponding to lazy random walks.
fast to stationary
P They proved that−2these walks converge
−1
5
distribution, defined by πu = du /( u∈V du ), after O(φG log(minu∈V πu )) steps. Andersen et
al. [1] gave an efficient local clustering algorithm
that relies on a lazy variation of PageRank, the
P
transition matrix of which is defined by ∞
α(1
− α)t W t , where α > 0 is a parameter. Their
t=0
local algorithm uses a truncated (finite summation) version of the preceding transition matrix.
b ≈ε
Recently,
et al. [4] initiated the study of computing spectral sparsifiers D − A
PN Cheng
D − D i=1 ξi (D −1 A)i of random walk Laplacian matrix polynomials, where ξ is a probability
PN
−1
i
distribution over [1 : N ], D − A is a Laplacian matrix and
i=1 ξi (D A) is a random walk
transition matrix. These matrix polynomials capture the long term behaviour of Markov chains.
Moreover, a sparsifier of a matrix polynomial yields a multiplicative approximation of the expected
generalized “escaping probability” [10, 14] of random walks. Cheng et al. [4] gave the first sequential algorithm that computes a spectral sparsifier of a random walk Laplacian matrix polynomial
e −2 · mN 2 · logc1 n · logc2 N ) for some small constants c1 , c2 .
and runs in time O(ε

2

Our Results

The lazy random walk length N in the regime of interest in [1, 10, 23] is of order N = Θ(poly(n)).
The quadratic runtime dependance on N makes the algorithm in [4] prohibitively expensive for
analysing the long term behaviour of Markov chains. In this paper, we overcome this issue for
e notation hides O(poly(log log n)) factors.
The O(·)
SDD is the class of symmetric and diagonally dominant matrices.
4
SDDM is the class of positive definite SDD matrices with non-positive off-diagonal entries.
P
5
Graph conductance φG , minµ(S)6µ(V )/2 φ(S), where φ(S) , |w(S, S)|/µ(S) and µ(S) , u∈S du .

2

3

2

“large” class of probability distributions γ over [0 : N ] that are induced by mixture of discrete
Binomial distributions (MDBD) with N = 2k for k ∈ N+ . Our results are summarized as follows.
In Subsection 2.1, we analyze the representational power of MDBD. In Subsection 2.2, we give
a sequential and a parallel algorithm for computing a spectral sparsifier of matrix polynomials
induced by MDBD. In Subsection 2.3, we propose the first parallel algorithm that runs in nearly
linear work and poly-logarithmic depth and analyzes the long term behaviour of Markov chains in
non-trivial settings. In Subsection 2.4, we strengthen the Spielman and Peng’s [21] parallel SDD
solver.

2.1

Representational Power of MDBD

Let B(p, N ) be Binomial distribution for some parameters p ∈ (0, 1) and N ∈ N+ . MDBD is a set
of pairs {(B(pi , N ), αi )}Ti=1 , denoted by BN,T (α, p), that satisfies the following two conditions:
1. (distinctness) pi ∈ (0, 1) and pi 6= pj for all i 6= j ∈ [1 : T ];
PT
2. (positive linear combination)
i=1 αi 6 1 and αi ∈ (0, 1) for all i ∈ [1 : T ].
We prove in Section 9 that for every function w in a “large” class of continuous p.d.f., there
exists MDBD that induces a component-wise approximation of w.
Theorem 2.1 (MDBD Yields a Component-Wise Approximation). Let w(x) be a four times differentiable p.d.f., εI > 0 a parameter and I = [0, 1] an interval. Suppose there is an integer N0 ∈ N
and reals µ ∈ (0, 1) and φw > 1 such that:
3) maxx∈I |w(x)| 6 φw ,
2) maxx∈I |w′ (x)| 6 12 φw · N0 ,
1) maxx∈I |w′′ (x)| 6 2φw · N02 ,
1
1
2
5) maxx∈I |b1 (x)| 6 2 µ · N0 ,
4) maxx∈I |b2 (x)| 6 2 µ · N0 ,
1
[−w(x) + (1 − 2x)w′ (x) + 21 x(1 − x)w′′ (x)] and
where the functions b1 , b2 are defined by b1 (x) = w(x)
1
b2 (x) = w(x)
[w(x)−3(1−2x)w′ (x)+(1−6x+6x2 )w′′ (x)+ 65 x(1−x)(1−2x)w′′′ (x)+ 81 x2 (1−x)2 w′v (x)].
p
Then for every N > N0 , any T > Ω(N φw /εI ) and all i ∈ [3 : N − 3] there is ηi ∈ [−µ, µ] such
that
T
w(i/N ) X
εI
(1 + ηi )
Fi (j/[T + 1]) 6 ,
−
(1)
N
N
j=1

where Fi (x) , (w(x)/[T + 1]) · BN,i (x) and BN,i (x) ,

N i
i x (1

− x)N −i .

For every function w that satisfies the hypothesis in Theorem 2.1 we associate MDBD BN,T (α, p)
that is defined by pj = j/(T + 1) and αj = w(pj )/(T + 1) for all j ∈ [1 : T ]. Moreover, in Section 10
we give an efficient parallel algorithm that on input a continuous p.d.f. w (satisfying the conditions
in Theorem 2.2) and integer N ∈ N+ , outputs MDBD that induces a discretized p.d.f. which
approximates component-wise a desired discretized p.d.f. w.
b

Theorem 2.2 (An Efficient Parallel Algorithm for Finding MDBD). Let w(x) = C ·f (x) be a p.d.f.
that satisfies the hypothesis of Theorem 2.1. Suppose also a) 0 6 f (x) 6 1, b) 12 [f (0) + f (1)] >
R1
Ω(1), c) 1 6 C 6 o(N ), and d) 0 f (2) (x) dx 6 o(N ). Then there is a parallel algorithm
AppDscrPDF
that on input w(x) p
as above, integer N ∈ N+ and parameter εI > 0, outputs
p
in O(N φw /εI ) work and O(log(N φw /εI )) depth MDBD BN,T (α, p) that induces a discretized
probability distribution γ/[1 − δw ] over [0 : N ] such that for all i ∈ [3 : N − 3]


2εI
γi
∈ (1 + 2ηi ) · w(i/N
b
)±
,
1 − δw
SN +1,N
3

where the target discretized p.d.f. is defined by w(i/N
b
) , w(i/N )/SN +1,N for all i ∈ [0 : N ],
PN
P
ST,T +1 /(T +1)
δw , 1 − SN+1,N /N , SN +1,N , j=0 w(j/N ) and ST,T +1 , Tk=1 w(j/[T + 1]). Moreover, it
holds that δw ∈ [0, o(1)].
In Appendix D, we illustrate the representational power of MDBD by applying Theorem 2.2 for
two canonical continuous p.d.f.: the Uniform distribution and the Exponential Families.
Exact Recovery Interestingly, in case when a discretized p.d.f. w
b is induced by MDBD BN,T (α, p)
with exactly T = N + 1 distinct Binomial distributions, we give in Section 11 an efficient parallel
algorithm that on input vectors p and w,
b outputs the vector α in O(N log2 N ) work and O(logc N )
depth for some constant c ∈ N+ .

Theorem 2.3 (Canonical Instances Admit Exact Recovery). Suppose p ∈ (0, 1)N +1 is a vector
such that 0 < pi 6= pj < 1 for every i 6= j and w
b ∈ (0, 1)N +1 is a discretized p.d.f. that is induced
PN +1
by MDBD BN,N +1 (α, p) that satisfies w(i)
b
=
j=1 αj · BN,i (pj ) for every i ∈ [0 : N ]. Then
there is a parallel algorithm that on input the vectors p and w,
b outputs the vector α ∈ (0, 1)N +1 in
2
c
O(N log N ) work and O(log n) depth, for some constant c ∈ N+ .

2.2

Spectral Sparsification of Matrix Polynomials induced by MDBD

A matrix B is T -matrix if it is either Laplacian or SDDM matrix. To highlight that an algorithm A
preserves the matrix type, we write that the algorithm A on input a T -matrix B outputs a matrix
B ′ that is also T -matrix.
Moreover, we say that a matrix X is a spectral sparsifier of a matrix Y if it satisfies (1 − ε)Y 
X  (1 + ε)Y , for short X ≈ε Y , where the partial relation X  0 stands for X is symmetric
positive semi-definite (SPSD) matrix.
We denote by nnz(A) or mA the number of non-zero entries of matrix A. When we write
“B = D − M is T -matrix” we assume that D is positive diagonal matrix and B ∈ Rn×n . All
algorithms presented in this paper output spectral sparsifiers with high probability.
Sequential Algorithms Cheng et al. [5, Theorem 1.5] gave an algorithm that on input a
Laplacian matrix L = D − A, even integer N ∈ N+ and parameter ε > 0, outputs in time
b ≈ε D − D(D −1 A)N of a matrix-monomial
O(ε−2 mL log3 n · log2 N ) a spectral sparsifier D − A
−2
b 6 O(ε n log n). In Section 4, we give for N = 2k and k ∈ N+ a O(log2 N )such that nnz(A)
factor faster algorithm that computes a spectral sparsifier of T -matrix monomials. Furthermore,
for any T -matrix D − M such that M is SPSD matrix, we prove that the initial sparsification step
dominates the algorithm’s runtime.
Theorem 2.4 (Power Method for Monomials). There is an algorithm PwrSS that on input T cN ≈ε
matrix B = D − M , N = 2k for k ∈ N+ and ε ∈ (0, 1), outputs a spectral sparsifier D − M
−1
N
−2
cN ) 6 O(ε n log n). The algorithm runs in time
D − D(D M ) that is T -matrix with nnz(M
(
e B log2 n + ε−4 · n log4 n · log5 N )
O(m
, if M is SPSD matrix;
3
4
5
−2
−4
e
O(ε mB log n + ε · n log n · log N ) , otherwise.

Using Theorem 2.4, we give in Section 5 an algorithm that runs by Θ(N 2 )-factor faster than [4,
Theorem 2] P
and computes a spectral sparsifier of a single Binomial T -matrix polynomials of the
−1
i
N
−1
form D − D N
i=0 BN,i (p) · (D M ) = D − DWp , where Wp = (1 − p)I + pD M and p ∈ (0, 1).
4

Theorem 2.5 (Single Binomial Matrix Polynomials). There is an algorithm LazySS that on input
T -matrix B = D − M , number N = 2k for k ∈ N+ , and parameters ε, p ∈ (0, 1), outputs a spectral
sparsifier
N
X
N
c
BN,i (p) · (D −1 M )i
D − Mp,N ≈ε D − DWp = D − D
i=0

O(ε−2 n log n)

that is T -matrix with at most
non-zero entries. The algorithm LazySS runs in time
(
e B log2 n + ε−4 · n log4 n · log5 N )
O(m
, if p ∈ (0, 1/2];
3
2
4
5
−2
−4
e
O(ε mB log n · log N + ε · n log n · log N ) , otherwise.

In Section 6, we give our main sequential algorithm that builds upon Theorem 2.5 and computes
a spectral sparsifier of T -matrix polynomials induced by MDBD.
Theorem 2.6 (Mixture of Binomial Matrix Polynomials). There is an algorithm SS MDBD
P that
on input T -matrix B = D − M , integer N = 2k , k ∈ N+ , MDBD BN,T (α, p) with δ = 1 − Ti=1 αi
c ≈ε D − D PN (γi /[1 − δ]) · (D −1 M )i
and parameter ε ∈ (0, 1), outputs a spectral sparsifier D − M
i=0
−2 n log n) non-zero entries, where γ is a discretized p.d.f. that satisfies
that is
T
-matrix
with
O(ε
P
γi = Tj=1 αj · BN,i (pj ) for all i. The algorithm runs in time
(
e B log2 n + ε−4 · nT · log4 n · log5 N )
O(m
e −2 mB log3 n + ε−4 · nT · log4 n · log5 N )
O(ε

, if M is SPSD matrix;
, otherwise.

We motivate now the first conclusion of Theorem 2.6. When B = D − DWp is a Laplacian
matrix (of a graph with self-loops) associated with a Markov Chain with transition matrix Wp that
corresponds to p-lazy random walk process, it holds for p ∈ (0, 1/2] (c.f. Lemma 5.1) that DWp is
SPSD matrix.
Given MDBD BN,T (α, p) that induces a vector γ, the algorithm in [4, Theorem 2] outputs a
b −2 mN 2 + N T )6 , where
spectral sparsifier of the corresponding T -matrix polynomial in time O(ε
the term O(N T ) accounts for computing the vector γ. In comparison, our improved algorithm
b −2 m + ε−4 nT ) for any MDBD with N = 2k for k ∈ N+ .
SS MDBD runs in time O(ε

Parallel Algorithms Building upon the seminal works of Spielman and Teng [25], Orecchia and
Vishnoi [20] and Spielman and Peng [21], we prove in Section 7 that algorithm SS MDBD can be
efficiently parallelized. To the best of our knowledge, this is the first efficient and parallel algorithm
that sparsifies T -matrix polynomials induced by MDBD with N = 2k for k ∈ N+ .
Theorem 2.7 (Efficient Parallel Spectral Sparsification of Matrix Polynomial induced by MDBD).
on input as in Theorem 2.6, outputs a spectral sparThere is a parallel algorithm
PN pSS MDBD that
−1
c
c) 6 O(ε−2 n logc n)
sifier D − M ≈ε D − D i=0 (γi /[1 − δ]) · (D M )i that is T -matrix with nnz(M
PT
for some constant c, where γ is a discretized p.d.f. such that γi = j=1 αj · BN,i (pj ) for all i.
e −2 mB logc1 +1 n · log2 N + ε−4 · nT · logc+c1+1 n · log5 N ) and depth
The algorithm runs in work O(ε
O(logc2 n · log N + log T ) for constants c1 , c2 ∈ N.
By combining Theorem 2.2 and Theorem 2.7, we develop an efficient parallel algorithm that
outputs a spectral sparsifier of a T -matrix polynomial whose coefficients approximate componentwise a target discretized p.d.f. w.
b
6

b notation hides poly(log n, log N ) factors.
O(·)

5

Corollary 2.8 (Approximating Target Transition Matrices). There is a parallel algorithm that
takes as input a continuous p.d.f. w satisfying the conditions of Theorem 2.2,
p T -matrix B =
b −2 mB + ε−4 nN φw /εI ) work and
D − M and parameters ε, εI ∈ (0, 1), and it outputs in O(ε
p
c ≈ε D − D PN (γi /[1 − δw ]) ·
O(logc2 n · log N + log(N φw /εI )) depth a spectral sparsifier D − M
i=0
c) 6 O(ε−2 n logc n) such that for all i ∈ [3 : N − 3] it holds
(D −1 M )i that is T -matrix with nnz(M
γi /[1 − δw ] ∈ [(1 + 2ηi ) · w(i/N
b
) ± 2εI /SN +1,N ].

2.3

Analyzing the Long Term Behaviour of Markov Chains

For many finite Markov chains [1, 10, 14, 23] there exists N ′ ∈ N+ such that for every N > N ′
certain phenomenon occurs with high probability - (local) mixing time, truncated PageRank, etc.
Therefore, to analyze the long term behaviour of a finite Markov chain, it suffices to select the
smallest N = 2k for k ∈ N+ that is larger or equal to N ′ .
Corollary 2.9 (Capturing The Long Term Behaviour of Markov Chains). Suppose LP= D − A is
dense Laplacian matrix with mL = Θ(n2 ), ε ∈ (0, 1), BN,T (α, p) is MDBD such that Ti=1 αi = 1,
the degree N = 2k 6 O(n) for k ∈ N+ and the number of Binomials T 6 O(n). Then algorithm
e −4 · mL · logc+c1 +6 n) and depth O(logc2 +1 n) a spectral sparsifier
O(ε
pSS MDBD outputs
PN in work
−1
b
D − A ≈ε D − D i=0 γi (D A)i that is Laplacian matrix with O(ε−2 n logc2 n) non-zero entries for
some constants c, c1 , c2 ∈ N and γ is a probability distribution induced by the MDBD BN,T (α, p).
Multiplicative Approximation P
of Generalized Escaping Probability Consider a Markov
−1
i
chain with transition matrix Gγ = N
i=0 γi (D A) that corresponds to a generalized random walk
process of length N . Perform a random walk of length N induced by Gγ that starts at vertex
v ∈ V . Then for any subset S ⊂ V , the corresponding generalized escaping probability is defined
by gEsc(v, S, Gγ ) = 1T
v Gγ 1S , where we denote
P by 1T the characteristic vector of a subset T ⊂ V .
We define the volume of S by µ(S) = u∈S du and let πS be a probability distribution over V
defined by πS (u) = du /µ(S) if u ∈ S and πS (u) = 0 otherwise. The expected generalized escaping
probability (E.G.E.P.) with respect to πS is defined by
Ev∼πS [gEsc(v, S, Gγ )] = πST Gγ 1S .

(2)

We show in Appendix A that a spectral sparsifier of a random walk Laplacian matrix polynomial,
yields a multiplicative approximation of E.G.E.P. for all subsets S ⊂ V .

bN ≈ε
Lemma
(Multiplicative Approximation of E.G.E.P.). For any spectral sparsifier D − A
P2.10
N
−1
i
D − D i=0 γi (D A) of a random walk Laplacian matrix polynomial such that γ is a probability
distribution over [0 : N ], it holds for every subset S ⊂ V that
p
bN )ξS ∈ [ (1 ± ε) · Ev∼π [gEsc(v, S, Gγ )] ], where ξS , 1S / µ(S).
ξST (D − A
S

Using Corollary 2.9 and Lemma 2.10, we propose the first efficient and parallel algorithm that
runs in nearly linear work and poly-logarithmic depth that yields a multiplicative approximation
of E.G.E.P. for Markov chains with transition matrices induced by MDBD.

2.4

Faster SDDM Solver

Spielman and Peng [21] gave the first parallel SDD solver that constructs in O(m logc1 n·log3 κ) work
and O(logc2 n · log κ) depth a sparse O(1)-approximate inverse chain that solves approximately to
any ε > 0 precision an SDD system in O((m+n logc n·log3 κ) log 1/ε) work and O(log n·log κ·log 1/ε)
depth. In Section 8, we give a simple parallel preprocessing step that strengthens their algorithm.
6

Theorem 2.11. There is an algorithm that on input an n-dimensional SDDM matrix M with m
non-zeros and condition number at most κ, produces with probability at least 1/2 a sparse O(1)approximate inverse chain that can be used to solve any linear equation in M to any precision ε > 0
in O(n logc n · log3 κ · log 1/ε) work and O(log n · log κ · log 1/ε) depth, for some constant c. The
algorithm runs in O(m logc1 n) work and O(logc2 n · log κ) depth for some other constants c1 , c2 .
For the current state-of-the-art result on parallel SDD solvers we refer the reader to the work
of Lee et al. [17].

3

Algorithmic Background on Spectral Sparsification

We write X ≈ε1 ⊕ε2 Y to indicate (1 − ε1 )(1 − ε2 )Y  X  (1 + ε1 )(1 + ε2 )Y . Our analysis uses
the following five basic facts (c.f. [2, 26]).
Fact 3.1. For positive semi-definite (PSD) matrices X, Y, W and Z it holds
a. if Y ≈ε Z then X + Y ≈ε X + Z;
b. if X ≈ε Y and W ≈ε Z then X + W ≈ε Y + Z;
c. if X ≈ε1 Y and Y ≈ε2 Z then X ≈ε1 ⊕ε2 Z;
d. if X and Y are invertible matrices such that X ≈ε Y then X −1 ≈2ε Y −1 , ∀ε ∈ (0, 21 );
e. for any matrix V if X ≈ε Y then V T XV ≈ε V T Y V .

3.1

Prior Algorithms

Our algorithms for computing spectral sparsifiers of matrix-polynomials use as a black-box several
spectral sparsification algorithms for Laplacian and SDDM matrices.
More precisely, our sequential algorithms build upon Theorem 3.2 that relies on Kelner and
Levin’s [13, Theorem 3] and Cohen et al.’s [6, Lemma 4], and Theorem 3.3 proposed by Spielman
and Peng [21, Corollary 6.4].
Theorem 3.2. [13] There is an algorithm SS takes as input parameter ε ∈ (0, 1), matrices D and
A such that D is positive diagonal and A is symmetric non-negative with Aii = 0 for all i such that
e L log2 n) time outputs a positive diagonal matrix D
e
L = D − A is Laplacian matrix. Then in O(m
−2
e such that nnz(A)
e 6 O(ε n log n), A
eii = 0 for all i, and
and symmetric non-negative matrix A
e
e
e
e
D − A ≈ε D − A. Moreover, D − A is Laplacian matrix.

Theorem 3.3. [21] There is an algorithm PS that takes as input SDDM matrix B = D − M
e and
and parameter ε ∈ (0, 1). Then in O(ε−2 mB log2 n) time outputs a positive diagonal matrix D
−2
f with nnz(M)
f 6 O(ε mB log n) and M
fii = 0 for all i, such that
symmetric non-negative matrix M
−1
e
f
e
e
f
D − M ≈ε D − M D M and D ≈ε D. Moreover, D − M is SDDM matrix.

Our parallel algorithm uses Theorem 3.4, which is the culmination of a research line conducted
by Spielman and Teng [25], Orecchia and Vishnoi [20] and Spielman and Peng [21].

Theorem 3.4. [21] There is an algorithm STOVP takes as input a Laplacian matrix D − M and
f ≈ε D − M with M
fii = 0 for all i
parameter ε ∈ (0, 1/2). Then it outputs a spectral sparsifier D − M
c
−2
f
and nnz(M ) 6 O(ε n log n) for some constant c. Moreover, this algorithm requires O(m logc1 n)
work and O(logc2 n) depth, for some other constants c1 and c2 .
Based on Theorem 3.4 Spielman and Peng [21] parallelized algorithm PS (c.f. Theorem 3.3).

7

Theorem 3.5. [21] There is a parallel algorithm that on input an SDDM matrix D − M and
f ≈ε D − M D −1 M with D
e ≈ε D, M
fii = 0
parameter ε ∈ (0, 1/2), outputs a spectral sparsifier D − M
c
f) 6 O(ε−2 n log n) for some constant c. Moreover, this algorithm requires
for all i and nnz(M
c1 +1
−2
O(ε m log
n) work and O(logc2 n) depth, for some other constants c1 and c2 .

3.2

Spectral Sparsification of T -Matrices

We show that the algorithms SS and PS can be amended to produce T -matrix sparsifiers that are
in normalized form, i.e. the sparsifiers are expressed in terms of the diagonal matrix D minus a
c. Our analysis relies on several results established by Peng et
symmetric non-negative matrix M
al. [3, 4, 21, 22].

Lemma 3.6. There is an algorithm mSS that takes as input a positive diagonal matrix D, symmetric non-negative matrix A (possibly Aii 6= 0) such that B = D − A is Laplacian matrix and
e B log2 n) time a spectral sparsifier D− A
b ≈ε D−A that
parameter ε ∈ (0, 1). Then it outputs in O(m
b is symmetric non-negative matrix with nnz(A)
b 6 O(ε−2 n log n).
is Laplacian matrix and satisfies A

The next result implicitly appears in [4]. For completeness we prove it in Appendix B.
e −A
e a sparsifier with
Lemma 3.7. Suppose D − A is Laplacian matrix (possibly Aii 6= 0) and D
eii = 0 for every i such that (1 − ε)(D − A)  D
e−A
e  (1 + ε)(D − A). Then the symmetric
A
1 e
1 e
b  (1 + 2ε)(D − A).
b
non-negative matrix A = (D − 1+ε D) + 1+ε A satisfies (1 − 2ε)(D − A)  D − A
We present now the proof of Lemma 3.6.

Proof of Lemma 3.6. Notice that D − A = D ′ − A′ , where D ′ is positive diagonal matrix and
A′ is symmetric non-negative matrix such that A′ii = 0 for all i. By Theorem 3.2 we obtain
b′ ≈ε D ′ − A′ , where
e′ − A
e′ ≈ε/2 D ′ − A′ . Then by Lemma 3.7 we have D ′ − A
a sparsifier D
b′ = (D ′ − 1 D
e ′) + 1 A
e′ is symmetric non-negative matrix. We define by DA = D − D ′ a
A
1+ε
1+ε
b = DA + A
b′ and observe that it is symmetric and non-negative
non-negative diagonal matrix. Set A
b = D′ − A
b′ .
matrix. Now the statement follows since D − A

T -Matrices Building upon the work of Spielman and Peng [21, Proposition 5.6] and Cheng et
al. [4, Proposition 25], we prove in Appendix B the following statement.

Lemma 3.8 (Closure). Suppose D − M is T -matrix. Then D − D(D −1 M )N is T -matrix for every
c ≈ε D − M is a spectral sparsifier, then D − D(D −1 M
c)N is T -matrix
N ∈ N+ . Moreover, if D − M
for every N ∈ N+ .
Normalized Algorithms We present now two algorithms that sparsify matrices of the form
D − D(D −1 M )N for N ∈ {1, 2} such that the resulting sparsifiers are in normalized form.
Lemma 3.9 (Normalized Spectral Sparsification). There is an algorithm mKLC that takes as
e B log2 n) time a
input T -matrix B = D − M and parameter ε ∈ (0, 1), then it outputs in O(m
c ≈ε D − M that is T -matrix and M
c is symmetric non-negative matrix
spectral sparsifier D − M
−2
c) 6 O(ε n log n).
with nnz(M

Proof. By definition B = D1 + L where D1 is non-negative diagonal matrix and L = D2 − M is
c ≈ε D2 − M that is Laplacian
Laplacian matrix. We obtain by Lemma 3.6 a sparsifier D2 − M
matrix. Now we consider two cases. If D1 = 0 then we are done. Otherwise D1 is PSD matrix
c ≈ε D − M . Since D − M is SDDM matrix and the operator ≈ε
and by Fact 3.1.a we have D − M
c is SDDM matrix.
preserves the kernel space, it follows that D − M

8

We proceed by stating an interesting structural result that implicitly appears in [21] (c.f. Section
“Efficient Parallel Construction”). For completeness we prove it in Appendix B.1.
T − M · 1 be a column vector,
Lemma 3.10. Suppose B = D − M is T -matrix. Let ηi = Mi,:
i,i
i
di = hMi,: , 1i and si = di − Mii numbers, and DNi = (si /di ) · diag(Ni ) positive diagonal matrix
for all i, where Ni = {Mij | Mij 6= 0}. Let Bij = (Mii /di + Mjj /dj ) · Mij be the (i, j)th entry of a
matrix with same dimensions as matrix M and DBP= diag(B · 1) be a diagonal matrix.
Then it holds that D − M D −1 M = D1 + LB + ni=1 LNi where D1 = diag([D − M D −1 M ]1) is
non-negative diagonal matrix, LB = DB − B is Laplacian matrix with at most mB non-zero entries
and every LNi = (si /di )DNi − ηi ηiT /di is Laplacian matrix corresponding to a clique with positively
weighted edges that is induced by the neighbour set Ni .

Spielman and Peng [21] gave algorithm PS (c.f. Theorem 3.3) for sparsifying matrices of the
form D − M D −1 M , where D − M is SDDM matrix. We extend their result to T -matrices and our
algorithm outputs a spectral sparsifier in normalized form.
Lemma 3.11 (Normalized 2-Hops Spectral Sparsification). There is an algorithm mPS that on
e −2 mB log3 n) time a spectral
input a T -matrix B = D − M and parameter ε ∈ (0, 1), outputs in O(ε
c ≈ε D − M D −1 M that is T -matrix and M
c is symmetric non-negative matrix with
sparsifier D − M
−2
c
nnz(M ) 6 O(ε n log n).

Proof. By Lemma 3.6 we have D − M D −1 M = D1 + L, where D1 is non-negative diagonal matrix
and L is sum of Laplacian matrices. Using similar arguments as in “Section 6 Efficient Parallel
e −M
f ≈ε/2 L. Moreover, we can compute the positive
Construction” [21] we find a sparsifier D
′
diagonal matrix D = diag(L) in O(mB ) time (c.f. Appendix B.1), and then by Lemma 3.7 we
c ≈ε L. Since D1 is PSD matrix the statement follows by Fact 3.1.a. 
obtain a sparsifier D ′ − M

4

Core Iterative Algorithm

Our goal now is to prove Theorem 2.4. We argue in a similar manner as in [5], but in contrast
our analysis shows that the initial sparsification step tolerates higher approximation error. This
observation yields an improved algorithm whose runtime is faster by a O(log2 N )-factor.
Moreover, we prove that for any T -matrix D − M such that M is SPSD matrix, one can
c2 ≈ D − M D −1 M by first computing D − M
c ≈ D − M and
construct a spectral sparsifier D − M
−1
c
c
c
then D − M2 ≈ D − M D M . This demonstrates that when M is SPSD matrix, the runtime is
dominated by the initial sparsification.
The rest of this section is organized as follows. In Subsection 4.1 we describe the initial phase
of algorithm PwrSS. Then in Subsection 4.2, we present the iterative construction of the desired
cN ≈ε D − D(D −1 M )N .
spectral sparsifier D − M

4.1

Initialization

We begin by extending [5, Lemma 4.3 and 4.4]. For completeness, we provide a prove in Appendix E
where in addition we generalize [5, Fact 4.2].
c ≈ε D − M is a spectral sparsifier. If
Lemma 4.1. Suppose B = D − M is T -matrix and D − M
−1
c
c
M is SPSD matrix then it holds that D − M D M ≈ε D − M D −1 M .

Based on Lemma 4.1, we give a faster sparsification algorithm for T -matrices D − M D −1 M
such that M is SPSD matrix.
9

Lemma 4.2. There is an algorithm fSS that takes as input T -matrix B = D − M such that M
e B log2 n + ε−4 n log4 n) time a
is SPSD matrix, and parameter ε ∈ (0, 1). Then it outputs in O(m
−1
c2 ≈ε D − M D M that is T -matrix with nnz(M
c2 ) 6 O(ε−2 n log n).
spectral sparsifier D − M
c ≈ε/4 D−M in O(m
e B log2 n)
Proof of Lemma 4.2. We apply Lemma 3.9 to obtain a sparsifier D−M
c 6 O(ε−2 n log n) such that D − M
c is T -matrix. Then by Lemma 4.1 we know
time with nnz(M)
cD −1 M
c ≈ε/4 D − M D −1 M . Now, by Lemma 3.8 D − M
cD −1 M
c is T -matrix. Then we
that D − M
e −4 n log4 n) time a sparsifier D − M
c2 ≈ε/4 D − M
cD −1 M
c with
apply Lemma 3.11 to obtain in O(ε
c) 6 O(ε−2 n log n) such that D − M
c2 is T -matrix. The claims follows by Fact 3.1.c.
nnz(M


4.2

Iterative Construction

Our analysis of the incurred approximation error after O(log N ) consecutive square sparsification
operations builds upon [5, Lemma 4.1]. In contrast, we prove that for the initial and the final
sparsifiers it suffices to have only an ε approximation, while all intermediate spectral sparsifiers
require finer ε′ = Ω(ε/ log N ) approximation. Due to this higher initial error tolerance, we improve
the runtime of their algorithm by a O(log2 N )-factor.
c2 be T -matrices such that D − M
c2 ≈ε
Lemma 4.3 (Accumulative Error). Let D − M and D − M
−1
−2
c
D − M D M and nnz(M2 ) 6 O(ε n log n). There is an algorithm IndSS that on input T c2 , integer N = 2k for k ∈ N+ and parameter 0 < ε′ 6 ε, outputs in time
matrix D − M
−4
4
′
e
cN with nnz(M
cN ) 6 O(ε−2 n log n) such
O(ε n log n · log N ) a symmetric non-negative matrix M
cN ≈ (log N−1) ′ 2 D − D(D −1 M )N is T -matrix.
that D − M
(⊕
ε )⊕ ε
Our goal now is to prove Lemma 4.3. We establish next a useful algebraic property that all
k
matrices of the form D(D −1 M )2 have in common.
k

Lemma 4.4. If M is symmetric matrix, then D(D −1 M )2 is SPSD matrix for every k ∈ N+ .
k

k

Proof. Let Y , D −1/2 M D −1/2 . Notice that D(D −1 M )2 = D 1/2 Y 2 D 1/2 = X T X, where X =
k−1

Y 2 D 1/2 . The statement follows since X T X is SPSD matrix.
We present now the main iterative procedure used in algorithm IndSS.
c2k ≈ε D−
c2k be T -matrices such that D−M
Lemma 4.5 (Iterative Procedure). Let D−M and D−M
k
−1
2
c2k
D(D M ) , for k ∈ N+ . There is an algorithm SqrSS that takes as input the T -matrix D − M
−2
e ′ nnz(M
c k ) log3 n) time a symmetric non-negative
and parameter ε′ ∈ (0, 1), then it outputs in O(ε
2
−2
′
c2k+1 ≈ε⊕ε′ D − D(D −1 M )2k+1 is
c2k+1 ) 6 O(ε n log n) such that D − M
c2k+1 with nnz(M
matrix M
T -matrix.
k

Proof. By Lemma 4.4, D(D −1 M )2 is SPSD matrix for any k ∈ N+ . By Lemma 3.8 both D −
k
c2k are T -matrices. Hence, by Lemma 4.1 we have that D −
c2k D −1 M
D(D −1 M )2 and D − M
k+1
c2k ≈ε D − D(D −1 M )2 . Now by Lemma 3.11 we have D − M
c2k D −1 M
c2k
c2k D −1 M
c2k+1 ≈ε′ D − M
M

and hence the statement follows by Fact 3.1.c.
Based on the preceding results we are ready to prove Lemma 4.3.

e ′ ·ε)−2 n log4 n) we can compute a spectral sparsiProof of Lemma 4.3. By Theorem 3.3 in time O((ε
−1
4
c4 ) 6 O(ε′ −2 n log n). Then we apply (log N −1) times
c4 ≈ε⊕ε′ D−D(D M ) with nnz(M
fier D− M
e ′ −4 n log4 n · log N ) time a spectral sparsifier D − M
cN ≈ (log N−1) ′
Lemma 4.5 to obtain in O(ε
(⊕
ε )⊕ε
cN ) 6 O(ε′ −2 n log n). The statement follows by applying Theorem 3.2
D − D(D −1 M )N with nnz(M
cN .
with ε to compute a refined spectral sparsifier of D − M

10

c2 ≈ε/4 D − M D −1 M
Proof of Theorem 2.4 In the initial phase we compute a sparsifier D − M
c2 ) 6 O(ε−2 n log n) using either Lemma 4.2 (when M is SPSD matrix) or Lemma 3.11.
with nnz(M
c2 .
The statement follows by applying Lemma 4.3 with ε′ = Ω(ε/ log N ) to the sparsifier D − M

5

Spectral Sparsification of Binomial T -Matrix Polynomials

In this section, we prove Theorem 2.5. We analyze first the properties of matrices of the form
Wp = (1 − p)I + pD −1 M for p ∈ (0, 1). It is convenient to associate with them matrix-polynomials
fp (x) = (1 − p) + px such that fp (D −1 M ) = Wp . Since the coefficients of the matrix-polynomial
P
−1
i
[fp (x)]N follow Binomial distribution B(N, p), it follows that WpN = N
i=0 BN,i (p) · (D M ) for

every p ∈ (0, 1), where BN,i (p) = Ni pi (1 − p)N −i .
When D − M is a Laplacian matrix, the matrix WpN corresponds to the transition matrix of a
associate to such a Markov chain a matrixp-lazy random walk process of length
P N (c.f. [23]). We
−1 M )i . We present now some useful algebraic
polynomial D − DWpN = D − D N
B
(p)
·
(D
i=0 N,i
k
properties of matrices of the form DWp2 and D − DWpN .
Lemma 5.1. Suppose D − M is T -matrix. Then D − DWpN is T -matrix for every N ∈ N+ . Also
k
DWp is SPSD matrix ∀p ∈ (0, 1/2] and DWp2 is SPSD matrix ∀p ∈ (0, 1) and ∀k ∈ N+ .
Proof. By definition of Wp , we have D − DWp = p(D − M ) is T -matrix. Suppose D − M is
Laplacian matrix, then D − DWp is Laplacian matrix and by Lemma 3.8, D − DWpN is Laplacian
matrix for every N ∈ N+ . Suppose now that D − M is SDDM matrix, then D − DWp is SDDM
matrix and by Lemma B.2 D − DWpN is SDDM matrix for every N ∈ N+ .
By definition DWp = (1 − p)D + pM and since D − M is diagonally dominant, it holds that
DWp is SPSD matrix for every p ∈ (0, 1/2]. Moreover, since DWp is symmetric matrix by Lemma
k
k

4.4 it holds that D(D −1 · DWp )2 = DWp2 is SPSD matrix for every k ∈ N+ .
Proof of Theorem 2.5. The statement follows by Lemma 5.1 and Theorem 2.4.

6



Spectral Sparsification of T -Matrix Polynomials Induced by MDBD

Here we prove Theorem 2.6. Our approach relies on the following key algorithmic idea.
c1 ≈ε
Lemma 6.1 (Preprocessing of 2-Hop Spectral Sparsification). Let D−M be a T -matrix, D− M
−1
c
D − M and D − M2 ≈ε D − M D M are spectral sparsifiers. Then for every p ∈ (0, 1) the sparse
cp,2 = (1 − p)2 D + 2(1 − p)pM
c1 + p2 M
c2 yields a spectral sparsifier D − M
cp,2 ≈ε D − DWp2
matrix M
that is T -matrix.
Proof. The statement follows by
cp,2
D−M

=

c1 ] + p2 [D − M
c2 ]
2p(1 − p)[D − M

≈ε 2p(1 − p)[D − M ] + p2 [D − M D −1 M ] = D − DWp2 .


Our algorithm SS MDBD builds upon Lemma 6.1 and Lemma 4.3. Due to the preprocessing
step in Lemma 6.1 we speed up the sparsification of each T -matrix polynomial D − DWpNj for all
j ∈ [1 : T ]. We present now the pseudo code of algorithm SS MDBD.
11

Algorithm 1
c) = SS MDBD(D, M, BN,T (α, p), ε)
(D, M
ctmp = 0 and δ = 1 − PT αi .
1. Let ε′ = ε/[3 log N ], M
i=1
c1 , M
c2 ) = InitSS(D, M, ε/3).
2. (D, M
3. For every j ∈ {1, . . . , T } do
cp,2 = (1 − p)2 D + 2p(1 − p)M
c1 + p2 M
c2 .
3.1 Set p = pj,T = j/(T + 1) and M
cp,N ) = IndSS(M
cp,2 , N, ε′ ) where D − M
cp,N ≈2ε/3 D − DW N (c.f. Lemma 4.3).
3.2 (D, M
p
ctmp = M
ctmp + αj · M
cp,N .
3.3 M
ctmp by algorithm mKLC (c.f. Lemma 3.9).
c ≈ε/3 D − 1 M
4. Sparsify D − M
1−δ
c).
5. Return (D, M
Algorithm 2
c1 , M
c2 ) = InitSS(D, M, ε)
(D, M
c1 ≈ε D − M by algorithm mKLC (c.f. Lemma 3.9).
1. Sparsify D − M
c2 ≈ε D − M D −1 M
2. Sparsify D − M
2.1 If M is SPSD matrix call algorithm fSS (c.f. Lemma 4.2),
2.2 otherwise call algorithm mPS (c.f. Lemma 3.11).
c1 , M
c2 ).
3. Return (D, M

P
Let δ = 1 − Ti=1 αi . We denote a T -matrix polynomial induced by MDBD BN,T (α, p) as
P
P
PT
−1
i
PB , Tj=1 αj (D − DWpNj ) = (1 − δ)D − D N
i=0 γi (D M ) , where γi =
j=1 αj BN,i (pj ). We are
now ready to prove Theorem 2.6.
Proof of Theorem 2.6. Let ε′ = ε/[4 log N ]. We perform first a preprocessing step. We apply
c ≈ε′ D − M . Then depending on whether M is SPSD
Lemma 3.9 to obtain a sparsifier D − M
c2 ≈ε′ D−M D −1 M . The
matrix we use either Lemma 3.11 or Lemma 4.2 to obtain a sparsifier D− M
3
2
2
4
−2
−4
e
e
run time is at most O(ε mB log n · log N ) or O(mB log n + ε n log n) respectively. Moreover,
c1 ), nnz(M
c2 ) 6 O(ε−2 n log n · log2 N ).
the sparsifiers satisfy nnz(M
cp ,N ≈ε D − DW N by
We combine Lemma 6.1 and Theorem 2.4 to find each sparsifier D − M
pj
j
PT
2
c
c
c
αj Mp ,N ,
initializing algorithm PwrSS with a sparsifier D−Mp ,2 ≈ε′ D−DWp . Let Mtmp =
j

j

j=1

j

ctmp ≈ε PB . This phase has O(ε
e −4 nT log4 n ·log5 N ) runtime
then by Fact 3.1.b it holds (1− δ)D − M
cp ,N ) 6 O(ε−2 n log n).
and each sparsifier satisfies nnz(M
i
c
c ≈ε
However, matrix Mtmp can be dense. We apply Lemma 3.9 to obtain a sparsifier D − M
2
3
2
1 c
1
2
−2
e
c) 6
D − 1−δ Mtmp ≈ε 1−δ PB in time O(min{n
log n, ε nT log n · log N }) such that nnz(M
−2
O(ε n log n).


7

Parallelization of Algorithm SS MDBD

In this section we parallelize algorithm SS MDBD. This gives the first efficient parallel algorithm
that computes a spectral sparsfier of any T -matrix polynomial with coefficients induced by MDBD.
Our goal now is to prove Theorem 2.7. By construction of algorithms SS MDBD and InitSS, it
suffices to show that we can efficiently parallelize algorithms mKLC, mPS and PwrSS. Then the
statement follows by noting that each T Binomial matrix-polynomial can be sparsified separately
and in parallel.
12

We parallelize now algorithms mKLC, mPS and IndSS.
Lemma 7.1. There is a parallel algorithm pKLC that on input T -matrix B = D − M and
c ≈ε D − M that is T -matrix such that M
c is
parameter ε ∈ (0, 1), outputs a spectral sparsifier D − M
c
c) 6 O(ε−2 n log n) for some constant c. The algorithm
symmetric non-negative matrix with nnz(M
c1
c2
runs in O(mB log n) work and O(log n) depth, for some other constants c1 , c2 .

Proof. We argue in a similar manner as in Lemma 3.9 to show that the statement holds for T 
matrices (Laplacian or SDDM matrices). Then the statement follows by Theorem 3.4.
Lemma 7.2. There is a parallel algorithm pPS that on input T -matrix B = D − M and parameter
f ≈ε D − M D −1 M that is T -matrix such that
ε ∈ (0, 1/2), outputs a spectral sparsifier D − M
−2
e
f
f
D ≈ε D, Mii = 0 for all i and nnz(M ) 6 O(ε n logc n) for some constant c. The algorithm runs
in O(ε−2 m logc1 +1 n) work and O(logc2 n) depth, for some other constants c1 and c2 .
Proof. Using similar arguments as in Lemma 3.11 we prove that the statement holds for T -matrices.

Then the statement follows by Theorem 3.5.

c2k ≈ε D − D(D −1 M )2k ,
c2k are T -matrices such that D − M
Lemma 7.3. Let D − M and D − M
c2k and
for k ∈ N+ . There is a parallel algorithm pSqrSS that on input the T -matrix D − M
k+1
−1
2
′
c
that is
parameter ε ∈ (0, 1), outputs a spectral sparsifier D − M2k+1 ≈ε⊕ε′ D − D(D M )
−2
c
′
c
T -matrix with nnz(M2k+1 ) 6 O(ε n log n) for some constant c. The algorithm runs in work
c2k ) logc1 +1 n) and depth O(logc2 n), for some other constants c1 , c2 .
O(ε′ −2 nnz(M

Proof. We use similar arguments as in Lemma 4.5, but we substitute Lemma 3.11 with Lemma 7.2.


c2 are T -matrices such that D − M
c2 ≈ε D − M D −1 M and
Lemma 7.4. Let D − M and D − M
c
−2
c
nnz(M2 ) 6 O(ε n log n) for some constant c. There is a parallel algorithm pIndSS that on input
c2 , integer N = 2k for k ∈ N and parameter 0 < ε′ 6 ε, outputs a spectral sparsifier
T -matrix D − M
c
cN ) 6 O(ε−2 n logc n). The
D − MN ≈(⊕(log N−1) ε′ )⊕2 ε D − D(D −1 M )N that is T -matrix and nnz(M
e ′ −4 n logc+c1 +1 n) and depth O(logc2 n · log N ) for some constants c1 , c2 .
algorithm runs in work O(ε

Proof. We argue in a similar manner as in Lemma 4.3. By Lemma 7.1 we compute a sparsifier
c4 ≈ε′ ⊕ε D − D(D −1 M )4 with nnz(M
c4 ) 6 O(ε′ −2 n logc n) in work O((ε · ε′ )−2 n logc+c1 +1 n)
D−M
c2
and depth O(log n). Then we apply (log N − 1) times Lemma 7.3 to obtain a spectral sparsifier
−1
N
′ −2 n logc n) in O(ε′ −4 n logc+c1 +1 n)
cN ≈ (log N−1) ′
c
D−M
(⊕
ε )⊕ε D − D(D M ) with nnz(MN ) 6 O(ε
c2
work and O(log · log N ) depth. The statement follows by applying Lemma 7.1 with ε to compute
cN .
a refined spectral sparsifier of D − M

We present now the proof of Theorem 2.7 which yields the parallel algorithm pSS MDBD.

Proof of Theorem 2.7. We sketch first our parallel algorithm pSS MDBD. We parallelize algorithm InitSS based on algorithms pKLC and pPS. Then, we sparsify separately and in parallel
each of the T distinct single Binomial T -matrix polynomials by algorithm pIndSS. The resulting
T sparsifiers are scaled and merged into a T -matrix polynomial induced by MDBD. Since this
matrix-polynomial might be dense, we sparsify it using algorithm pKLC.
The correctness of algorithm pSS MDBD follows by Theorem 2.6. We analyze now the work
and the depth of algorithm pSS MDBD. By Lemma 7.1 and Lemma 7.2 the initial phase is dominated by O(ε−2 m logc1 +1 n · log2 N ) work and O(logc2 n) depth. Moreover, each of the sparsifiers
c ≈ε D − M and D − M
c2 ≈ε D − M D −1 M has at most O(ε−2 n logc n) non-zero entries.
D−M
13

Let ε′ = ε/[12 log N ]. By Lemma 7.4 for each j ∈ [1 : T ] we apply algorithm pIndSS to
cp ,N ≈ε/6 D − DWpN with nnz(M
cp ,N ) 6 O(ε−2 n logc n). This
compute a spectral sparsifier D − M
j
j
j
e ′ −4 nT logc+c1 +1 n) and depth O(logc2 n · log N ). Furthermore, the linear
phase runs in work O(ε
ctmp = PT αj · M
cp,N (in Step 3.3) can be computed in depth O(log T ).
combination M
j=1
P
1
1 c
−1
i
Mtmp ≈ε/6 D − D N
Therefore, we approximate by D − 1−δ
i=0 1−δ γi (D M ) the desired T ctmp might be dense, by Lemma 7.1 we
matrix polynomial induced by MDBD. Since matrix M
ctmp with nnz(M
c) 6 O(ε−2 n logc2 n). Algorithm
c ≈ε/12 D− 1 M
compute a spectral sparsifier D−M
1−δ

pKLC runs in work O(min{ε−2 nT logc n, n2 } · logc1 n) and depth O(logc2 n).

8



Faster SDDM Solver

In this section we prove Theorem 2.11. We argue in a similar manner as in [21], but in contrast our
improved analysis relies on the refined initialization phase developed in Section 4 and its consecutive
parallelization in Section 7. Spielman and Peng’s [21] proof involves two major steps: the first is
to construct a sparse O(1)-approximate inverse chain, and the second is to apply this chain as a
preconditioner into an algorithm known as “Preconditioned Richardson Iteration” [21, Lemma 4.4].
We give now an improved construction for a sparse ε-approximate inverse chain. This directly
implies the desired statement of Theorem 2.11.
c ≈ε/16 D − M = B
Proof of Theorem 2.11. We apply Lemma 7.1 to compute a sparsifier D − M
c) 6 O(ε−2 n logc n) in work O(mB logc1 n) and depth O(logc2 n). Then by Fact 3.1.d
with nnz(M
c)−1 ≈ε/8 (D − M )−1 .
we have (D − M
b = D−M
c.
Our goal now is to construct a sparse ε-approximate inverse chain of the sparsifier B
1+ε/8
′
b satisfies κ b 6
The condition number of matrix B
b ). Spielman
b . Let ε = ε/(16 log tB
B
1−ε/8 κB = tB
and Peng [21] proved that O(log tBb ) iterations suffice for the following iterative procedure to output
a sparse ε-approximate inverse chain.
c2 ≈ε/8 D − M
cD −1 M
c with nnz(M
c2 ) 6
By Lemma 7.2 we compute a spectral sparsifier D − M
c
c+c
+1
c
O(ε−2 n log n) in work O(ε−4 n log 1 n) and depth O(log 2 n). By Lemma 3.8 and Lemma 7.3
c2k with the at most
c2k D −1 M
c2k+1 ≈ε′ D − M
for each consecutive call we compute a sparsifier D − M
c
2
c+c
+1
2
O(ε−2 n log n · log κB ) non-zero entries in work O(ε−4 n log 1 n · log κB ) and depth O(logc2 n).
We combine now Fact 3.1 and apply recursively O(log tBb ) times the relation

≈ε′
=

c2k D −1 )]/2
c2k+1 ]−1 (I + M
c2k )[D − M
[D −1 + (I + D −1 M
c2k D −1 )]/2
c2k ]−1 (I + M
c2k D −1 M
c2k )[D − M
[D −1 + (I + D −1 M
c2k )−1 .
(D − M

c k D −1 M
c k can be replaced with D for
Spielman and Peng showed in [21, Corollary 5.5] that D − M
2
2

k = O(log tBb ) and maintain the desired approximation. The statement follows by Fact 3.1.

9

Representational Power of MDBD

In this section we prove Theorem 2.1. Our analysis relies on the following three influential works.
Hald [12] analyzed mixed Binomial distributions in continuous case. Cruz-Uribe and Neugebauer [7,
8] gave sharp guarantees for approximating integrals using the Trapezoid method. Doha et al. [9]
proved a simple closed formula for higher order derivatives of Bernstein basis.
14

Our goal now is to prove Theorem Theorem 2.1. Hald [12] proved the following result on mixed
Binomial distributions.
Theorem 9.1. [12] Let w(x) be a probability density function that is four times differentiable.
Then for every N ∈ N and i ∈ [0, N ] the Bernstein basis BN,i (p) satisfies



Z 1
w(i/N )
1
b1 (i/N ) b2 (i/N )
w(p) · BN,i (p)dp =
· 1+
+
+O
(3)
N
N
N2
N3
0
1
[−w(x) + (1 − 2x)w′ (x) + 12 x(1 − x)w′′ (x)] and
where the functions are defined by b1 (x) = w(x)
1
[w(x)−3(1−2x)w′ (x)+(1−6x+6x2 )w′′ (x)+ 56 x(1−x)(1−2x)w′′′ (x)+ 81 x2 (1−x)2 w′v (x)].
b2 (x) = w(x)

We distinguish two types of approximation errors. The error term (1 + ηi ) (c.f. Equation 1 in
Theorem 2.1) is caused by the error introduced in Equation 3. The second error type is due to
the integral discretization with finite summation. The later approximation error in analyzed by
Cruz-Uribe and Neugebauer [7, 8]. We summarize below their result.
Theorem 9.2. [7, 8] Suppose f be continuous and twice differentiable function,
PT −1 T ∈ N is number,
1
and the discrete approximator of f is defined by AT (f ) = [ 2 (f (0) + f (1)) + i=1 f (i/T )]/T . Then
R1
PT
the approximation
error
is
given
by
the
expression
E
(f
)
=
|A
(f
)
−
f
(t)dt|
=
|
T
T
i=1 Li |, where
0
R
 1
 ′′
1 xi
2
Li = 2 xi−1 4T 2 − (t − ci ) f (t)dt and ci = (xi−1 + xi )/2.
The rest of this section is devoted to prove Theorem 2.1. We use the following two results
established by Cruz-Uribe and Neugebauer, and Doha et al.
R1
Lemma 9.3. [7, 8] The Bernstein basis satisfies 0 BN,i (x)dx = N 1+1 for every i ∈ [0 : N ].

Lemma 9.4. [9] The pth derivative of a Bernstein basis satisfies for every i ∈ [0 : N ] that
dp BN,i (x)
N!
=
p
dx
(N − p)!

min{i,p}

X

k+p

(−1)

k=max{0,i+p−N }

 
p
·
· BN −p,i−k (x).
k

We propose an upper bound on the integral of pth order derivative of Bernstein basis.
Corollary 9.5. For every p ∈ [1 : N − 1] and i ∈ [p + 1 : N − (p + 1)] such that i + p 6 N , the
Bernstein basis satisfies
Z 1 p
d BN,i (x)
2p
N!
(t)
·
.
dt
6
dxp
(N − p)! N + 1
0
Proof. We combine Lemma 9.3 and Lemma 9.4 to obtain
Z 1 p
p   Z 1
X
d BN,i (x)
p
N!
2p
N!
B
(t)dt
=
·
(t)
·
.
dt
6
N
−p,i−k
k
dxp
(N − p)!
(N − p)! N + 1
0
0
k=0



We are now ready to prove Theorem 2.1.
Proof of Theorem 2.1. Recall that Fi (x) = w(x)BN,i (x). By Theorem 9.2 we have
T
X
k=1

T

Lk

1X
=
2

Z

xk

k=1 xk−1



 2
Z 1 2
d Fi (x)
1
1
d Fi (x)
2
6
−
(t
−
c
)
(t)dt
(t) dt.
·
k
2
2
2
4T
dx
8T 0
dx2
15

2

Fi (x)
′
′′ , we consider following three cases:
| = w′′ · BN,i + 2 · w′ · BN,i
+ w · BN,i
Since | d dx
2
′′
Case 1: We combine maxx∈[0,1] |w (x)| 6 2φw · N 2 and Lemma 9.3 to obtain

Z

1

0

|w′′ (t) · BN,i (t)|dt 6 2φw · N.

Case 2: Using maxx∈[0,1] |w′ (x)| 6 12 φw · N and Corollary 9.5 it holds
Z

0

1

′
|w′ (t) · BN,i
(t)|dt 6 φw · N.

Case 3: Combining maxx∈[0,1] |w(x)| 6 φw and Corollary 9.5 yields
Z

1
0

′′
|w(t) · BN,i
(t)|dt 6 φw ·

Z

1
0

′′
BN,i
(t)dt 6 4φw · N.



The desired result follows from the preceding three cases and Theorem 9.1.

10

Approximating Discretized PDF

In this section we prove Theorem 2.2. We begin our discussion by presenting the pseudo code of
algorithm AppDscrPDF.
Algorithm 3 Approximate Discretized PDF by MDBD
BN,T (α, p) = AppDscrPDF(w, φw , N, εI )
p
P
PT
φw /εI ⌉.
1. Compute SN +1,N = N
i=0 w(i/N ) and ST,T +1 =
i=1 w(i/[T + 1]), where T = ⌈N
2. Compute αj = w(pj ) · N/[(T + 1) · SN +1,N ], where pj = j/[T + 1] for all j ∈ [1 : T ].
3. Return BN,T (α, p).
Before we prove Theorem 2.2, we analyze the class of continuous probability density functions
that admit a discretized approximation by MDBD.
Lemma 10.1. Let w(x) = C · f (x) be a twice differentiable p.d.f. such that
R 1 for N ∈ N+ it holds
a) 0 6 f (x) 6 1, b) 12 [f (0) + f (1)] > Ω(1), c) 1 6 C 6 o(N ), and d) 0 f (2) (x) dx 6 o(N ).
p
Then, for T = ⌈N φw /εI ⌉ with φw /εI > 1 it holds
1−

T

N

k=1

k=0

X
X
ST,T +1 /(T + 1)
w(k/N ).
w(k/[T + 1]) and SN +1,N ,
= o(1), where ST,T +1 ,
SN +1,N /N

Proof. By Theorem 9.2 for the discrete approximator of f
"
#
M
−1 
X
i
1 1
f
(f (0) + f (1)) +
,
AM (f ) =
M 2
M
i=1

it holds that
AT +1 (f ) −

Z

0

1

1
f (t)dt 6
8 (T + 1)2
16

Z

0

1

f

(2)

εI
·o
(t) dt .
φw



1
N



,

and similarly
AN (f ) −
By definition,

R1

Z

1
0

1
f (t)dt 6
8N 2

Z

1

f

(2)

(t) dt . o

0

f (t)dt = C −1 ∈ (1/o(N ), 1] and thus
 

1
εI
1
±
·o
, and
AT +1 (f ) ∈
C
φw
N



1
N



.

0



1
±o
AN (f ) ∈
C



1
N



.

Let d , [f (0) + f (1)]/2. Straightforward checking shows that




ST,T +1
SN +1,N
d
d
= C · AT +1 (f ) −
= C · AN (f ) +
and
.
T +1
T +1
N
N
We prove now the upper bound. By assumption d ∈ [Ω(1), 1] and since C 6 o(N ) we have

εI
1
1
d
d
AT +1 (f ) − T +1
ST,T +1 /(T + 1)
C − T +1 ± φw · o N

=
=
Λ,
d
1
1
SN +1,N /N
AN (f ) + Nd
C + N ±o N
p 


1 + 21 εkI · Nd − 1 + εkI · o N1
1
Ω(1)
d

· 1
6 1−
61−
1 = 1 − o(1).
1
1
N
C + N
C + N −o N
We can prove the lower bound Λ > 1 − o(1) using similar arguments.



We present now the proof of Theorem 2.2.
Proof of Theorem 2.2. By definition ST,T +1 =
desired discretized p.d.f. is

PT

k=1 w(j/[T

+ 1]), SN +1,N =

PN

j=0 w(j/N )

and the

w(i/N )
w(i/N )
=
.
w(i/N
b
) = PN
SN +1,N
j=0 w(j/N )

By Theorem 2.1, it holds for all i ∈ [3 : N − 3] that
T

(1 + ηi )

εI
w(i/N ) X w(j/[T + 1])
−
· BN,i (j/[T + 1]) 6 .
N
T +1
N

(4)

j=1

We construct now MDBD BN,T (α, p) as follows: for all j ∈ [1 : T ] we set
αj =

w(pj ) · N
,
SN +1,N · (T + 1)

where pj =

Moreover, BN,T (α, p) induces a vector γ that satisfies γi =
multiplying Equation 4 with N/SN +1,N we obtain

Furthermore, since
T
X
j=1

αj =

j
.
T +1

PT

j=1 αj

· BN,i (pj ) for all i ∈ [0 : N ]. By

|(1 + ηi ) w(i/N
b
) − γi | 6 εI /SN +1,N .

(5)

T
X
ST,T +1 /(T + 1)
N
w(pj ) =
,
(T + 1) · SN +1,N
SN +1,N /N
j=1

17

by Lemma 10.1 there is a small positive number δw = o(1) such that
T

δw = 1 −

X
ST,T +1 /(T + 1)
αj .
=1−
SN +1,N /N
j=1

Hence, we have
N
X
i=0

γi =

T
N X
X
i=0 j=1

αj · BN,i (pj ) =

T
X
j=1

αj

N
X
i=0

BN,i (pj ) =

T
X
j=1

αj = 1 − δw .

Since δw = o(1), by Equation 5 it follows that γ/[1 − δw ] is a discretized probability distribution
over [0 : N ] that approximates component-wise the desired discretized p.d.f. w.
b
We note that the summations SN +1,N and ST,T +1 can be computed in O(N + T ) work and
O(log N + log T ) depth. Hence, the statement follows.


11

Efficient Parallel Solver for Transpose Bernstein-Vandermonde
Systems

Problem 1. Suppose a vector γ ∈ (0, 1)N +1 is induced by a convex combination of exactly N + 1
discrete Binomial distributions B(pi , N ) such that 0 < pi 6= pj < 1 for all i 6= j. Find the unique
P +1
vector α ∈ (0, 1)N +1 such that γi = N
j=1 αj · BN,i (pj ) for all i ∈ [0 : N ].
The Bernstein basis is a well studied
primitive in the literature for polynomial interpolations [7,

N k
8]. It is defined by BN,k (p) = k p (1 − p)N −k for any k ∈ [0 : N ]. Let BN,T (α, p) be MDBD with
T = N + 1. Then the Bernstein basis matrix is defined by [BN (p)]ji = BN,i (pj ), ∀i, j ∈ [1 : N + 1],
and it has a full rank (c.f. Appendix C). Moreover, the vector α is the unique solution of the linear
system BN (p)T α = γ.
In this section, we give an efficient parallel algorithm that solves Problem 1 and works in nearly
linear work and poly-logarithmic depth. Our goal now is to prove Theorem 2.3. We reduce a
transpose Bernstein-Vandermonde system to a transpose Vandermonde system that can be solved
efficiently and in parallel by a variation of an algorithm proposed by Gohberg and Olshevsky [11].
We present now their main algorithmic result.
Theorem 11.1. [11] There is an algorithm that on input two vectors α, p ∈ RN +1 such that
0 < pi 6= pj < 1 for all i 6= j, outputs the vector γ = V(p)T α in O(N log2 N ) time.
Theorem 11.1 follows by [11, Algorithm 2.1] which relies on a non-trivial matrix decomposition of
V(p)T to compute in O(N log2 N ) time the desired matrix-vector product. It can be easily verified
that [11, Algorithm 2.1] can be amended to compute the vector α = [V(p)T ]−1 γ in O(N log2 N )
time. Furthermore, straightforward checking shows that this modified algorithm can be easily
parallelized. We summarize below the resulting parallel algorithm.
Theorem 11.2. [11] There is a parallel algorithm that on input two vectors γ, p ∈ RN +1 as in
Theorem 11.1, outputs the vector α = [V(p)T ]−1 γ in O(N log2 N ) work and O(logc n) depth, for
some constant c ∈ N+ .
We prove now that the Bernstein basis matrix BN (p) admits the following decomposition.

18

Lemma 11.3. Suppose p ∈ (0, 1)N +1 is vector such that 0 < pi 6= pj < 1 for all i 6= j, V(p)
p
+1
is Vandermonde matrix defined by [V(p)]ji = ( 1−pj j )i , Dp = diag({(1 − pj )N }N
j=1 ) and DCN =

diag({ Ni }N
i=0 ) are positive diagonal matrices. Then it holds that BN (p) = Dp · V(p) · DCN .

p
Proof. By definition [Dp · V(p) · DCN ]j,i = (1 − pj )N ( 1−pj j )i Ni = BN,i (pj ) = [BN (p)]ji .

We are ready now to prove Theorem 2.3.

Proof of Theorem 2.3. By Lemma C.1 the Bernstein matrix BN (p) is invertible. Given a vector
γ ∈ (0, 1)N +1 we want to find the vector α = [BN (p)T ]−1 γ. By Lemma 11.3 we have [BN (p)T ]−1 =
[Dp ]−1 · [V(p)T ]−1 · [DCN ]−1 . Moreover, we can compute a vector γ ′ = [DCN ]−1 γ in O(n) time.
Using Theorem 11.2, we obtain a vector γ ′′ = [V(p)T ]−1 γ ′ in O(N log2 N ) work and O(logc n)
depth. The desired vector α = [Dp ]−1 γ ′′ takes further O(n) work and O(1) depth to compute. 
Acknowledgements We are grateful to Arijit Ghosh and Kunal Dutta for the helpful discussions on mixture of Binomial distributions, and to Arnur Nigmetov and Shay Moran for pointing
us to the Bernstein interpolation polynomials. We would also like to thank Kurt Mehlhorn for the
insightful comments and suggestions to the early version of the manuscript.
This work has been funded by the Cluster of Excellence “Multimodal Computing and Interaction” within the Excellence Initiative of the German Federal Government.

References
[1] R. Andersen, F. R. K. Chung, and K. J. Lang. Local graph partitioning using pagerank vectors. In 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS 2006),
21-24 October 2006, Berkeley, California, USA, Proceedings, pages 475–486, 2006. URL
http://dx.doi.org/10.1109/FOCS.2006.44.
[2] M. W. Bern, J. R. Gilbert, B. Hendrickson, N. Nguyen, and S. Toledo. Support-graph preconditioners. SIAM Journal on Matrix Analysis and Applications, 27(4):930–951, 2006. doi:
10.1137/S0895479801384019. URL http://dx.doi.org/10.1137/S0895479801384019.
[3] D. Cheng, Y. Cheng, Y. Liu, R. Peng, and S. Teng. Scalable parallel factorizations of SDD
matrices and efficient sampling for gaussian graphical models. CoRR, abs/1410.5392, 2014.
URL http://arxiv.org/abs/1410.5392.
[4] D. Cheng, Y. Cheng, Y. Liu, R. Peng, and S. Teng. Efficient sampling for gaussian
graphical models via spectral sparsification. In Proceedings of The 28th Conference on
Learning Theory, COLT 2015, Paris, France, July 3-6, 2015, pages 364–390, 2015. URL
http://jmlr.org/proceedings/papers/v40/Cheng15.html.
[5] D. Cheng, Y. Cheng, Y. Liu, R. Peng, and S. Teng. Spectral sparsification of random-walk
matrix polynomials. CoRR, abs/1502.03496, 2015. URL http://arxiv.org/abs/1502.03496.
[6] M. B. Cohen, Y. T. Lee, C. Musco, C. Musco, R. Peng, and A. Sidford. Uniform sampling for
matrix approximation. In Proceedings of the 2015 Conference on Innovations in Theoretical
Computer Science, ITCS 2015, Rehovot, Israel, January 11-13, 2015, pages 181–190, 2015.
doi: 10.1145/2688073.2688113. URL http://doi.acm.org/10.1145/2688073.2688113.

19

[7] D. Cruz-Uribe and C. Neugebauer.
Sharp error bounds for the trapezoidal rule
and Simpson’s rule.
JIPAM. J. Inequal. Pure Appl. Math., 3(49), 2002.
URL
http://www.emis.de/journals/JIPAM/images/031_02_JIPAM/031_02.pdf.
[8] D. Cruz-Uribe and C. J. Neugebauer. An elementary proof of error estimates for the
trapezoidal rule. Mathematics Magazine, 76(4):pp. 303–306, 2003. ISSN 0025570X. URL
http://www.jstor.org/stable/3219088.
[9] E. Doha, A. Bhrawy, and M. Saker. On the derivatives of Bernstein polynomials: an application
for the solution of high even-order differential equations. Bound. Value Probl., 2011:16, 2011.
ISSN 1687-2770/e. doi: 10.1155/2011/829543.
[10] S. O. Gharan and L. Trevisan. Approximating the expansion profile and almost optimal local
graph clustering. In 53rd Annual IEEE Symposium on Foundations of Computer Science,
FOCS 2012, New Brunswick, NJ, USA, October 20-23, 2012, pages 187–196, 2012. URL
http://dx.doi.org/10.1109/FOCS.2012.85.
[11] I. Gohberg and V. Olshevsky. Fast algorithms with preprocessing for matrix-vector multiplication problems. Journal of Complexity, 10(4):411 – 427, 1994. ISSN 0885-064X. URL
http://www.sciencedirect.com/science/article/pii/S0885064X84710211.
[12] A. Hald. The mixed binomial distributin and the posterior distribution of p for a continuous
prior distribution. Journal of the Royal Statistical Society. Series B (Methodological), 30(2):
pp. 359–367, 1968. ISSN 00359246. URL http://www.jstor.org/stable/2984516.
[13] J. A. Kelner and A. Levin. Spectral sparsification in the semi-streaming setting. Theory Comput. Syst., 53(2):243–262, 2013. URL http://dx.doi.org/10.1007/s00224-012-9396-1.
[14] P. Kolev and H. Sun. Dirichlet eigenvalues, local random walks, and analyzing clusters in graphs. In Algorithms and Computation - 25th International Symposium, ISAAC
2014, Jeonju, Korea, December 15-17, 2014, Proceedings, pages 621–632, 2014. URL
http://dx.doi.org/10.1007/978-3-319-13075-0_49.
[15] I. Koutis, G. L. Miller, and R. Peng. A nearly-m log n time solver for sdd linear systems. In
Proceedings of the 2011 IEEE 52Nd Annual Symposium on Foundations of Computer Science,
FOCS ’11, pages 590–598, Washington, DC, USA, 2011. IEEE Computer Society. ISBN 9780-7695-4571-4. URL http://dx.doi.org/10.1109/FOCS.2011.85.
[16] I. Koutis, A. Levin, and R. Peng. Improved spectral sparsification and numerical algorithms for
SDD matrices. In 29th International Symposium on Theoretical Aspects of Computer Science,
STACS 2012, February 29th - March 3rd, 2012, Paris, France, pages 266–277, 2012. URL
http://dx.doi.org/10.4230/LIPIcs.STACS.2012.266.
[17] R. Kyng, Y. T. Lee, R. Peng, S. Sachdeva, and D. A. Spielman. Sparsified cholesky
and multigrid solvers for connection laplacians. CoRR, abs/1512.01892, 2015. URL
http://arxiv.org/abs/1512.01892.
[18] G. L. Miller and R. Peng. Approximate maximum flow on separable undirected graphs. In
Proceedings of the Twenty-Fourth Annual ACM-SIAM Symposium on Discrete Algorithms,
SODA 2013, New Orleans, Louisiana, USA, January 6-8, 2013, pages 1151–1170, 2013. URL
http://dx.doi.org/10.1137/1.9781611973105.83.

20

[19] H. Minc. Nonnegative Matrices. John Wiley and Sons, New York, 1973.
[20] L. Orecchia and N. K. Vishnoi. Towards an sdp-based approach to spectral methods: A
nearly-linear-time algorithm for graph partitioning and decomposition. In Proceedings of the
Twenty-second Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’11, pages 532–
545. SIAM, 2011. URL http://dl.acm.org/citation.cfm?id=2133036.2133078.
[21] R. Peng and D. A. Spielman. An efficient parallel solver for SDD linear systems. In Symposium
on Theory of Computing, STOC 2014, New York, NY, USA, May 31 - June 03, 2014, pages
333–342, 2014. URL http://doi.acm.org/10.1145/2591796.2591832.
[22] R. Y. Peng. Algorithm Design Using Spectral Graph Theory. PhD thesis, Carnegie Mellon
University, Pittsburgh, August 2013. CMU CS Tech Report CMU-CS-13-121.
[23] A. Sinclair and M. Jerrum.
Approximate counting, uniform
rapidly mixing markov chains.
Inf. Comput., 82(1):93–133,
http://dx.doi.org/10.1016/0890-5401(89)90067-9.

generation and
1989.
URL

[24] D. A. Spielman and N. Srivastava.
Graph sparsification by effective resistances.
In Proceedings of the 40th Annual ACM Symposium on Theory of Computing, Victoria, British Columbia, Canada, May 17-20, 2008, pages 563–568, 2008.
URL
http://doi.acm.org/10.1145/1374376.1374456.
[25] D. A. Spielman and S. Teng. Spectral sparsification of graphs. SIAM J. Comput., 40(4):
981–1025, 2011. doi: 10.1137/08074489X. URL http://dx.doi.org/10.1137/08074489X.
[26] D. A. Spielman and S.-H. Teng. Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems. SIAM Journal on Matrix Analysis and
Applications, 35(3):835–885, 2014. URL http://dx.doi.org/10.1137/090771430.

A

Generalized Escaping Probability

b preserves
Proof of Lemma 2.10. By definition, for every vector x the spectral sparsifier D − A
approximately the quadratic form
b ∈ [(1 ± ε) · xT (D − DGγ )x].
xT (D − A)x

Hence, the statement follows by applying the identities
ξST (D − DGγ )ξS

= πST (I − Gγ )1S = 1 − πST Gγ 1S = πST Gγ 1S
= Ev∼πS [gEsc(v, S, Gγ )] .



B

Spectral Sparsification of T -Matrices

Our proof of Lemma B.2 is based on the Perron-Frobenius Theorem [19] for non-negative matrices.
Theorem B.1. [19, Perron-Frobenius] Suppose A is symmetric nonnegative matrix. Then it has
a nonnegative eigenvalue λ which is greater than or equal to the modulus of all other eigenvalues.
21

Lemma B.2. Suppose D−M is SDDM matrix. Then D−D(D −1 M )N is SDDM matrix ∀N ∈ N+ .
Proof. Since D − M is SDDM we have D ≻ M . By Fact 3.1.e it holds I ≻ D −1/2 M D −1/2 , X.
radius ρ(X) < 1, i.e.
Hence, the largest eigenvalue λ(X) < 1. By Theorem B.1 the spectral
P
|λi (X)|P< 1 for all i. Since X is symmetric it has the form X = i λi ui uT
i . Moreover, we have
k
k
T
k
X = i λi ui ui for every k ∈ N+ . Thus the spectral radius of X satisfies ρ(X k ) = ρ(X)k < 1.
Notice that Bk = D − D(D −1 M )k is symmetric non-negative matrix for every k ∈ N+ . By
definition D − M is diagonally dominant and thus D −1 M 1  1 component-wise. This implies that
Bk is diagonally dominant matrix. Notice that Bk = D 1/2 [I − X k ]D 1/2 and since ρ(X k ) < 1 it
follows that Bk is positive definite and hence SDDM matrix.

Proof of Lemma 3.8 We combine Lemma B.2 with the following two statements.
Fact B.3. Spielman and Peng [21, Proposition 5.6] showed that if D − M is SDDM matrix, then
D − M D −1 M is SDDM matrix. Also Cheng et al. [4, Proposition 25] showed that if D − M is
Laplacian matrix then D − D(D −1 M )N is Laplacian matrix for every N ∈ N+ .
Based on Lemma 3.9 and Fact B.3 we establish the following result.
c ≈ε D − M is a spectral sparsifier. Then
Lemma B.4. Suppose D − M is T -matrix and D − M
−1
N
c
D − D(D M ) is T -matrix for every N ∈ N+ .
Proof of Lemma 3.7 We use the following result that appears in Peng’s thesis [22].

e −A
e a sparsifier
Lemma B.5. [22] Suppose D − A is Laplacian matrix (possibly Aii 6= 0), and D
e
e
e
with Aii = 0 for every i such that (1 − ε)(D − A)  D − A  D − A. Then the symmetric
b = (D − D)
e +A
e satisfies (1 − ε)(D − A)  D − A
b  (D − A).
non-negative matrix A

f1 − A
f1  D − A
− A)  D
b = (D − D
f1 ) + A
f1 satisfies 1−ε (D − A) 
and by Lemma B.5 the symmetric non-negative matrix A
1+ε
b  D − A. Since 1−ε > 1 − 2ε for every ε ∈ (0, 1 ) the statement follows.
D−A

1+ε
2
f1 =
Proof of Lemma 3.7. Let D

B.1

1 e
1+ε D

f1 =
and A

1 e
1+ε A.

Then

1−ε
1+ε (D

Structural Result

Suppose D − M is T -matrix. We show that the matrix D − M D −1 M can be expressed as a sum
of a non-negative main diagonal matrix and a sum of Laplacian matrices.
Proof of Lemma 3.10. Let M ∈ Rn×n and nnz(M ) = m. We decomposePthe entries of matrix
2 /D for all i.
M D −1 M into three types. We set type 1 to be the entries (M D −1 M )ii = nk=1 Mik
k
We note that all entries of type 1 can be computed in O(m) time. We consider next the off-diagonal
entries
n
X
Mik Mjk /Dkk
(Mii /di + Mjj /dj ) · Mij +
−1
{z
}
|
.
(M D M )ij =
k6={i,j}
|
{z
}
type 2
type 3

Observe that the number of type 2 entries is at most m. Now for a fixed k we note that the
corresponding entries that appear in type 1 and type 3 form a weighted clique (with self-loops)
whose adjacency matrix is defined by d1k ηk ηkT .
P
Straightforward checking shows that M D −1 M = B + i d1i ηi ηiT . By Lemma 3.8 D − M D −1 M
is T -matrix and thus diagonally dominant. Hence, the Laplacian matrices LB and LNi for all i
22

exist. Moreover, we can compute in O(m) time the positive diagonal matrices DB and DNi =
(si /di )diag(Ni ) for all i. To see this, observe that si and di can be computed in O(m) time for all
i, and the number of elements in the disjoint union | ⊔i Ni | 6 m.


C

Bernstein Basis Matrix

We prove below that the Bernstein basis matrix in Problem 1 has full rank.
Lemma C.1. Suppose a vector p ∈ (0, 1)N +1 satisfies 0 < pi 6= pj < 1 for all i 6= j. Then the
Bernstein basis matrix BN (p) has a full rank.
N +1 such
Proof. Suppose for contradiction that rank(BN (p)) < N + 1. Then
PNthere is a vector λ ∈ R
that the linear combination of the columns of BN (p) satisfies i=0 λj [BN (p)]:,i = 0. Let fλ (x)
 i
P
PN
N
N −i . Notice that
be a polynomial defined by fλ (x) , N
i=0 λi · BN,i (x) =
i=0 λi · i x (1 − x)
fλ (pj ) = 0 for every j ∈ [1 : N + 1], i.e. fλ (x) has N + 1 roots. However, since the polynomial
fλ (x) has degree N it follows that fλ (x) ≡ 0. Therefore, we obtained the desired contradiction. 

D

Approximating Two Canonical PDFs

Here, we illustrate the representational power of MDBD. We show that there are MDBD satisfying
the hypothesis in Theorem 2.2 and approximate two canonical continuous p.d.f.: the Uniform
distribution and the Exponential Families. More precisely, we prove that the Uniform distribution
and the Exponential families admit a multiplicative and an additive approximation, respectively.

D.1

Uniform Distribution

Lemma D.1 (Uniform Distribution). Let w(x) = 1, N ∈ N+ and ε > 0. If T > Ω(N ε−1/2 ) then
it holds that
T

1 X
BN,i
T +1
j=1



j
T +1






1
∈ (1 ± ε)
, for all i ∈ [3 : N − 3].
N +1

Proof. By Theorem 9.2 we have that
T
X

T

Li 6

i=1

1X
2
i=1

Z

xi
xi−1

2
1
1
2 d BN,i (x)
−
(t
−
c
)
(t) dt 6
i
2
2
4N
dx
8N 2

Z

0

1

d2 BN,i (x)
(t) dt.
dx2

By combining Lemma 9.3 and Corollary 9.5 for every i ∈ [3 : N − 3] it holds
T

1
1 X
BN,i
−
N +1 T +1
j=1

We note that

N
2T 2

6

ε
N +1



j
T +1



6

T
X
i=1

1
Li 6
8T 2

Z

0

1

d2 BN,i (x)
N
(x) dx 6
.
2
dx
2T 2

since T > Ω(N ε−1/2 ), and hence the statement follows.

that w(x) = 1 · 1, i.e., f (x) = 1.
Remark D.2. All conditions of Theorem 2.2 hold. Note
R 1 (2)
1
a) C = 1, b) f (x) = 1, c) 2 [f (0) + f (1)] = 1 and d) 0 |f (x)|dx = 0, since f (2) (x) = 0.
23



D.2

Exponential Families

√
Lemma D.3 (Exponential Families). Let N ∈ N, k ∈ [1, N ] and w(x) = 1−ek−k · exp{−k · x} is
p
a probability density function. For any ε > 0 if T > Ω(N k/ε) then for every i ∈ [3 : N − 3] it
holds for the function Fi (x) = w(x) · BN,i (x) that
T

1 X
w(i/N )
Fi
−
(1 + ηi )
N
T +1
j=1



j
T +1



6

ε
,
N

where

1
|ηi | 6 .
4

Proof. The pth derivative of function w(x) satisfies w(p) (x) = (−k)p · w(x). Let I = [0, 1] be an
interval. Straightforward checking shows that
max |w(p) (x)| = kp · max |w(x)| < 2 · kp+1 .
x∈I

x∈I

(6)

By the definition of function b1 (x) (c.f. Theorem 9.1) we have


k2
k2 2
· x − (1 + k),
b1 (x) = − · x + 2k +
2
2
and we can show that maxx∈I b1 (x) 6 1 + k2 /8 6 1 + N/8. The function b2 (x) satisfies
5
1
b2 (x) = 1 + 3(1 − 2x)k + (1 − 6x + 6x2 )k2 − x(1 − x)(1 − 2x)k 3 + x2 (1 − x)2 k4 ,
6
8
and we can show that maxx∈I b2 (x) ≪ k4 /8 6 N 2 /8. By Theorem 2.1 it suffices to upper bound
the following four cases.
Case 1: By Theorem 9.1 µ 6 14 , since maxx∈I |b1 (x)| 6 1 + N/8 and maxx∈I |b2 (x)| ≪ N 2 /8.
Case 2: We combine Equation 6 and Lemma 9.3 to obtain
Z 1
Z 1
2 · k3
|BN,i (t)|dt =
.
|w′′ (t) · BN,i (t)|dt < 2 · k3 ·
N +1
0
0
Case 3: By combining Equation 6 and Lemma 9.5 it holds
Z 1
Z 1
′
′
′
2
|BN,i
(t)|dt < 4 · k2 .
|w (t) · BN,i (t)|dt < 2 · k ·
0

0

Case 4: We use again Equation 6 and Lemma 9.5 to obtain
Z 1
Z 1
′′
′′
|BN,i
(t)|dt < 8 · k · N.
|w(t) · BN,i (t)|dt < 2 · k ·
0

0

2

Fi (x)
′
′′ and
Recall that Fi (x) = w(x)BN,i (x). By combining | d dx
| = w′′ · BN,i + 2w′ · BN,i
+ w · BN,i
2
Theorem 9.2 we have
Z 1 2
T
X
d F (x)
k·N
1
(t) dt .
.
Li 6
2
2
8T 0
dx
T2
i=1



Hence, the statement follows.

24

Remark D.4. The hypothesis in Theorem 2.2 holds. Note that w(x) = 1−ek−k · exp{−k · x} and
√
φw = 1−ek−k . Thus C = 1−ek−k , f (x) = exp{−k · x}, f (2) (x) = k2 · f (x) and k ∈ [1, N ].
Furthermore,
a) C = o(N ), b) 0 6 f (x) 6 1, c) 12 6 12 [f (0) + f (1)] < 1 and for d) we have
Z

0

E

1

|f

(2)

(x)|dx = k

2

Z

1
0

|f (x)|dx =

k2
= (1 − e−k )k = o(N ).
C

Schur Complement

In this section we prove Lemma 4.1. We use the following result proposed by Peng et al. [5].
c  (1+ε)(D −M )
Lemma E.1. [5, Lemma 4.3] If M is SPSD matrix and (1−ε)(D −M )  D − M
c  (1 + ε)(D + M ).
then it holds that (1 − ε)(D + M )  D + M

We extend next two technical results on Schur complement that appeared in [4, 18, 22].


P1 −M
where P1 and P2 are symmetric positive definite matrices
Claim E.2. Suppose X ,
−M P2
 T  
u
u
−1
T
and M is symmetric matrix. Then v [P2 − M P1 M ]v = minu
X
for every v.
v
v


T 
 
u
P1 −M
u
Proof. Suppose fv (u) ,
= uT P1 u − 2uT M v + v T P2 v. Notice
v
−M P2
v
that f is minimized when u = P1−1 M v, since ∇fv (u) = 2P1 u − 2M v. Hence, it follows that
 T  
u
u
X
minu
= v T P2 v − v T M P1−1 M v = v T [P2 − M P1−1 M ]v.

v
v
Lemma E.3. (Schur Complement) Suppose
D1 , D2are positive
main diagonal matrices and M, Q


D1 −M
D1 −Q
are symmetric matrices. If XM ,
≈ε
, XQ then it holds that
−M D2
−Q D2
D2 − M D1−1 M ≈ε D2 − QD1−1 Q.
Proof. Here we show the upper bound, but the lower bound follows by analogy. Let w be a vector

 

 T

T
w
u
u
w
such that
XQ
XQ
. We apply twice Claim E.2 to obtain the
= minu
v
v
v
v
following chain of inequalities
T

v [D2 −

M D1−1 M ]v

= min
u



6 (1 + ε)

u
v


T
w
v

XM
T



XQ

u
v



w
XM
6
v

 T
 
w
u
u
= (1 + ε) min
XQ
u
v
v
v





w
v

T



= (1 + ε)v T [D2 − QD1−1 Q]v.



We prove Lemma 4.1 by arguing in a similar manner to in [5, Lemma 4.4]. We present the
proof here for completeness.

25




D −X
Proof of Lemma 4.1. For any symmetric matrix X, we denote by PX =
. We prove
−X D
the upper bound, but the lower bound follows by analogy. Straightforward checking shows that


T
 
u
u
cu − uT M
cv + v T Dv
PM
= uT Du − v T M
c
v
v
i
1h
c)(u + v) + (u − v)T (D + M
c)(u − v)
(u + v)T (D − M
2

=

c ≈ε D + M and thus we have
By Lemma E.1 it holds that D + M


u
v

T

PM
c



u
v



i
1h
c)(u + v) + (u − v)T (D + M
c)(u − v)
(u + v)T (D − M
2

1
6 (1 + ε) (u + v)T (D − M )(u + v) + (u − v)T (D2 + M )(u − v)
2
 T
 
u
u
PM
= (1 + ε)
.
v
v
=

Hence, it follows that PM
c =

"

c
D −M
c
−M
D

cD −1 M
c ≈ε D − M D −1 M .
holds that D − M

#

≈ε



26

D
−M

−M
D



= PM . Now, by Lemma E.3 it


