Similarity-Aware Spectral Sparsification by Edge Filtering
Zhuo Feng

arXiv:1711.05135v2 [cs.DS] 15 Nov 2017

Department of Electrical and Computer Engineering
Michigan Technological University
Email: zhuofeng@mtu.edu

Abstract—In recent years, spectral graph sparsification techniques that
can compute ultra-sparse graph proxies have been extensively studied
for accelerating various numerical and graph-related applications. Prior
nearly-linear-time spectral sparsification methods first extract low-stretch
spanning tree of the original graph to form the backbone of the sparsifier,
and then recover small portions of spectrally-critical off-tree edges to the
spanning to significantly improve the approximation quality. However, it
is not clear how many off-tree edges should be recovered for achieving a
desired spectral similarity level within the sparsifier. Motivated by recent
graph signal processing techniques, this paper proposes a similarityaware spectral graph sparsification framework that leverages an efficient
off-tree edge filtering scheme to construct spectral sparsifiers with guaranteed spectral similarity (relative condition number) level. An iterative
graph densification framework and a generalized eigenvalue stability
checking scheme are introduced to facilitate efficient and effective filtering
of off-tree edges even for highly ill-conditioned problems. The proposed
method has been validated using various kinds of graphs obtained from
public domain sparse matrix collections relevant to VLSI CAD, finite
element analysis, as well as social and data networks frequently studied
in many machine learning and data mining applications.

I. I NTRODUCTION
Spectral methods are playing increasingly important roles in many
graph and numerical applications [25], such as scientific computing
[23], numerical optimization [2], data mining [14], graph analytics
[9], machine learning [6], graph signal processing [17], and VLSI
computer-aided design [7], [26]. For example, classical spectral graph
partitioning (data clustering) algorithms embed original graphs into
low-dimension space using the first few nontrivial eigenvectors of
graph Laplacians and subsequently perform graph partitioning (data
clustering) on the low-dimensional graphs to obtain high-quality
solution.
To further push the limit of spectral methods for large graphs,
mathematics and theoretical computer science researchers have extensively studied many theoretically-sound research problems related
to spectral graph theory [1], [3], [10], [14], [19], [22]. Recent spectral
graph sparsification research allows computing nearly-linear-sized
subgraphs (sparsifiers) that can robustly preserve the spectrum (i.e.,
eigenvalues and eigenvectors) of the original graph’s Laplacian,
which immediately leads to a series of “theoretically nearly-lineartime” numerical and graph algorithms for solving sparse matrices,
graph-based semi-supervised learning (SSL), spectral graph partitioning (data clustering), and max-flow problems [2], [8], [10], [11],
[18], [19], [22], [23]. For example, sparsified circuit networks allow
for developing more scalable computer-aided (CAD) design algorithms for designing large VLSI systems [7], [26]; sparsified social
(data) networks enable to more efficiently understand and analyze
large social (data) networks; sparsified matrices can be immediately
leveraged to accelerate the solution computation of large linear
system of equations [27]. To this end, a nearly-linear time spectral
sparsification algorithm leveraging an edge sampling scheme that sets
sampling probabilities proportional to edge effective resistances (of
the original graph) has been proposed in [19]. However, it becomes
a chicken-and-egg problem since calculating effective resistances
requires solving the original Laplacian matrix multiple times and thus
can be extremely expensive for large graphs.

A practically-efficient, nearly-linear complexity spectral graph
sparsification algorithm has been recently introduced in [7], which
first extracts a “spectrally critical” spanning tree subgraph as a
backbone of the sparsifier, and subsequently recovers a small portion
of dissimilar “spectrally critical” off-tree edges to the spanning tree.
This algorithm has been shown to be highly efficient for solving
symmetric diagonally dominant (SDD) matrices in many circuit and
finite element analysis tasks. However, in many scientific computing
and graph-related applications, it is important to generate spectral
graph sparsifiers of desired spectral similarity level with minimum
number of edges: introducing too few edges may lead to poor
approximation of the original graph, whereas too many edges can
result in high computational complexity. For example, considering
a preconditioned conjugate gradient (PCG) solver to be used for
solving multiple right-hand-side (RHS) vectors, it is hoped the PCG
solver should convergence to a good solution as fast as possible,
which usually requires the sparsifier (preconditioner) to have a high
spectral similarity; on the other hand, in many graph partitioning
tasks [20], only the Fielder vector (the first nontrivial eigenvector) of
graph Laplacian is needed, which usually only needs a much sparser
subgraph with relatively low spectral similarity.
This paper introduces a similarity-aware spectral graph sparsification framework that leverages an efficient off-tree edge filtering
scheme to construct spectral sparsifiers with guaranteed spectral
similarity. An incremental graph densification procedure and a generalized eigenvalue stability checking scheme are also proposed
for efficiently and effectively filtering off-tree edges. The proposed
methods will allow sparsifying real-world large graphs according
to desired spectral similarity levels, making it possible to flexibly
trade off computational complexity and solution quality by effectively
controlling the graph sparsity and spectral similarity of the spectral
sparsifiers. The major contributions of this work include:
1) We present a similarity-aware spectral graph sparsification
framework leveraging an iterative edge filtering scheme motivated by recent graph signal processing problems [16], [17].
2) An iterative graph densification procedure is introduced to
incrementally improve the approximation quality of the graph
sparsifier, which allows flexibly trading off the sparsity and
spectral similarity of the spectral sparsifier in various numerical
and graph-related applications.
3) Extensive experiments have been conducted to validate the
proposed method in various numerical and graph-related applications, such as solving sparse SDD matrices, and spectral
graph partitioning, and reduction of large social, and data
networks.
The rest of this paper is organized as follows. Section II provides
a brief introduction to prior work on spectral graph sparsification
problems. In Section III, a similarity-aware method for spectral
graph sparsification is described in details. Section IV demonstrates
extensive experiment results for a variety of real-world, large-scale
sparse Laplacian matrix and graph problems, which is followed by
the conclusion of this work in Section V.

