arXiv:1608.01426v1 [cs.CC] 4 Aug 2016

Solving Laplacian Systems in Logarithmic Space
François Le Gall
Graduate School of Informatics
Kyoto University
legall@i.u-kyoto.ac.jp

Abstract
We investigate the space complexity of solving linear systems of equations. While all
known deterministic or randomized algorithms solving a square system of n linear equations
in n variables require Ω(log2 n) space, Ta-Shma (STOC 2013) recently showed that on a quantum computer an approximate solution can be computed in logarithmic space, giving the first
explicit computational task for which quantum computation seems to outperform classical
computation with respect to space complexity. In this paper we show that for systems of
linear equations in the Laplacian matrix of graphs, the same logarithmic space complexity
can actually be achieved by a classical (i.e., non-quantum) algorithm. More precisely, given
a system of linear equations Lx = b, where L is the (normalized) Laplacian matrix of a
graph on n vertices and b is a unit-norm vector, our algorithm outputs a vector x̃ such that
k x̃ − x k ≤ 1/poly(n) and uses only O(log n) space if the underlying graph has polynomially
bounded weights. We also show how to estimate, again in logarithmic space, the smallest
non-zero eigenvalue of L.

1 Introduction
Background. Ta-Shma showed a few years ago that several fundamental tasks in linear algebra,
such as computing the eigenvalues of an n × n matrix A or finding a solution to a linear system
of equations Ax = b, can be solved approximately with polynomial precision on a quantum computer using O(log n) space [16]. In comparaison, no o(log2 n)-space classical (i.e., non-quantum)
algorithm is known for these problems: the best (in the space complexity setting) known classical
algorithms that solve them exactly use O(log2 n) space [1, 2, 5] and nothing better is known for
only approximating the solutions. This breakthrough was the first example of concrete computational tasks for which a quantum algorithm outperforms the best known classical algorithms
in the standard space complexity setting. The power of space-bounded quantum algorithms for
matrix problems comes from their ability to represent an n-dimensional vector using O(log n)
quantum bits of memory and consequently estimate the entries of Ak , given an n × n matrix A
and an integer k ≤ poly(n), in logarithmic quantum space. In comparison, the best known
algorithm for this task uses O(log2 n) space [2].
Doron and Ta-Shma [6] recently showed how to “dequantize” the algorithm [16] for classes of
matrices where Ak can be computed space efficiently, and in particular when A is the adjacency
matrix of a graph, in which case entries of Ak can be estimated by a classical algorithm running k
steps of the random walk corresponding to k. They focused on one task, computation of the
eigenvalues, and succeeded in constructing a O(log n)-space classical algorithm that computes
the eigenvalues of such a matrix A, but only with constant precision (a much weaker accuracy
that the polynomial precision obtained in the quantum case).
Our results. In this paper we investigate the classical (i.e., non-quantum) space complexity of
matrix problems for matrices associated with graphs. More precisely, we focus on the Laplacian
matrices of undirected weighted graphs (defined in Section 2).
The following definition will be useful to state our results concisely: We say that an undirected
weighted graph G on n vertices has polynomially bounded weights if the weight of any edge is
upper bounded by a polynomial of n and lower bounded by a polynomial of 1/n. Unweighted
graphs, where the weight of each edge is one, are examples of such graphs.
Our first result is stated in the following theorem.
Theorem 1.1. Let G be an undirected weighted graph on n vertices with polynomially bounded weights.
For any ǫ, γ ∈ (0, 1), an ǫ-additive approximation of the solution of a Laplacian system corresponding
to G can be computed with probability at least 1 − γ in O(log(n/ǫ) + log log(1/γ)) space.
Theorem 1.1 shows that if G has polynomially bounded weights, linear systems of equations
involving the Laplacian matrix of G can be solved approximately with polynomial precision in
O(log n) space, which gives a classical algorithm with the same space complexity and the same
precision as Ta-Shma’s quantum algorithm [16]. Note that Laplacian systems are a natural and
very well-studied subclass of linear systems, and have a multitude of algorithmic applications
(see [17] for a survey). While in the time complexity setting extremely fast algorithms for Laplacian systems have been obtained in the past decade [8, 11, 12, 13, 14, 15], we are not aware of any
work on the space complexity of this problem.
Roughly speaking, Theorem 1.1 is proved by showing that solving a Laplacian system reduces
to computing powers of the normalized adjacency matrix of the graph, and then computing
2

space-efficiently the entries of the powers by using random walks. We note that Chung and
Simpson [4] considered a fairly similar approach based on random walks to design time-efficient
algorithms solving Laplacian systems under boundary conditions. In the present paper the focus
is different since we consider space complexity. Indeed, the main difficulty here is to show that
the reduction to random walks can be implemented space efficiently.
The values of the eigenvalues of the Laplacian of a graph are related to many important
properties of the graph — these relations are precisely the main subject of the field of spectral
graph theory. The smallest non-zero eigenvalue is especially a fundamental quantity. If the graph
is connected, then this eigenvalue is called the algebraic connectivity of the graph and controls for
example expanding properties of the graph. Our second result, stated in the following theorem,
shows that the smallest non-zero eigenvalue of the Laplacian matrix of an undirected weighted
graph with polynomially bounded weights can be approximated with constant multiplicative
precision in logarithmic space.
Theorem 1.2. Let G be an undirected weighted graph on n vertices with polynomially bounded weights.
For any δ, γ ∈ (0, 1), a δ-multiplicative approximation of the smallest non-zero eigenvalue of the Laplacian
matrix of G can be computed with probability at least 1 − γ in O( 1δ log n + log log(1/γ)) space.
This result is obtained by developing a space-efficient version of a well known technique,
called the power method and used for estimating the largest eigenvalue of a matrix, and showing
how to implementing it space-efficiently using random walks.
To our knowledge the only related prior work on the space complexity of approximating
the algebraic connectivity of a graph is the recent work by Doron and Ta-Shma [6] mentioned
above. Their techniques can be used to estimate the algebraic connectivity in logarithmic space
with constant additive precision. Our techniques give constant multiplicative precision, which
is a much stronger result (the smallest non-zero eigenvalue is always smaller than two; in many
applications this eigenvalue is actually close to zero).