0.04

0.03

0.02

0.01

0

-0.01
0.04

-0.02
-0.03

-0.02

0.02
0

0
-0.02

Fig. 1. Spectral drawing of two σ-spectrally similar (σ = 5) airfoil graphs.

II. S PECTRAL G RAPH S PARSIFICATION
Consider a graph G = (V, E, w) with V denoting the vertex (data
point) set of the graph, E denoting the edge set of the graph, and w
denoting a weight (similarity) function that assigns positive weights
to all edges. The graph Laplacian of G is a symmetric diagonally
dominant (SDD) matrix that can be constructed as follows:

−w(p, q)
if (p, q) ∈ E


 P
w(p,
t)
if (p = q)
LG (p, q) =
(1)
(p,t)∈E



0
otherwise .
Spectral graph sparsification [22] aims to preserve original graph
spectrum within ultra-sparse subgraphs (sparsifiers), which allows
preserving not only eigenvalues and eigenvectors of the original graph
Laplacian, but also distances (e.g. effective resistances) between
vertices, low-dimensional graph embedding, as well cuts in the graph.
Two graphs G and P are σ−spectrally similar if their quadratic
forms satisfy the following conditions for all real vectors x ∈ RV :
xT LP x
≤ xT LG x ≤ σxT LP x.
(2)
σ
Define the relative condition number to be κ(LG , LP ) =
λmax /λmin , where λmax and λmin denote the largest and smallest
nonzero generalized eigenvalues satisfying:
LG u = λLP u,

(c) similarity-aware sparsifier improvement (graph densification) by
gradually adding small portions of “filtered” off-tree edges to the
spanning tree. Fig. 1 shows the spectral drawing [9] of an airfoil
graph [5] as well as its spectrally-similar subgraph computed by the
proposed similarity-aware spectral sparsification algorithm.
In the following sections, assume that G = (V, E, w) is a
weighted, undirected and connected graph, whereas P = (V, Es , ws )
is its graph sparsifier. To simplify the our analysis, we assume the
edge weights in the sparsifier remain the same as the original ones,
though edge weight re-scaling schemes [22] also can be applied to
further improve the approximation quality. The descending eigenvalues of L+
P LG are denoted by λmax = λ1 ≥ λ2 ≥ · · · ≥ λn ≥ 1,
where L+
P denotes the Moore-Penrose pseudoinverse of LP .
B. Off-Tree Edge Embedding via Spectral Perturbation Analysis
1) Spectral Distortion of Spanning-Tree Sparsifiers.: It has been
shown that there are not too many large generalized eigenvalues
for spanning tree sparsifiers [24]: L+
P LG has at most k generalized
eigenvalues greater than stPk(G) , where stP (G) is the total stretch of
the spanning-tree subgraph P with respect the original graph G that
measures the spectral distortion due to the tree approximation. In fact,
recent research shows that every graph has a low-stretch spanning tree
(LSST) such that [18]:
O(m log n log log n) ≥ stP (G) = T r(L+
P LG ) =

III. S IMILARITY-AWARE S PECTRAL S PARSIFICATION OF G RAPHS
A. Overview of Our Approach
The overview of the proposed method for similarity-aware spectral
sparsification of large graphs has been summarized as follows. For
a given input graph, the following key procedures are involved in
the proposed algorithm flow: (a) low-stretch spanning tree extraction
based on its original graph Laplacian; (b) spectral (generalized
eigenvalue) embedding and criticality ranking of off-tree edges
leveraging the recent spectral perturbation analysis framework [7];

λi ≥ σ 2 ,

i=1

(4)
where m = |E| and n = |V |. As a result, it is possible to construct an
ultra-sparse yet spectrally similar sparsifier by recovering only a small
portion of “spectrally critical” off-tree edges to the spanning tree: for
log n
example, σ-similar spectral sparsifiers with O( m log nσlog
) off2
tree edges can be constructed in nearly-linear time [7].
2) Identification of Spectrally Critical Off-Tree Edges.: To identify
“spectrally critical” off-tree edges, the following generalized eigenvalue perturbation analysis is considered:
LG (ui + δui ) = (λi + δλi ) (LP + δLP ) (ui + δui ) ,

(3)

with u being the generalized eigenvector of λ. It can be further
shown that κ(LG , LP ) = σ 2 , which indicates that a smaller relative
condition number or σ 2 corresponds to a higher spectral similarity.
Spectral graph sparsification allows to construct spectrally-similar
subgraphs for mimicking the original graph so that they can be
leveraged for developing much faster numerical and graph algorithms.
For complete graphs, Ramanujan graphs are the best spectral sparsifiers [18], whereas for general graphs the linear-sized sparsifiers
(e.g. Twice-Ramanujan sparsifier) are the best but will be very costly
to compute in practice [1], [12]. Recent nearly-linear-time spectral
sparsification methods aim to find spectral sparsifier through the
following two phases [7], [10], [11], [18], [21], [23]: phase (a)
extracts a low-stretch spanning tree (LSST) P from the original graph
G, while phase (b) adds a few “spectrally critical” off-tree edges to
the tree to form the final nearly-linear-sized spectral sparsifier.

n
X

(5)

where a perturbation δLP is applied for the inclusion of extra offtree edges into LP , resulting in perturbed generalized eigenvalues
and eigenvectors λi + δλi and ui + δui for i = 1, ..., n, respectively.
The key to effective spectral sparsification is to identify the most
“spectrally-critical” off-tree edges that will result in the greatest
reduction in δλi . To this end, the following scheme for embedding
generalized eigenvalues into each off-tree edge is adopted in this work
[7]:
n
P
Step 1: Start with an initial random vector h0 =
αi ui , where ui
i=1

for i = 1, ..., n are the LP -orthogonal generalized eigenvectors of
T
L+
P LG that satisfy ui LP uj = 1 for i = j, and 0 for i 6= j;
Step 2: Perform t-step generalized power iterations with h0 to obtain
n
t
P
ht = L+
h0 =
αi λti ui ;
P LG
i=1

Step 3: Compute the Laplacian quadratic form for δLP,max =
LG − LP with ht :
n
2
P
QδLP,max (ht ) = ht T δLP,max ht =
αi λti (λi − 1)
=

P
(p,q)∈E\Es

wp,q

n
P
i=1

i=1
2

αi2 λ2t
ui T ep,q
i

=

P

wp,q v p,q 2 ,

(p,q)∈E\Es

(6)
where δLP,max includes all off-tree edges, and ep,q ∈ Rn is a
vector with only the p-th element being 1, the q-th element being
−1 and others being 0. In fact, (6) reflects the spectral similarity
between graphs P and G: larger QδLP,max (ht ) values indicates

C. “Spectrally-Unique” Off-Tree Edges
To simplify the following analysis, a “spectrally unique” off-tree
edge ei is defined to be the edge that connects to vertices pi and
qi , and only impacts a single large generalized eigenvalue λi , though
each off-tree edge can usually influence more than one eigenvalue or
eigenvector according to (6). Then the following truncated expansion
of the Laplacian quadratic form of (6) including the top k dominant “spectrally unique” off-tree edges for fixing the top k largest
eigenvalues of L+
P LG can be written as:

v1

“Resonant Frequency”

greater σ 2 and thus lower spectral similarity. More importantly,
(6) allows embedding generalized eigenvalues into the Laplacian
quadratic form of each off-tree edge and thus ranking off-tree edges
according to their “spectral criticality” levels: recovering the off-tree
edges with largest edge Joule heat wp,q v p,q 2 or |v p,q | values will
most significantly decrease the largest generalized eigenvalues. In
practice, using a small number (e.g. 0 < t < 3) of power iterations
will suffice for embedding purpose.

W1

v2

0

W 2 ! W1

W3 ! W2

v3

v4

W4 ! W3

Fig. 2. The eigenvectors (vi ) of increasing eigenvalues (τi ) for a line graph.

translation, modulation, dilation, and down-sampling to the graph
setting has been provided in [17].
The proposed similarity-aware spectral sparsification method aims
to maintain a simplest subgraph sufficient for preserving the slowlyk
k
2
P
P
varying or “low-frequency” signals on graphs, which therefore can be
QδLP,max (ht ) ≈
wpi ,qi αi2 λ2t
ui T epi ,qi =
αi2 λ2t
i
i (λi − 1).
i=1
i=1
regarded as a “low-pass” graph filter. In other words, such spectrally
(7) sparsified graphs will be able to preserve the eigenvectors associate
Since each off-tree edge can only impact one generalized eigenvalue, with low eigenvalues more accurately than high eigenvalues, and
we can assume epi ,qi = γi LP ui according to (6), leading to:
thus perform well for representing “low-frequency” graph signals