2 Preliminaries
General notations. In this paper, log always denote the natural logarithm. We use Ps to represent
the Poisson distribution with parameter s (for any real number s > 0). This discrete probability
−s s
distribution is defined as Ps (k) = e k!k for any integer k ≥ 0.
For any real number a < b, the notation ( a, b), respectively ( a, b], represents the set of real
numbers x such that a < x < b, respectively a < x ≤ b. Given two real numbers a, b and a
positive real number δ, we say that a is a δ-additive approximation of b if | a − b| ≤ δ, and say
that a is a δ-multiplicative approximation of b if | a − b| ≤ δ| a − b|.
We will generally work in the vector space R n for some positive integer n (representing the
number of vertices of the graph considered). Given a vector v ∈ R n , we write kvk its Euclidean
norm, and vt its transpose. Given two vectors u, v ∈ R n , we write hu, vi their inner product.
Given any s ∈ {1, . . . , n}, we use es to denote the vector in R n with s-th coordinate 1 and all other
coordinates zero.
The Laplacian and its eigenvalues. We now define the Laplacian of a graph and introduce
more specific notations used through this paper. It will be more convenient for us to work with
the normalized version of the Laplacian of a graph, due to its natural connections with random
walks. We refer to [3] for details on the normalized Laplacian and all the notions below.
3

Let G = (V, E) be an undirected weighted graph (possibly with loops) on n vertices, i.e., a
graph with a weight function w : V × V → R satisfying w(i, j) = w( j, i ) for all pairs of vertices
(i, j) ∈ V × V, w(i, j) = 0 if (i, j) ∈
/ E and w(i, j) > 0 if (i, j) ∈ E. The degree of a vertex i ∈ V,
denoted di , is defined as di = ∑ j∈V w(i, j). The degree matrix of the graph G is the n × n diagonal
matrix DG = diag(d1 , . . . , dn ). The volume of G is defined as vol( G ) = ∑nℓ=1 dℓ . Finally, if the
graph has no isolated vertex, i.e., the degree of each vertex is positive, let us define the quantity
d=

maxi∈{1,...,n} di
.
mini∈{1,...,n} di

The (normalized) Laplacian of the graph G is

w ( i,i)


 1 − di
w ( i,j)
L G [i, j] =
−√
di d j



0

(1)

the n × n symmetric matrix L G such that
if i = j and di 6= 0,

if (i, j) ∈ E,
otherwise,

for all i, j ∈ {1, . . . , n}. Let λ1 ≤ λ2 · · · ≤ λn denote its eigenvalues. It is known that λ1 ≥ 0
and λn ≤ 2. Moreover, λ2 6= 0 if and only if G is a connected graph. We have the general upper
bound λ2 ≤ n/(n − 1). When G is connected, we obtain the lower bound
λ2 ≥

1
,
diam( G )vol( G )

(2)

where diam( G ) denotes the diameter of G (the maximal distance between two vertices in G). Let
{u1 , . . . , un } be an orthonormal basis for L G , where ui is an eigenvector associated with eigen√
√ t
value λi . We can thus write L G = ∑ni=1 λi ui uti . It is easy to check that u1 = √ 1
d1 , . . . , d n .
vol( G )

Image of L G and its pseudo inverse. We define Im( L G ) as the image of L G , that is, the linear
span of the vectors uc , uc+1 , . . . , un , where c is the smallest integer such that λc 6= 0. The pseudo
inverse (also called the Moore-Penrose pseudo inverse) of L G , denoted L†G , is the matrix
L†G =

n

∑ λ−i 1 ui uti ,

i= c

which corresponds to the inverse of L G when restricted to the subspace Im( L G ). Note that if the
graph G is connected then c = 2. In that case Im( L G ) is simply the subspace of R n consisting of
all vectors orthogonal to u1 .
The transition matrix and random walks. Another useful matrix is the transition matrix of G,
which represents one step of a random walk on the graph G. Assume that the graph has no
isolated vertex. Then the transition matrix of G is the n × n matrix PG defined as

w(i, j)/di if (i, j) ∈ E,
PG [i, j] =
0
otherwise
−1/2
for all i, j ∈ {1, . . . , n}. Note that L G = I − D1/2
.
G PG D G
Laplacian systems. Given a Laplacian L G of an undirected weighted graph G and a vector
b ∈ Im( L G ), the corresponding Laplacian system is the equation L G x = b. The goal is to solve
this equation, i.e., compute L†G b or to compute an (additive or approximative) approximation of

4

it. A standard assumption is that G is connected since one can always reduce to this case dealing
individually with each connected component (note that this reduction is space-efficient since one
can identify the connected components of a graph in logarithmic space). For this reason, in
Sections 3 and 4 of this paper we will assume that G is connected. For convenience, and without
loss of generality, we will also assume that b is a unit-norm vector.
Remark on numerical precision. The algorithms we present in this paper perform arithmetic
operations on real numbers. When working with real numbers numerical precision is always a
delicate issue, and especially in a space-bounded setting. To illustrate this point, let√
us consider
as a simple example the multiplication of an integer a by the irrational numbers 2. One of
the most rigorous ways to√proceed is to work in the bit complexity model, and compute only a
fixed number of bits of a 2, which necessarily introduce a (small) approximation error. Analyzing rigorously a complicated algorithm in this model is typically extremely tedious, since the
approximation errors of essentially all the arithmetic steps of the algorithm have to been considered. In this paper we (implicitly) use a slightly more abstract, but still fairly standard, model
where basic arithmetic operations involving “reasonable” numbers can be implemented exactly.
More precisely, we assume that the standard arithmetic operations (addition, substraction, multiplication, division) on O(log m)-bit reals numbers (i.e., real numbers of absolute value between
′
m−c and mc for some constant c and c′ ) can be done exactly in O(log m) space. This model
enables us to focus on the most interesting and important algorithmic aspects of our approach,
without having to deal with a multitude of minor technical details. Note that all the statements
of the technical results of Sections 3 and 4 hold with only minor modifications in the bit complexity model as well. The statements of Theorems 1.1 and 1.2 in the introduction, described in
the more convenient setting of graphs with polynomially bounded weighted, hold without any
modification even in the bit complexity model.

3 Space-efficient Laplacian Solver
In this section we describe how to compute in logarithmic space an approximate solution of the
Laplacian system L G x = b, and prove Theorem 1.1. We assume that G is connected.
The theoretical foundation of our algorithm comes from the formula of the following theorem,
which shows that the solution L†G b to the Laplacian system can be approximated using powers
of the transition matrix PG of the graph.
Theorem 3.1. Let ǫ be any positive real number and λ ∈ (0, 2] be a lower bound on all the non-zero
log (6/( ǫλ))
, N ≥ 6T
eigenvalues of L G (i.e., a lower bound on λ2 ). For any positive integers T ≥
λ
ǫ and
K ≥ max(6T, log(6T/ǫ)), the inequality


L†G −

T
N

N K −1


ǫ k bk
1/2 k −1/2
P
(
k
)
D
P
D
b ≤
∑ ∑ jT/N
G G
G
2
j=1 k=0

holds for any vector b ∈ Im( L G ).

k −1/2 =
Proof. Observe that D1/2
∑ni=1 (1 − λi )k uti ui for any integer k ≥ 0. We show below that
G PG D G
the inequality
T N K −1
1
− ∑ ∑ P jT/N (k)(1 − λi )k ≤ ǫ/2,
λi
N j=1 k=0

5

holds for all i ∈ {2, . . . , n}, which will imply the statement
R ∞ −tλ in the theorem.
1
For any i ∈ {2, . . . , n}, first observe that λi = 0 e i dt. For any T ≥ 0, we thus have
1
−
λi

Z T
0

e−tλi dt ≤

Z ∞
T

e−tλi dt =

e− Tλi
.
λi

(3)

Let us approximate the integral by a right Riemann sum: for any integer N ≥ 1,
Z T
0

e−tλi dt −

T
N

N

∑ e− j

T
N λi

j=1

≤

T
T
× (1 − e− Tλi ) ≤ .
N
N

(4)

For any s ≥ 0 we have
∞

e−sλi = e−s e(1−λi )s = e−s

∞
((1 − λi )s)k
=
∑
∑ Ps (k)(1 − λi )k .
k!
k=0
k=0

For any integer K ≥ 1 we thus have
e

− sλi

−

K −1

∞

∑ Ps (k)(1 − λi )

k

k=0

≤

∑ P s ( k) ≤

k= K



s
1−
K+1

 −1

Ps (K ) ≤



s
1−
K+1

 −1

√

sK e− s
2πKK K e−K

The first inequality uses λi ∈ (0, 2]. The second inequality uses a well-known upper bound on
the tail probability of the Poisson distribution [7]. The third inequality uses a standard lower
bound on K! from Stirling approximation. For any positive integer K ≥ e2 s, where e = 2.718 . . .
is Euler’s number, and in particular for any K ≥ 6s, we get
e

− sλi

−

K −1

∑

k=0

Ps (k)(1 − λi )

k

≤



1
1− 2
e

 −1

√

e− s
2πKeK

≤

1
.
eK

(5)

Combining Inequalities (3), (4) and (5), for any positive integer K ≥ 6T we obtain:
1
T
−
λi
N

N K −1

∑∑

j=1 k=0

P jT/N (k)(1 − λi )k ≤

e− Tλi
T
T
+ + K.
λi
N e

Taking values of T, N and K as in the statement of the theorem guarantees that the right side of
the above inequality is at most ǫ/2.
We now present a lemma that shows how to approximate space-efficiently the Poisson distribution appearing in the formula of Theorem 3.1. This is done by using the well-known property
that the Poisson distribution can be expressed as the limit distribution of binomial random variables, and showing that the convergence is fast enough.
Lemma 3.2. For any δ ∈ (0, 1) and any ζ > 0, there exists a O(log (ks/δ) + log log(1/ζ ))-space
algorithm that outputs a δ-additive approximation of Ps (k) with probability at least 1 − ζ.

6

.

Proof. Let n be any integer such that n ≥ 2(k2 + s2 )/δ. Let Xn be the binomial random variable
with parameters n and s/n, i.e., Xn is the number of successes in n repeated trials of a binomial
experiment with success probability s/n. Standard computations (see, e.g., page 99 of [10]) show
that




1
s2
k
Ps (k)(1 − k/n) 1 −
≤ Pr[Xn = k] ≤ Ps (k)
.
n
1 − sk/n
Using the inequality 1/(1 − x) ≤ 1 + 2x valid for any x ∈ (0, 1/2) and the inequality (1 − x)k ≥
1 − kx valid for any k ≥ 0, we get
k2 + s2
P s ( k) 1 −
n






2sk
≤ Pr[Xn = k] ≤ Ps (k) 1 +
n



.

With our value of n, and since Ps (k) ≤ 1, we obtain

|Ps (k) − Pr[Xn = k]| ≤ δ/2.

(6)