but poorly for highly-oscillating (signal) components due to the
γi , i = j,
T
uj epi ,qi =
(8) significant amount of missing edges.
0, i 6= j.
In practice, preserving the spectral (structural) properties of the
Consequently, the effective resistance of edge ei in P becomes:
original graph within the spectral sparsifier is becoming key to design
2 T
2
(9) of many fast numerical and graph-related algorithms [2], [11], [18],
Reefi f = eTpi ,qi L+
P epi ,qi = γi ui LP ui = γi ,
[19], [23]. For example, when using spectral sparsifier as a preconwhich immediately leads to:
ditioner in preconditioned conjugate gradient (PCG) iterations, the
k
k
convergence rate (the number of PCG iterations n < O(σ log )) is
X
X
ef f
αi2 λ2t
=
αi2 λ2t
QδLP,max (ht ) ≈
only governed by the spectral similarity (or relative condition number)
i wpi ,qi Rei
i (λi − 1).
i=1
i=1
for achieving a desired accuracy level , while in spectral graph
(10)
partitioning and data clustering tasks only the first few eigenvectors
Since the stretch of off-tree edge ei is computed by stP (ei ) =
associated with the smallest nontrivial eigenvalues of graph Laplacian
wpi ,qi Reefi f , (10) also indicates that stP (ei ) = λi − 1 holds for
are important [14], [20].
“spectrally-unique” off-tree edges. Consequently, the most “spectrally
critical” off-tree edges identified by (6) or (10) will have the largest E. Off-Tree Edge Filtering with Joule Heat
stretch values and therefore most significantly impact the largest
To only include the edges that are most critical for achieving the
eigenvalues of L+
P LG . In fact, (10) also can be considered as a
+
desired
spectral similarity level, we propose the following scheme for
randomized version of the trace T r(LP LG ) that is further scaled
truncating “spectrally-unique” off-tree edges based on each edge’s
for
each
λ

1.
up by a factor of λ2t
i
i
Joule heat. The theorem in [24] indicates that for a spanning-tree
D. Spectral Sparsification: A Graph Filter Perspective
preconditioner, there will be at most k generalized eigenvalues that
Although (6) and (10) provide a “spectral-criticality” ranking for are greater than stP (G)/k. Therefore, the following simple yet
each off-tree edge, it is not clear how many off-tree edges should be nearly worst-case generalized eigenvalue distribution can be assumed
selected for achieving a desired spectral similarity level. To this end, when the total stretch stP (G) is dominated by the largest eigenvalue
we introduce a simple yet effective off-tree edge filtering strategy λmax = λ1 = stP (G)/2:
motivated by recent graph signal processing problems [16], [17]. To
stP (G)
2λmax
more efficiently analyze signals on general undirected graphs, graph
λi =
=
,i ≥ 1
(11)
i
+
1
i+1
signal processing techniques have been extensively studied recently
[16], [17]. There is a clear analogy between traditional signal process- To most economically select the top-k “spectrally-unique” off-tree
ing based on classical Fourier analysis and graph signal processing: edges that will dominantly impact the top-k largest eigenvalues, the
1) the signals at different time points in classical Fourier analysis following sum of quadratic forms (Joule heat levels) can be computed
correspond to the signals at different nodes in an undirected graph; using (10) by performing t-step generalized power iterations using r
2) the more slowly oscillating functions in time domain correspond to multiple random vectors ht,1 , ..., ht,r , for λi  1:
the graph Laplacian eigenvectors associated with lower eigenvalues

2t+1
r X
k
X
2λmax
and more slowly varying (smoother) components across the graph.
QδLP,max (ht,1 , ..., ht,r ) ≈
(αi,j )2
. (12)
For example, the first few nontrivial eigenvectors associated with the
i+1
j=1 i=1
smallest non-zero eigenvalues of a line graph Laplacian have been
illustrated in Fig. 2, where the increasing eigenvalues correspond to Consider that we are ready to select the k-th “spectrally-unique” offincreasing oscillation frequencies in the line graph. A comprehensive tree edge for fixing the k-th largest eigenvalue λk so that the final
review of fundamental signal processing operations, such as filtering, spectral similarity will become σ 2 = λ̃λ̃max with λ̃max and λ̃min
min

denoting the largest and smallest eigenvalues of L+
P LG after adding
top-k “spectrally-unique” off-tree edges, then we have:
k = 2λmax /λ̃max − 1.

(13)

When using multiple random vectors for computing (12), it is
r
r
P
P
2
2
expected that
αk,j
≈
α1,j
, which allows us to define the
j=1

j=1

normalized edge Joule heat θk for the k-th “spectrally-unique” offtree edge through the following simplifications:
P

r
2
αk,j

 2
2t+1

 j=1
 λ̃max 2t+1
heatλk
σ λ̃min

=
≈
.
θk =
r
P
 λmax
heatλ1
λmax
2
α1,j

original and sparsified Laplacians, λmin can be easily computed. In
our method, we will only assign value 1 or 0 to each node in graphs
G and P for minimizing the following ratio between quadratic forms:
P
wpq
T
x(p)6=x(q),(p,q)∈E
x LG x
P
,
= min
λmin ≈ min
|x|6=0
|x|6=0
wpq
xT LP x
x(i)∈[0,1]

x(i)∈[0,1] x(p)6=x(q),(p,q)∈Es

(17)
which allows approximately computing λmin . To this end, we initialize all nodes with 0 values and will identify only a small portion of
nodes to be colored with 1 values so that the ratio between quadratic
forms can be minimized. An extreme case is to only select a single
node p to be colored with 1 such that:

j=1