Our algorithm for estimating Ps (k) is as follows. We will use the value
n = m 2(k2 + s2 )/δ .
l
2 log(2/ζ )
times the folThe algorithm creates a counter C initialized to zero, and repeat m =
δ2
lowing: perform n trials of a binomial experiment with success probability s/n and increment C
by one if exactly k among these n trials succeeded. The algorithm finally outputs C/m.
The expected value of the output of this algorithm is precisely the quantity Pr[ Xn = k] considered above. From Chernoff bound, the output of the algorithm is thus a δ/2-additive ap2
proximation of this quantity with probability at least 1 − 2e−mδ /2 ≥ 1 − ζ. Combining this with
Inequality (6) and the triangular inequality, we conclude that the output of the algorithm is a
δ-additive approximation of Ps (k) with the same probability.
Finally, observe that the space complexity of the algorithm is linear in log(1/δ), log k, log s
and log log(1/ζ ).
We are now ready to present our algorithm for estimating L†G b. The algorithm, denoted
Algorithm A, is described in Figure 1 and analyzed in the following theorem.
Theorem 3.3. Let G be an undirected weighted connected graph and λ ∈ (0, 2] be a lower bound on the
second smallest zero eigenvalue λ2 of G. Algorithm A outputsan ǫ-additive approximation of the i-th
nd
entry of L†G b with probability at least 1 − γ, and uses O(log ǫλ
+ log log(1/γ)) space.

Proof. Remember that for any s ∈ {1, . . . , n}, es denotes the 1 × n vector with s-th coordinate 1
and all other coordinates zero. Note that the i-th coordinate of L†G b, which is the quantity we want
to estimate, is ei , bt L†G since L†G is symmetric. From Theorem 3.1 and the triangular inequality,
we know that any 2ǫ -additive approximation of
T
N

N K −1

∑∑

j=1 k=0

D

k −1/2
P jT/N (k) ei , bt D1/2
G PG D G

E

T
=
N

N K −1 n

∑∑

∑ P jT/N (k)bℓ

j=1 k=0 ℓ=1

s

E
dℓ D
ei , eℓ PGk
di

is an ǫ-additive approximation of ei , bt L†G . We show below that Algorithm A precisely outputs
an ǫ/2-additive approximation of this quantity.

7

Input: an undirected, weighted and connected graph G on n vertices,
an integer i ∈ {1, . . . , n}, a unit vector b ∈ Im( L G )
a precision parameter ǫ ∈ (0, 1], an error parameter γ ∈ (0, 1]
1
2

T←

l



log (6/( ǫλ))
λ

δ ← 6TK

q

m

;N←
 −1

∑nℓ=1

dℓ
di



; K ← max(6T, log( 6T
ǫ )) ;
m
l
log (2/ζ )
; ζ ← NK(γ1+n) ; r ←
;
2
2δ
 6T 
ǫ

3 R ← 0;
4 for j from 1 to N
5
for k from 0 to K − 1
6
Compute an approximation a of P jT/N (k) using Lemma 3.2 with δ and ζ;
7
for ℓ from 1 to n
8
S ← 0;
9
repeat r times:
10
Run the walk PG starting on vertex ℓ for k steps;
11
If the walk q
ends on vertex i then S ← S + 1;
abℓ S
12
R ← R + r × ddℓi ;
13 Output RT/N.
Figure 1: Algorithm A computing an ǫ-additive approximation of the i-th entry of L†G b with
probability at least 1 − γ.
Note that the probability of a walk PG starting on vertex ℓ reaches vertex i after exactly k steps
is hei , eℓ PGk i, the i-th coordinate of eℓ PGk . At the end of Steps 9-11 we thus have
h
i
2
Pr |S/r − hei , eℓ PGk i| ≤ δ ≥ 1 − 2e−2rδ ≥ 1 − ζ,
(7)

from Chernoff bound. Lemma 3.2 also shows that a is a δ-additive approximation of P jT/N (k)
with probability at least 1 − ζ. Let us continue our analysis under the assumption that all these
approximations are correct (we discuss the overall success probability at the end of the proof).
At Step 12 we thus have
E
E
D

D
aS
δ ≤ 3δ.
− P jT/N (k) ei , eℓ PGk ≤ δ + P jT/N (k) + ei , eℓ PGk
r

The output of the algorithm at Step 13 then satisfies
s
s
s
E
n
N K −1 n
RT
T
T N K −1 n
dℓ D
d
dℓ
ℓ
3δ|bℓ |
ei , eℓ PGk ≤
= TK ∑ 3δ|bℓ |
− ∑ ∑ ∑ P jT/N (k)bℓ
∑
∑
∑
N
N j=1 k=0 ℓ=1
di
N j=1 k=0 ℓ=1
di
di
ℓ=1
s
n
d
≤ 3δTK ∑ ℓ ≤ ǫ/2.
d
ℓ=1 i
The space complexity of the algorithm is O(log(nd/(λǫ)) + log log(1/γ)), from Lemma 3.2
and the observation that only registers of this size are needed to implement the algorithm. Finally,
8

let us discuss the success probability of this algorithm. Errors can only occur at Steps 6 or 10.
From Inequality (7) and Lemma 3.2, and using the union bound, we know that the overall success
probability is at least 1 − ζNK (1 + n) = 1 − γ.
Theorem 3.3 implies Theorem 1.1 by observing that d/λ can be upper bounded by a polynomial in n when the weights are polynomially bounded, as shown in Equations (1) and (2) of
Section 2 (if G is not connected we can simply apply Theorem 4.4 on each connected component).

4 Space-efficient Approximation of the Spectral Gap
Let us consider the following matrix:
MG =


1
−1/2
I + D1/2
P
D
.
G G
G
2