(14)
In practice, (14) can be leveraged as a truncation threshold for achieving a desired spectral similarity level (relative condition number)
during spectral sparsification. Although the above filtering scheme
is derived for “spectrally-unique” off-tree edges, general off-tree
edges can also be filtered using similar strategies. Since adding
“spectrally-critical” off-tree edges will mainly impact the largest
generalized eigenvalues but not the smallest ones, we can always
assume λ̃min ≈ λmin and thus use the following edge truncation
scheme for filtering general off-tree edges to achieve a final spectral
similarity level σ: the edge (p, q) will be included into the sparsifier if
its normalized Joule heat value is greater than the truncation threshold
formulated below:
 2
2t+1
heat(p,q)
σ λmin
> θσ ≈
,
(15)
θ(p,q) =
heatmax
λmax
where heatmax is the maximum Joule heat of all off-tree edges
computed by (6) with multiple initial random vectors.
F. Estimation of Extreme Eigenvalues λmax and λmin
To achieve the above edge filtering scheme, we need to compute
(15) that requires to estimate the extreme eigenvalues λmax and λmin
of L+
P LG . In this work, we propose the following efficient methods
for computing these extreme generalized eigenvalues.
1) Estimating λmax Using Power Iterations : It has been shown
that the largest eigenvalues of L+
P LG are well separated from
each other [24], which allows computing rather accurate largest
generalized eigenvalue (λmax or λ1 ) by performing only a small
number of generalized power iterations with an initial random vector.
Since generalized power iterations can converge at a geometric rate
determined by λ2 /λ1 , the error of the estimated largest generalized
eigenvalue will quickly drop to |λ2 /λ1 |k e0 after k iterations for an
initial error e0 . Our results show that only performing five to ten
iterations will suffice for obtaining a good estimation of λ1 .
2) Estimating λmin Using Node Coloring: Since the smallest
eigenvalues of L+
P LG are crowded together [24], using (shifted)
inverse power iterations may not be feasible due to the extremely slow
convergence rate considering the very poor separation of smallest
eigenvalues. To the extent of our knowledge, none of existing
eigenvalue decomposition methods can efficiently compute λmin . To
more efficiently estimate the smallest generalized eigenvalue without
loss of accuracy, this work exploits the following Courant-Fischer
theorem for generalized eigenvalue problems:
T

λmin = min

|x|6=0

x LG x
,
xT LP x

(16)

where x is also required to be orthogonal to the all-one vector to
avoid zero-value denominator. (16) indicates that if we can find a
vector x that minimizes the ratio between the quadratic forms of the

λmin ≈ min
p∈V

D(p)
,
Ds (p)

(18)

where D and Ds are the diagonal vectors satisfying D(p) =
LG (p, p) and Ds (p) = LP (p, p). (18) motivates to sort nodes
and subsequently color
according to the node degree ratio DD(p)
s (p)
a few (e.g. 10% of total) nodes that have the smallest degree ratios
with value 1. The resultant vector x can be used to compute λmin
using (16). The proposed approach only involves a partial sort of
nodes according to their node degree ratios and thus will be highly
scalable for large problems. Our experimental results on real-world
graphs show that the proposed method is highly efficient and accurate
for estimating the smallest generalized eigenvalues when compared
with existing generalized eigenvalue methods [15].
G. Iterative Graph Densification
To achieve more effective edge filtering for similarity-aware spectral graph sparsification, we propose to iteratively recover off-tree
edges to the sparsifier through an incremental graph densification
procedure. Each densification iteration adds a small portion of “filtered” off-tree edges to the latest spectral sparsifier, while the spectral
similarity is estimated to determine if more off-tree edges are needed.
The i-th graph densification iteration includes the following steps:
1) Estimate the spectral similarity by computing λmax and λmin
using the methods in Section III-F;
2) If the spectral similarity is not satisfactory, continue with the
following steps; otherwise, stop the densification procedure.
3) Perform t-step generalized power iterations with a few (e.g.
O(log |V |)) random vectors to compute the sum of Laplacian
quadratic forms similar to (12);
4) Rank each off-tree edge according to its “spectral criticality”
based on (6);
5) Use the truncation threshold (15) to recover top critical yet
dissimilar off-tree edges to the sparsifier.
IV. E XPERIMENTAL RESULTS
The proposed similarity-aware spectral sparsification method has
been implemented in C++ 1 . Test cases demonstrated in this paper
cover a great variety of graph Laplacian problems obtained from
the sparse matrix collection of the University of Florida [5]. All
experiments are performed using a single CPU core of a computing
platform running 64-bit RHEW 7.2 with a 2.67GHz 12-core CPU.
A. Effectiveness of Similarity-Aware Spectral Graph Sparsification
The results of edge ranking according to Joule heat levels computed
by one-step generalized power iteration using (6) have been shown
in Fig. 3 for two SDD matrices in [5]. The off-tree edges critical
for different spectral similarity levels can be easily identified using
1 Available

at https://sites.google.com/mtu.edu/zhuofeng-graphspar



G2_matrix

Thermal1

TABLE I
R ESULTS OF PCG SOLVER FOR GRAPH L APLACIANS IN [5].

V 2 | 500

V 2 | 500

Graphs
G3 circuit
thermal2
ecology2
tmt sym
paraboli fem

V 2 | 100

V 2 | 100

Not too many large generalized eigenvalues!

Fig.
3. Spectral edge
ranking and filtering by Joule heat of off-tree edges
5HVXOWVRI6SHFWUDO
*UDSK6SDUVLILFDWLRQ
for G2 circuit (left) and T hermal1 (right) test cases [5] with top off-tree
 'PHVKJULG Q  ZLWKXQLWHGJHZHLJKW
edges
highlighted in red rectangles.
±  ZUWRIQRGHV H[WUDHGJHVDGGHGWRWKH³KDLUFRPE´VSDQQLQJWUHH
±

\

OPD[ LVUHGXFHGIURP(WR(;UHGXFWLRQ
,QLWLDO-RXOH
+HDW0DS



\

×10 -3




1

6



 



[






4


0.5

2



[
7RWDO6WUHWFK2 Q Q

)LQDO-RXOH 0
+HDW0DS200

0

200



200

100

100
0






0

|E|
3.0E6
3.7E6
2.0E6
2.2E6
1.6E6

|E50 |
|V |

1.11
1.14
1.14
1.21
1.22

N50
21
20
20
19
18

T50
20s
23s
16s
16s
16s

|E200 |
|V |

1.05
1.06
1.06
1.14
1.09

N200
37
36
40
38
38

T200
8s
9s
5s
4s
3s

converging to the required solution accuracy level in much fewer PCG
iterations, but need to retain more edges in the subgraphs and thus
longer time to compute (set up).
C. A Nearly-Linear Time Spectral Graph Partitioning Algorithm



Q u Q 'PHVK

|V |
1.6E6
1.2E6
1.0E6
0.7E6
0.5E6





200

100

100
0

0

Fig. 4. Joule heat distributions before (middle) and after (rightmost) adding
off-tree edges to the “hair comb” spanning tree (leftmost).

the proposed edge filtering scheme. It is also observed in Fig. 3 that
there are not many large eigenvalues of L+
P LG in both cases, which
matches the theoretical results [24].
We demonstrate the results of spatial Joule heat levels of off-tree
edges computed by (6) for a random vector using one-step generalized
power iteration (t = 1) for a “hair comb” spanning tree shown in
Fig. 4. It is not difficult to show that such a spanning tree can not
well match the top part of the original 2D grid, so the Joule heat
levels of off-tree edges at the top should be much greater than the
ones at the bottom part.
√ In fact,
√ it has been shown that a “hair comb”
spanning
tree
for
a
n
×
n 2D mesh will have a total stretch of
√
O(n n) that is mainly contributed by the off-tree edges with largest
stretch values near the top part of the mesh grid [13]. Using the
proposed similarity-aware spectral graph sparsification framework,
we are able to efficiently identify and recover the most “spectrally
critical” off-tree edges, thereby dramatically reducing the largest
generalized eigenvalues. For example, the “hair comb” spanning tree
of a 200 × 200 mesh grid can be dramatically improved in terms
of spectral similarity by recovering top 400 “spectrally critical” offtree edges: σ 2 is reduced from about 64, 000 to about 100 (640×
reduction), which can also be indicated by the Joule heat distributions
before and after adding these off-tree edges as illustrated in Fig. 4.
B. A Nearly-Linear Time Sparse SDD Matrix Solver
The spectral sparsifier obtained by the proposed similarity-aware
algorithm is also leveraged as a preconditioner in a PCG solver. The
RHS input vector b is generated randomly and the solver is set to
converge to an accuracy level ||Ax − b|| < 10−3 ||b|| for all test
cases. “|V |” and “|E|” denote the numbers of nodes and edges in
the original graph, whereas “|Eσ2 |”, “Nσ2 ” and “Tσ2 ” denote the
number of edges in the sparsifier, the number of PCG iterations
required for converging to the desired accuracy level, and the total
time of graph sparsification for achieving the spectral similarity of
σ 2 , respectively. As observed in all test cases, there are very clear
trade-offs between the graph density, computation time, and spectral
similarity for all spectral sparsifiers extracted using the proposed
method: sparsifiers with higher spectral similarities (smaller σ 2 ) allow

It has been shown that by applying only a few inverse power
iterations, approximate Fiedler vector (uf ) that correspond to the
smallest nonzero eigenvalue of the (normalized) graph Laplacian
matrix can be obtained for spectral graph partitioning purpose with
high solution quality [23]. Therefore, using the spectral sparsifiers
computed by the proposed spectral graph sparsification algorithm can
immediately accelerate the PCG solver for inverse power iterations,
leading to nearly-linear runtime algorithms for graph partitioning
problems [23]. In fact, if the spectral graph sparsifier can already
well approximate the original graph (with a small relative condition
number), its Fiedler vector can be directly used for partitioning the
original graph with high fidelity.
To validate the accelerated graph partitioning algorithm, sparse
matrices from [5] and a few larger 2D mesh graphs synthesized
with random edge weights have been tested. As shown in Table II,
the graphs associated with sparse matrices are partitioned based on
the sign of each node according to the approximate Fiedler vectors
computed by a few steps of inverse power iterations. The direct solver
[4] and the preconditioned iterative solver are invoked within each
inverse power iteration for updating the approximate Fiedler vectors
|V+ |
” denotes the ratio of nodes assigned
uf and ũf , respectively. “ |V−
|
with positive and negative signs according to the approximate Fiedler
vector, and “Rel.Err.” denotes the relative error of the proposed solver
|Vdif |
compared to the direct solver computed by |V
, where |Vdif |
|
denotes the number of nodes with different signs in uf and ũf . “ TD ”
(“ TI ”) and “MD ” (“ MI ”) denote the total solution time (excluding
spectral sparsification time) and memory cost of the direct (iterative)
algorithm. To gain high runtime efficiency and solution quality for
all test cases, we extract sparsifiers with σ 2 ≤ 200.
It can be observed that the proposed preconditioned spectral graph
partitioner only results in a very small portion of nodes (0.07% to
4%) assigned with different signs when comparing with the original
spectral graph partitioner, while achieving significant runtime and
memory savings (4−10×). The approximate Fiedler vector computed
by our fast solver for the test case “mesh 1M” is also illustrated in
Fig. 5, showing rather good agreement with the true solution.
D. Spectral Sparsification of Other Complex networks
As shown in Table III, a few finite element, protein, data and
social networks have been spectrally sparsified to achieve σ 2 ≈ 100
using the proposed similarity-aware method. “Ttot ” is the total time
for extracting the sparsifier, “ λλ̃1 ” denotes the reduction ratio of the
1
largest generalized eigenvalues before and after adding off-tree edges
o
s
into the spanning tree sparsifier, and Teig
(Teig
) denotes the time for
computing the first few eigenvectors of the original (sparsified) graph
Laplacians using Matlab’s eigs function. Since spectral sparsifiers
can well approximate the spectral (structural) properties of the
original graph, the sparsified graphs can be leveraged for accelerating

Fig. 5. The approximate Fiedler vector (left) and its magnitude error (right)
for “mesh 1M”.
TABLE II
R ESULTS OF SPECTRAL GRAPH PARTITIONING WITH THE PROPOSED
SIMILARITY- AWARE SPECTRAL GRAPH SPARSIFICATION ALGORITHM .
Test Cases

|V |

|V+ |
|V− |

TD (MD )

TI (MI )

Rel.Err.

G3 circuit
thermal2
ecology2
tmt sym
paraboli fem
mesh 1M
mesh 4M
mesh 9M

1.6E6
1.2E6
1.0E6
0.7E6
0.5E6
1.0E6
4.5E6
9.0E6

1.35
1.00
1.03
0.99
0.98
1.01
0.99
0.99

52.3s (2.3G)
13.0s (0.9G)
12.1s (0.7G)
10.2s (0.6G)
8.8s (0.4G)
10.2s (0.7G)
49.6s (3.0G)
138.5s (6.9G)

7.6s (0.3G)
3.0s (0.2G)
3.4s (0.2G)
1.9s (0.1G)
2.4s (0.1G)
1.7s (0.2G)
8.2s (0.7G)
13.3s (1.5G)

2.2E-2
6.8E-4
8.9E-3
2.1E-2
3.9E-2
3.3E-3
7.5E-3
7.8E-4

many numerical and graph-related tasks. For example, the original
“RCV-80NN” 80-nearest-neighbor graph can not be clustered on our
server with 50GB memory whereas the sparsified one only requires
a few minutes without sacrificing solution quality.
V. C ONCLUSIONS
For the first time, this paper introduces a similarity-aware spectral
graph sparsification framework that can be immediately leveraged
to develop fast numerical and graph-related algorithms. Motivated
by recent graph signal processing concepts, an iterative graph densification procedure based on spectral filtering of off-tree edges
has been proposed for producing ultra-sparse yet spectrally similar
graph sparsifiers, which enables highly flexible trade-offs between
computation complexity and solution quality in various numerical
and graph-related applications. Extensive experiment results have
confirmed the effectiveness and scalability of an iterative matrix
solver and a spectral graph partitioning algorithm for a variety of
large-scale, real-world graph problems, such as VLSI and finite
element analysis problems, as well as data and social networks.
R EFERENCES
[1] J. Batson, D. Spielman, and N. Srivastava. Twice-Ramanujan Sparsifiers.
SIAM Journal on Computing, 41(6):1704–1721, 2012.
[2] P. Christiano, J. Kelner, A. Madry, D. Spielman, and S. Teng. Electrical
flows, laplacian systems, and faster approximation of maximum flow in
undirected graphs. In Proc. ACM STOC, pages 273–282, 2011.

TABLE III
R ESULTS OF SIMILARITY- AWARE NETWORK SPARSIFICATION .
Test Cases

|V |

|E|

Ttot

|Es |

|E|
|Es|

λ1
λ̃1

o
s
Teig
(Teig
)

Reuters911
pdb1HYS
fe tooth
coAuthorsDBLP
RCV-80NN

1.3E4
3.6E4
7.8E4
3.0E5
1.9E5

1.5E5
2.2E6
4.5E5
1.0E6
1.2E7

0.2s
3.8s
3.0s
7.2s
46.5s

1.6E4
7.1E4
9.8E4
3.6E5
3.3E5

9×
31×
5×
3×
36×

21×
8, 981×
8, 050×
1, 364×
28, 410×

2.1s (0.4s)
10.5s (1.3s)
10.3s (1.6s)
2,047s (36s)
N/A (170s)

[3] M. B. Cohen, J. Kelner, J. Peebles, R. Peng, A. B. Rao, A. Sidford,
and A. Vladu. Almost-linear-time algorithms for markov chains and
new spectral primitives for directed graphs. In Proceedings of the 49th
Annual ACM SIGACT Symposium on Theory of Computing, pages 410–
419. ACM, 2017.
[4] T. Davis.
CHOLMOD: sparse supernodal Cholesky
factorization and update/downdate.
[Online]. Available:
http://www.cise.ufl.edu/research/sparse/cholmod/, 2008.
[5] T. Davis and Y. Hu. The university of florida sparse matrix collection.
ACM Trans. on Math. Soft. (TOMS), 38(1):1, 2011.
[6] M. Defferrard, X. Bresson, and P. Vandergheynst. Convolutional neural
networks on graphs with fast localized spectral filtering. In Advances in
Neural Information Processing Systems, pages 3844–3852, 2016.
[7] Z. Feng. Spectral graph sparsification in nearly-linear time leveraging
efficient spectral perturbation analysis. In Design Automation Conference
(DAC), 2016 53nd ACM/EDAC/IEEE, pages 1–6. IEEE, 2016.
[8] W. Fung, R. Hariharan, N. Harvey, and D. Panigrahi. A general
framework for graph sparsification. In Proc. ACM STOC, pages 71–
80, 2011.
[9] X. Hu, A. Lu, and X. Wu. Spectrum-based network visualization
for topology analysis. IEEE Computer Graphics and Applications,
33(1):58–68, 2013.
[10] A. Kolla, Y. Makarychev, A. Saberi, and S. Teng. Subgraph sparsification
and nearly optimal ultrasparsifiers. In Proc. ACM STOC, pages 57–66,
2010.
[11] I. Koutis, G. Miller, and R. Peng. Approaching Optimality for Solving
SDD Linear Systems. In Proc. IEEE FOCS, pages 235–244, 2010.
[12] Y. T. Lee and H. Sun. An SDP-based Algorithm for Linear-sized
Spectral Sparsification. In Proceedings of the 49th Annual ACM SIGACT
Symposium on Theory of Computing, STOC 2017, pages 678–687, New
York, NY, USA, 2017. ACM.
[13] R. Peng. Algorithm Design Using Spectral Graph Theory. PhD thesis,
Carnegie Mellon University, 2013.
[14] R. Peng, H. Sun, and L. Zanetti. Partitioning well-clustered graphs:
Spectral clustering works. In Proceedings of The 28th Conference on
Learning Theory (COLT), pages 1423–1455, 2015.
[15] Y. Saad. Numerical Methods for Large Eigenvalue Problems: Revised
Edition, volume 66. Siam, 2011.
[16] A. Sandryhaila and J. M. Moura. Big data analysis with signal processing
on graphs: Representation and processing of massive data sets with
irregular structure. IEEE Signal Processing Magazine, 31(5):80–90,
2014.
[17] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst. The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular
domains. IEEE Signal Processing Magazine, 30(3):83–98, 2013.
[18] D. Spielman. Algorithms, graph theory, and linear equations in Laplacian
matrices. In Proceedings of the International Congress of Mathematicians (ICM), volume 4, pages 2698–2722, 2010.
[19] D. Spielman and N. Srivastava. Graph sparsification by effective
resistances. SIAM Journal on Computing, 40(6):1913–1926, 2011.
[20] D. Spielman and S. Teng. Spectral partitioning works: Planar graphs and
finite element meshes. In Foundations of Computer Science (FOCS),
1996. Proceedings., 37th Annual Symposium on, pages 96–105. IEEE,
1996.
[21] D. Spielman and S. Teng. Nearly-linear time algorithms for graph
partitioning, graph sparsification, and solving linear systems. In Proc.
ACM STOC, pages 81–90, 2004.
[22] D. Spielman and S. Teng. Spectral sparsification of graphs. SIAM
Journal on Computing, 40(4):981–1025, 2011.
[23] D. Spielman and S. 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.
[24] D. Spielman and J. Woo. A note on preconditioning by low-stretch
spanning trees. arXiv preprint arXiv:0903.2816, 2009.
[25] S.-H. Teng. Scalable algorithms for data and network analysis. Foundations and Trends R in Theoretical Computer Science, 12(1–2):1–274,
2016.
[26] Z. Zhao and Z. Feng. A spectral graph sparsification approach to scalable
vectorless power grid integrity verification. In Proceedings of the 54th
Annual Design Automation Conference 2017, page 68. ACM, 2017.
[27] Z. Zhao, Y. Wang, and Z. Feng. SAMG: Sparsified Graph Theoretic
Algebraic Multigrid for Solving Large Symmetric Diagonally Dominant
(SDD) Matrices. In Proceedings of the 36th International Conference
on Computer-Aided Design (ICCAD). ACM, 2017.