Note that MG is a symmetric matrix. Its eigenvalues are 0 ≤ 1 − λn /2 < · · · < 1 − λ2 /2 <
1 − λ1 /2 = 1. The eigenvectors of MG are the same as the eigenvectors of L G . In particular, the
eigenvector of MG corresponding to the eigenvalue 1 is u1 , and the eigenvector of MG corresponding to the eigenvalue 1 − λ2 /2 is u2 . A well-known approach for approximating the largest
eigenvalue of a matrix is the power method (see, e.g., [17]). Our idea is to apply this method on
k+1
k vk on a random vector v ∈ Im( L )
MG restricted to Im( L G ), and compute the ratio k MG
vk/k M G
G
— it is easy to show that with high probability this ratio is close to 1 − λ2 for large enough k. In
this section we will develop a space-efficient version of this approach, and prove Theorem 1.2
When using the power method to estimate 1 − λ2 /2, we need to use vectors orthogonal to
u1 that have a non-negliglible ”overlap” with the eigenvector u2 . We say that a vector v ∈ R n
is good if the following three conditions are satisfied: kvk = 1, v ∈ Im( L G ), and |hv, u2 i| ≥
√ 1 . While a random unit-norm vector in Im( L G ) is a good vector with high probability, several
2nd
technical difficulties arise when considering space-efficient vector sampling. Instead of using
such probabilistic arguments, we introduce below a set Σ of vectors that necessarily contains at
least one good vector.
Let Σ ⊂ R n be the set containing the n(n − 1)/2 vectors defined as follows. Each of these vectors corresponds to choosing two distinct indexes i, j ∈ {1, . . . , n} and taking the n-dimensional
vector with i-th coordinate − √ 1
, j-th coordinate √ 1
, and all other coordinates being
1+ di /d j

1+ d j /di

zero. The following easy lemma shows that Σ indeed contains at least one good vector.
Lemma 4.1. There exists a good vector in Σ.

Proof. Any vector in Σ is a unit-norm vector orthogonal to u1 . We show below that there exists a
vector in Σ that also satisfies the third condition of the definition of good vectors.
Let us write u2 = ( x1 , . . . , xn ). Let S+ ⊆ {1, . . . , n} be the set of indices ℓ such that xℓ ≥ 0,
and S−√⊆ {1, . . . , n} be √
the set of indices ℓ such that xℓ < 0. Since u2 is orthogonal to u1 , we have
∑ℓ∈S+ dℓ | xℓ | = ∑ℓ∈S− dℓ | xℓ |. Since u2 is a unit vector we have
1=

∑+ x2ℓ + ∑− x2ℓ ≤ ∑+ |xℓ | + ∑− |xℓ |.

ℓ∈S

ℓ∈S

ℓ∈S

9

ℓ∈S

We conclude that

∑+

ℓ∈S

p

√

minℓ∈{1,...,n}
dℓ | xℓ | ≥
2

dℓ

and

∑−

ℓ∈S

p

minℓ∈{1,...,n}
dℓ | xℓ | ≥
2

√

dℓ

which implies that there exist i ∈ S− and j ∈ S+ such that | xi | and | x j | are at least 2n1√d . Let v be
and j-th coordinate √ 1
. The inner product
the vector in Σ with i-th coordinate − √ 1
1+ di /d j

1+ d j /di

of v and u2 is thus at least

1
1
p
+p
1 + di /d j
1 + d j /di

1
√
2n d
as claimed.

!

dj

d
p
+p i
d j + di
di + d j

1
= √
2n d

!

≥√

1
2nd

,

The following proposition is our version of the power method.
Proposition 4.1. Let δ be any real number such that 0 < δ ≤ 1, and ζ be any real number such that
0 < ζ ≤ δλ2 /12.
(i) For any integer k ≥ 0, any non-zero vector v ∈ Im( L G ) and any ζ-multiplicative approximations
k vk and k M k+1 vk, respectively, the inequality (1 − δ ) λ ≤ 2 1 − C /C holds.
C1 and C2 of k MG
(
2
2
1)
G
C1 and C2 of

k MGk vk

and

k MGk+1 vk,

√
3 log( 2nd)
δλ2

− 1 and any ζ-multiplicative approximations
respectively, the inequality 2 (1 − C2 /C1 ) ≤ (1 + δ)λ2 holds

(ii) Let v be a good vector. For any integer k ≥
Proof. Let us first prove part (i). We have
C2
≤
C1

k+1
MG
v
kv
MG

1+ζ
1+ζ
×
≤ (1 − λ2 /2) ×
=
1−ζ
1−ζ



2ζ
1+
1−ζ



(1 − λ2 /2)

and thus






2ζλ2
4ζ
C2
4ζ
5ζ
+
≥ 1−
λ2 ≥ (1 − δ ) λ2 ,
2 1−
≥ λ2 −
λ2 ≥ 1 −
C1
1−ζ
1−ζ
λ2 (1 − ζ )
λ2
where the third inequality was obtained from ζ ≤ 1/5 (from the assumption ζ ≤ δλ2 /12).
Let us now prove part (ii). For any unit vector v, Hölder’s inequality gives
k
v
MG

2

≤

∑
i

+2
v2i λ2k
i

! k/(k+1)

∑
i

v2i

!1/(k+1)

= MGk+1 v

2k/( k+1)

If v is good then

k MGk+1 vk
k MGk+1 vk
(1 − λ2 /2)
≥
= k MGk+1 vk1/(k+1) ≥ √
.
k
k
+
1
k M G vk
k MG vkk/(k+1)
( 2nd)1/(k+1)

10

.

Taking k + 1 ≥

√
3 log ( 2nd)
δλ2

gives

√ 1
( 2nd)1/(k +1)

statement of the proposition. We have

≥ e−δλ2 /3 ≥ 1 − δλ2 /3. Let C1 and C2 be as in the

k MGk+1 vk 1 − ζ
C2
×
≥
≥ (1 − λ2 /2)(1 − δλ2 /3) × (1 − 2ζ )
C1
1+ζ
k MGk vk

δλ2
λ2
δλ2 2ζδλ2
+
−
+ λ2 ζ + 2
≥ 1 − 2ζ −
3
3
2
3
δλ2 λ2
≥ 1 − 2ζ −
− ,
3
2



1
−ζ
2



where the last inequality uses ζ ≤ 1/2 (which is guaranteed from the assumption ζ ≤ δλ2 /12).
We thus obtain 2 (1 − C2 /C1 ) ≤ (1 + 2δ/3 + 4ζ/λ2 ) λ2 ≤ (1 + δ)λ2 , as claimed.
k vk and k M k+1 vk to apProposition 4.1 requires good multiplicative approximations of k MG
G
proximate λ2 . Using random walks, we nevertheless will only be able to obtain additive approximations. To convert additive approximations into a good multiplicative approximations, we will
need lower bounds on these two quantities. We will also need upper bounds in order to control
the running time (and the space complexity) of our algorithm. The following lemma shows the
bounds we will use.
k vk < τ
Lemma 4.2. Let τ be any real number such that τ ∈ (0, 1]. For any unit vector v ∈ Im( L G ), k MG

for all integers k >

2 log (1/τ )
.
λ2

k+1
Additionally, if v is good and n ≥ 4 then k MG
vk ≥ 2τ for all integers

√
log(1/τ ) − log(2 2nd)
− 1.
k≤
2λ2

(8)

Proof. For any unit vector v ∈ Im( L G ) we have
k
v ≤ (1 − λ2 /2)k ≤ e−kλ2 /2 ,
MG

which is upper bounded by τ for k >
k+1
v
MG

2 log (1/τ )
.
λ2

If v is good we further have

(1 − λ2 /2)k+1
1
√
≥
≥√
2nd
2nd



e

−1

λ2
−
4

(k+1)λ2 /2

,

1
where we used the formula e−1 − 2a
≤ (1 − 1/a) a valid for any a ≥ 1 (see, e.g., [9]) with a = 2/λ2 .
Note that
3
n
>
e−1 − λ2 /4 ≥ e−1 −
4( n − 1)
100
k+1
v ≥ 2τ whenever
for n ≥ 4. We get MG

√
√
2
log(1/τ ) − log(2 2nd)
log(2 2τnd)
2
=
×
×
.
k+1 ≤
λ2
log(3/100)
log(100/3)
λ2
Finally, note that 2/ log(100/3) > 1/2.
11

k v spaceIn order to use the the theory developed above, we need to be able to estimate MG
efficiently. This can be done using an approach based on quantum walks, similarly to what we
did in Section 3. The description the procedure based on this idea, and its analysis, are given
in the appendix. This procedure, denoted Estimate-norm( G, k, v, ǫ, γ), computes an ǫ-additive
k v with probability at least 1 − γ, for any k ≥ 0, any v ∈ Σ and any
approximation of MG
ǫ, γ ∈ (0, 1]. We state the main result of the appendix as the following theorem.

Theorem 4.3. Procedure Estimate-norm( G, k, v, ǫ, γ) uses O(log(nkd/ǫ) + log log(1/γ)) space and
k v with probability at least 1 − γ.
outputs an ǫ-additive approximation of MG
Our algorithm for estimating λ2 , denoted Algorithm B , is given in Figure 2 and analyzed in
Theorem 4.4 below.

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

Input: an undirected, weighted and connected graph G on n vertices, where n ≥ 4,
a precision parameter δ ∈ (0, 1], an error parameter γ ∈ (0, 1]


log (1/τ ) −1
4γ
×
1
+
;
;
ζ
←
τ ← √ 1 1+8/δ ; ǫ ← δλτ
12
λ2
n ( n − 1)
2( 2nd)

Rmax ← 0;
for all v ∈ Σ
k ← 1;
C1 ← Estimate-norm( G, 1, v, ǫ, ζ );
C2 ← Estimate-norm( G, 2, v, ǫ, ζ );

kv
# C1 will store an approximation of MG
k+1
# C2 will store an approximation of MG
v

repeat until C2 < 3τ/2
if C2 /C1 > Rmax
Rmax ← C2 /C1 ;
k ← k + 1;
C1 ← C2 ;
C2 ← Estimate-norm( G, k + 1, v, ǫ, ζ );
Output 2(1 − Rmax ).

Figure 2: Algorithm B computing an δ-multiplicative approximation of λ2 with probability at
least 1 − γ.
Theorem 4.4. Let G be an undirected weighted connected graph and λ ∈ (0, 2] be a lower bound on λ2 .
Algorithm B outputs a δ-multiplicative approximation of λ2 with probability at least 1 − γ, and uses
O( 1δ log(nd/λ) + log log(1/γ)) space.
Proof. Let us first analyze Algorithm B under the assumption that at Steps 5, 6 and 12, Procedure
k v k.
Estimate-norm( G, k, v, ǫ, ζ ) always correctly outputs an ǫ-additive approximation of k MG
k vk
Then during the execution of the algorithm, C1 and C2 are ǫ-additive approximations of k MG
k+1
and k MG
vk, respectively.
k vk < τ holds for all integers
Lemma 4.2 guarantees that for each v ∈ Σ, the inequality k MG
k > 2 log(1/τ )/λ2 , in which case we have C1 < τ + ǫ < 3τ/2 (the same inequality holds for C2 )
since ǫ < τ/2. For each v ∈ Σ, the loop of Steps 7-12 is thus repeated at most 2 log(1/τ )/λ2
times.
Whenever the ratio C2 /C1 is computed at Step 8-9, we have C1 > 3τ/2 and C2 > 3τ/2. Since
k vk > τ and k M k+1 vk > τ. In this case the quantities C and C are
ǫ < τ/2, this means that k MG
2
1
G
12

k vk and k M k+1 vk, respectively. From part (i) of
thus also τǫ -multiplicative approximations of k MG
G
Proposition 4.1 with ζ = ǫ/τ, we conclude that 2(1 − Rmax ) ≥ (1 − δ)λ2 .
Let v be a good vector. Observe that with the choice of τ made at Step 1 we have
√
√
√
log(1/τ ) − log(2 2nd)
4 log( 2nd)
3 log( 2nd)
=
≥
+1
2λ2
δλ2
δλ2

for n ≥ 4. There thus necessarily exists at least one integer k satisfying both the condition of
k+1
Part (ii) of Proposition 4.1 and Inequality (8). From Lemma 4.2, for such a k we have k MG
vk ≥
2τ, which gives the lower bound C2 ≥ 2τ − ǫ ≥ 3τ/2 (the same inequality holds for C1 ), and
implies that the ratio C2 /C1 is computed at Steps 8-9. Part (ii) of Proposition 4.1 and Lemma 4.1
thus imply 2(1 − Rmax ) ≤ (1 + δ)λ2 , as claimed.
From Theorem 4.3, each call of Procedure Estimate-norm requires at most
1
O(log(nkd/ǫ) + log log(1/ζ )) = O( log(nd/λ) + log log(1/γ))
δ
space. This bound is also an upper bound on the space complexity of all the other computational
steps of Algorithm B . Let us conclude by discussing the success probability of this algorithm.
Each application of Estimate-norm errs with probability at most ζ, from Theorem 4.3. There
are at most 2 + 2 log(1/τ )/λ2 calls to this procedure. The success probability is thus at least
1 − |Σ| (2 + 2 log(1/τ )/λ) ζ ≥ 1 − γ.
Theorem 4.4 implies Theorem 1.2 by observing again that d/λ can be upper bounded by a
polynomial in n when the weights are polynomially bounded (if G is not connected then we
apply Theorem 4.4 on each connected component and taking the minimum of the estimations
obtained).

Acknowkedgments
The author is grateful to Richard Cleve, Hirotada Kobayashi, Harumichi Nishimura, Suguru
Tamaki and Ryan Williams for helpful comments. This work is supported by the Grant-in-Aid for
Young Scientists (A) No. 16H05853, the Grant-in-Aid for Scientific Research (A) No. 16H01705,
and the Grant-in-Aid for Scientific Research on Innovative Areas No. 24106009 of the Japan
Society for the Promotion of Science and the Ministry of Education, Culture, Sports, Science and
Technology in Japan.

References
[1] Stuart J. Berkovitz. On computing the determinant in small parallel time using a small
number of processors. Information Processing Letters, pages 147–150, 1984.
[2] Allan Borodin, Joachim von zur Gathen, and John E. Hopcroft. Fast parallel matrix and
GCD computations. Information and Control, 52(3):241–256, 1982.
[3] Fan R. K. Chung. Spectral Graph Theory. American Mathematical Society, 1997.

13

[4] Fan R. K. Chung and Olivia Simpson. Solving local linear systems with boundary conditions
using heat kernel Pagerank. Internet Mathematics, 11:4–5, 2015.
[5] Laszlo Csanky. Fast parallel matrix inversion algorithms. In Proceedings of the 16th Annual
Symposium on Foundations of Computer Science, pages 11–12, 1975.
[6] Dean Doron and Amnon Ta-Shma. On the problem of approximating the eigenvalues of
undirected graphs in probabilistic logspace. In Proceedings of the 42nd International Colloquium, pages 419–431, 2015.
[7] Peter W. Glynn. Upper bounds on Poisson tail probabilities. Operation Research Letters,
6(1):9–14, 1987.
[8] Jonathan A. Kelner, Lorenzo Orecchia, Aaron Sidford, and Zeyuan Allen Zhu. A simple,
combinatorial algorithm for solving SDD systems in nearly-linear time. In Proceedings of the
45th Symposium on Theory of Computing, pages 911–920, 2013.
[9] Dragoslav Mitrinović. Analytic Inequalities. Springer, 1970.
[10] Michael Mitzenmacher and Eli Upfal. Probability and Computing. Cambridge University
Press, 2005.
[11] Richard Peng and Daniel A. Spielman. An efficient parallel solver for SDD linear systems.
In Proceedings of the 46th Symposium on Theory of Computing, pages 333–342, 2014.
[12] Daniel A. Spielman and Shang-Hua Teng. Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems. In Proceedings of the 36th Annual
Symposium on Theory of Computing, pages 81–90, 2004.
[13] Daniel A. Spielman and Shang-Hua Teng. Spectral sparsification of graphs. SIAM Journal
on Computing, 40(4):981–1025, 2011.
[14] Daniel A. Spielman and Shang-Hua Teng. A local clustering algorithm for massive graphs
and its application to nearly linear time graph partitioning. SIAM Journal on Computing,
42(1):1–26, 2013.
[15] Daniel A. Spielman and Shang-Hua 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.
[16] Amnon Ta-Shma. Inverting well conditioned matrices in quantum logspace. In Proceedings
of the 45th Symposium on Theory of Computing, pages 881–890, 2013.
[17] Nisheeth K. Vishnoi. Lx = b — Laplacian Solvers and their Algorithmic Applications. Now
publishers, 2013.

14

Appendix: Estimating MGk v
k v for any vector v ∈ Σ, and
In this appendix we explain how to space-efficiently estimate MG
prove Theorem 4.3.
The following lemma first shows how to space-efficiently estimate the quantity 2−k (ks).

Lemma 4.5. Let δ and ζ be any real numbers such that δ ∈ (0, 1) and ζ > 0. There exists a
O(log (ks/δ) + log log(1/ζ ))-space algorithm that, when given as input two integers k ≥ 1 and s ∈
{0, . . . k}, outputs an δ-additive approximation of 2−k (ks) with probability 1 − ζ.
Proof. Our algorithm
m is as follows. The algorithm creates a counter C initialized to zero, and
l
log(2/ζ )
times the following: take s bits uniformly at random and increment C by
repeat m =
2δ2
one if exactly k among these s bits are one. The algorithm finally outputs C/m. Observe that this
algorithm can be implemented in space linear in log(1/δ), log k, log s and log log(1/ζ ).
The expected value of the output of this algorithm is precisely 2−k (ks). From Chernoff bound,
the output of the algorithm is thus a δ-additive approximation of this quantity with probability
2
at least 1 − 2e−2mδ ≥ 1 − ζ.
k v is described in Figure 3.
The procedure Estimate-norm( G, k, v, ǫ, γ) estimating MG
We now prove Theorem 4.3 stated in Section 4.

Proof of Theorem 4.3. Observe that
k
MG

n

=

∑ (1 −

λi /2)k uti ui

=

i=1

Moreover

k
MG

is symmetric. Let us write v =

procedure. We thus have
k
v
MG

where

(k)
Γi = ∑ sk
s =0 2
k

s

= v

t

D1/2
G

D1/2
G
r

(ks) s
∑ k PG
s =0 2
k

(ks) s
∑ k PG
s =0 2
k

dℓ2
d ℓ1 + d ℓ2 e ℓ1

!

dℓ1 dℓ2
he , e Ps i −
( dℓ1 + dℓ2 ) di i ℓ1 G

−1/2
DG

s

−

!

r

=

−1/2
DG
.

dℓ1
d ℓ1 + d ℓ2 e ℓ2

s

as in Step 2 of the

n

∑ Γ2i

i=1

dℓ1 dℓ2
he , e Ps i
( dℓ1 + dℓ2 ) di i ℓ2 G

!

.

Note that |Γi | ≤ 1 for any i ∈ {1, . . . , n}.
At the end of Steps 9-11 we have
2

Pr [|S/r − hei , eℓ1 PGs i| ≤ δ] ≥ 1 − 2e−2rδ ≥ 1 − ζ,

(9)

using the same argument as in the analysis in Theorem 3.3. The same bound, with ℓ1 replaced
by ℓ2 , holds at the end of Steps 14-16 as well. Lemma 4.5 also shows that a is a δ-additive
approximation of 2−k (ks) with probability at least 1 − ζ at Step 7. Let us continue our analysis

15

Input: an undirected, weighted and connected graph G on n vertices,
a vector v ∈ Σ, a positive integer k,
a precision parameter ǫ ∈ (0, 1], an error parameter γ ∈ (0, 1]
1

δ←

√
ǫ2 2
;
54( k+1) nd

ζ←

γ
;
3n ( k+1)

r←

l

m

log (2/ζ )
;
2δ2
r

2 Let (ℓ1 , ℓ2 ) be the indices such that v =

dℓ2
d ℓ1 + d ℓ2 e ℓ1

−

r

dℓ1
d ℓ1 + d ℓ2 e ℓ2 ;

3 R ← 0;
4 for i from 1 to n
5
Q ← 0;
6
for s from 0 to k
7
Compute an estimate a of 2−k (ks) using Lemma 4.5 with δ and ζ;
8
S ← 0;
9
repeat r times:
10
Run the walk PG starting on vertex ℓ1 for s steps;
11
If the walk
r ends on vertex i then S ← S + 1;
dℓ1 dℓ2
;
12
Q ← Q + aS
r
(d +d )d
ℓ1

13
14
15
16
17
18
19

ℓ2

i

S ← 0;
repeat r times:
Run the walk PG starting on vertex ℓ2 for s steps;
If the walk
r ends on vertex i then S ← S + 1;
Q ← Q−

aS
r

dℓ1 dℓ2
( dℓ1 + dℓ2 ) di ;

2
R←
√R + Q ;
Output R.

Figure 3: Procedure Estimate-norm( G, k, v, ǫ, γ) computing an ǫ-additive approximation of
k v with probability at least 1 − γ.
MG
under the assumption that all these approximations are correct (we discuss the overall success
probability at the end of the proof). At Step 12 we thus have
!
k
k
(
)
aS ( s )
− k hei , eℓ1 PGs i ≤ δ + sk + hei , eℓ1 PGs i δ ≤ 3δ,
r
2
2
and the same bound, with ℓ1 replaced by ℓ2 , holds at Step 17 as well. For any i ∈ {1, . . . , n} we
thus have
s
dℓ1 dℓ2
6
≤ √ (k + 1)δd.
|Q − Γi | ≤ 6(k + 1)δ max
( dℓ1 + dℓ2 ) di
ℓ1 ,ℓ2 ∈{1,...,n}
2
at the end of the loop of Steps 5-17, which implies
18
Q2 − Γ2i = | Q − Γi | × | Q + Γi | ≤ √ (k + 1)δd.
2
16

since |Γi | ≤ 1 and | Q| ≤ | Q − Γi | + |Γi | ≤ 2. We thus have
18
≤ √ (k + 1)nδd ≤ ǫ2 /3.
2
√
k v by
at the end of of the algorithm. Let us show that R is an ǫ-additive estimation of MG
p
considering two cases. In the case | R| ≤ ǫ/3 we get
r
p
√
√
√
k
k
k v 2 ≤ ǫ/3 +
R − MG v ≤ R + MG v ≤ R + R + R − MG
ǫ2 /9 + ǫ2 /3 = ǫ.
k
v
R − MG

Now in the case

p

2

| R| ≥ ǫ/3 we get
√

2

k
v
R − MG

kv
R − MG
= √
kv
R + MG

≤

ǫ2 /3
= ǫ.
ǫ/3

The space complexity of this algorithm is O(log(nk/ǫ) + log log γ), from Lemma 4.5 and the
observation that all other computational steps can be implemented with registers of this size.
Let us conclude by discussing the success probability of this algorithm. Errors can only occur at
Steps 6 or 10. From Inequality (9) and Lemma 4.5, and using the union bound, we know that the
overall success probability is at least 1 − 3ζn(k + 1) = 1 − γ.

17

