Constructing Linear-Sized Spectral Sparsification in Almost-Linear Time

Introduction

Graph sparsification is the procedure of approximating a graph $G$ by a sparse graph $G^{\prime}$ such that certain quantities between $G$ and $G^{\prime}$ are preserved. For instance, spanners are defined between two graphs in which the distances between any pair of vertices in these two graphs are approximately the same ; cut sparsifiers are reweighted sparse graphs of the original graphs such that the weights of every cut between the sparsifiers and the original graphs are approximatedly the same . Since both storing and processing large-scale graphs are expensive, graph sparsification is one of the most fundamental building blocks in designing fast graph algorithms, including solving Laplacian systems , designing approximation algorithms for the maximum flow problem , and solving streaming problems . Beyond graph problems, techniques developed for spectral sparsification are widely used in randomized linear algebra , sparsifying linear programs , and various pure mathematics problems .

where $L_{G}$ and $L_{G^{\prime}}$ are the respective graph Laplacian matrices of $G$ and $G^{\prime}$.

Spielman and Teng presented the first algorithm for constructing spectral sparsification. For any undirected graph $G$ of $n$ vertices, their algorithm runs in $O(n\log^{c}n/\varepsilon^{2})$ time, for some big constant $c$, and produces a spectral sparsifier with $O(n\log^{c^{\prime}}n/\varepsilon^{2})$ edges for some $c^{\prime}\geqslant 2$. Since then, there has been a wealth of work on spectral sparsification. For instance, Spielman and Srivastava presented a nearly-linear time algorithm for constructing a spectral sparsifier of $O(n\log n/\varepsilon^{2})$ edges. Batson, Spielman and Srivastava presented an algorithm for constructing spectral sparsifiers with $O(n/\varepsilon^{2})$ edges, which is optimal up to a constant. However, all previous constructions either require $\Omega\left(n^{2+\varepsilon}\right)$ time in order to produce linear-sized sparsifiers , or $O(n\log^{O(1)}n/\varepsilon^{2})$ time but the number of edges in the sparsifiers is sub-optimal.

In this paper we present the first almost-linear time algorithm for constructing linear-sized spectral sparsification for graphs. Our result is summarized as follows:

Given any integer $q\geqslant 10$ and $0<\varepsilon\leqslant 1/120$. Let $G=(V,E,w)$ be an undirected and weighted graph with $n$ vertices and $m$ edges. Then, there is an algorithm that outputs a $(1+\varepsilon)$-spectral sparsifier of $G$ with $O\left(\frac{qn}{\varepsilon^{2}}\right)$ edges. The algorithm runs in $\widetilde{O}\left(\frac{q\cdot m\cdot n^{5/q}}{\varepsilon^{4+4/q}}\right)$ time.

Graph sparsification is known as a special case of sparsifying sums of rank-1 positive semi-definite (PSD) matrices , and our algorithm works in this general setting as well. Our result is summarized as follows:

Given any integer $q\geqslant 10$ and $0<\varepsilon\leqslant 1/120$. Let $I=\sum_{i=1}^{m}v_{i}v_{i}^{\intercal}$ be the sum of $m$ rank-1 PSD matrices. Then, there is an algorithm that outputs scalers $\{s_{i}\}_{i=1}^{m}$ with $|\{s_{i}:s_{i}\neq 0\}|=O\left(\frac{qn}{\varepsilon^{2}}\right)$ such that

The algorithm runs in $\widetilde{O}\left(\frac{qm}{\varepsilon^{2}}\cdot n^{\omega-1+3/q}\right)$ time, where $\omega$ is the matrix-multiplication constant.

A key ingredient in our algorithm is a novel combination of two techniques used in literature for constructing spectral sparsification: Random sampling by effective resistance of edges , and adaptive construction based on barrier functions . We will present an overview of the algorithm, and the intuitions behind it in Section 2.

Algorithm

We study the algorithm of sparsifying the sum of rank-1 PSD matrices in this section. Our goal is to, for any vectors $v_{1},\cdots v_{m}$ with $\sum_{i=1}^{m}v_{i}v_{i}^{\intercal}=I$, find scalars $\{s_{i}\}_{i=1}^{m}$ satisfying

We will use this algorithm to construct graph sparsifiers in Section 3.

Our construction is based on a probabilistic view of the algorithm presented in Batson et al. . We refer their algorithm BSS for short, and give a brief overview of the BSS algorithm at first.

always holds, . To guarantee this, Batson et al. introduces a potential function

The original BSS algorithm is deterministic, and in each iteration the algorithm finds a rank-1 matrix which maximizes certain quantities. To informally explain our algorithm, let us look at the following randomized variant of the BSS algorithm: In each iteration, we choose a vector $v_{i}$ with probability $p_{i}$, and add a rank-1 matrix

to the current matrix $A$. See Algorithm 1 for formal description.

where we used the fact that $vv^{\intercal}\preceq(v^{\intercal}B^{-1}v)B$ for any vector $v$ and PSD matrix $B$. Similarly, we have that

Our algorithm follows the same framework as Algorithm 1. However, to construct a spectral sparsifier in almost-linear time, we expect that the sampling probability $\{p_{i}\}_{i=1}^{m}$ of vectors (i) can be approximately computed fast, and (ii) can be further “reused” for a few iterations.

For fast approximation of the sampling probabilities, we adopt the idea proposed in : Instead of defining the potential function by (2.2), we define the potential function by

To “reuse” the sampling probabilities, we re-compute $\{p_{i}\}_{i=1}^{m}$ after every $\Theta\left(n^{1-1/q}\right)$ iterations: We show that as long as the sampling probability satisfies

for some constant $C>0$, we can still sample $v_{i}$ with probability $p_{i}$ and get the same guarantee on the potential function. The reason is as follows: Assume that $\Delta_{A}=\sum_{i=1}^{T}\Delta_{A,i}$ is the sum of the sampled matrices within $T=O\left(n^{1-1/q}\right)$ iterations. If a randomly chosen matrix $\Delta_{A,i}$ satisfies $\Delta_{A,i}\preceq\frac{1}{Cq}\left(uI-A\right)$, then by the matrix Chernoff bound $\Delta_{A}\preceq\frac{1}{2}\left(uI-A\right)$ holds with high probability. By scaling every sampled rank-1 matrix $q$ times smaller, the sampling probability only changes by a constant factor within $T$ iterations. Since we choose $\Theta(n/\varepsilon^{2})$ vectors in total, our algorithm only recomputes the sampling probabilities $\Theta\left(n^{1/q}/\varepsilon^{2}\right)$ times. Hence, our algorithm runs in almost-linear time if $q$ is a large constant.

2 Algorithm Description

The algorithm follows the same framework as Algorithm 1, and proceeds by iterations. Initially, the algorithm sets

holds for any $j$. In iteration $j$, the algorithm computes the relative effective resistance of vectors $\{v_{i}\}_{i=1}^{m}$ defined by

We remark that, although exact values of $N_{j}$ and relative effective resistances are difficult to compute in almost-linear time, we can use approximated values of $R_{i}$ and $N_{j}$ instead. It is easy to see that in each iteration an over estimate of $R_{i}$, and an under estimate of $N_{j}$ with constant-factor approximation suffice for our purpose.

Analysis

We analyze Algorithm 2 in this section. To make the calculation less messy, we assume the following:

We always assume that $0<\varepsilon\leqslant 1/120$, and $q$ is an integer satisfying $q\geqslant 10$.

We will show how the potential function evolves after each iteration in Section 3.1. Combing this with the ending condition of the algorithm, we will prove in Section 3.2 that the algorithm outputs a linear-sized spectral sparsifier. We will prove Theorem 1.1 and Theorem 1.2 in Section 3.3.

Assume that the number of samples satisfies

By the description of the sampling procedure, it holds that

and $\lambda_{\max}(z_{i}z_{i}^{\intercal})\leqslant\frac{\varepsilon}{q}$. Moreover, it holds that

it holds by the Matrix Chernoff Bound (cf. Lemma 3.2) that

where the last inequality follows from the condition on $N$. Hence, with probability at least

which implies that $\mathbf{0}\preceq\sum_{i=1}^{N}z_{i}z_{i}^{\intercal}\preceq\frac{1}{2}\cdot I$ and $\mathbf{0}\preceq W\preceq\frac{1}{2}\cdot(uI-A)$. ∎

Now we analyze the change of the potential function after each iteration, and show that the expected value of the potential function decreases over time. By Lemma 3.3, with probability at least $1-\frac{\varepsilon^{2}}{100qn}$, it holds that

Lemma 3.4 below shows how the potential function changes after each iteration, and plays a key role in our analysis. This lemma was first proved in for the case of $q=1$, and was extended in to general values of $q$. For completeness, we include the proof of the lemma in the appendix.

Let $w_{1}w_{1}^{\intercal},\cdots,w_{N_{j}}w_{N_{j}}^{\intercal}$ be the matrices picked in iteration $j$, and define for any $0\leqslant i\leqslant N_{j}$ that

We study the change of the potential function after adding a rank-1 matrix within each iteration. For this reason, we use

Assuming $\mathbf{0}\preceq W_{j}\preceq\frac{1}{2}(u_{j}I-A_{j})$, we claim that

for any $1\leqslant i\leqslant N_{j}$. Based on this, we apply Lemma 3.4 and get that

Putting (3.4) and (3.5) together, we have that

So, it suffices to prove the claim (3.3). Since $vv^{\intercal}\preceq(v^{\intercal}B^{-1}v)B$ for any vector $v$ and PSD matrix $B$, we have that

By the assumption of $W_{j}\preceq\frac{1}{2}(u_{j}I-A_{j})$, it holds that

This proves the first statement of the claim.

2 Analysis of the Approximation Guarantee

In this subsection we will prove that the algorithm produces a linear-sized $(1+O(\varepsilon))$-spectral sparsifier. We assume that the algorithm finishes after $k$ iterations, and will prove that the output $A_{k}$ is a $(1+O(\varepsilon))$-spectral sparsifier. It suffices to show that the condition number of $A_{k}$ is small, which follows directly from our setting of parameters.

The output matrix $A_{k}$ has condition number at most $1+O(\varepsilon)$.

Since the condition number of $A_{k}$ is at most

Now we prove that the algorithm finishes in $O\left(\frac{qn^{3/q}}{\varepsilon^{2}}\right)$ iterations, and picks $O\left(\frac{qn}{\varepsilon^{2}}\right)$ vectors in total.

With probability at least $4/5$, the algorithm finishes in $\frac{10qn^{3/q}}{\varepsilon^{2}}$ iterations.

With probability at least $4/5$, the algorithm chooses at most $\frac{10qn}{\varepsilon^{2}}$ vectors.

Since the algorithm finishes within $k$ iterations if

where the last inequality follows from the fact that

By Lemma 3.3, every picked matrix $W_{j}$ in iteration $j$ satisfies

with probability at least $1-\frac{\varepsilon^{2}}{100qn}$, and with probability $9/10$ all matrices picked in $k=\frac{10qn}{\varepsilon^{2}}$ iterations satisfy the condition above. Also, by Lemma 3.5 we have that

since the initial value of the potential function is at most 1. Therefore, it holds that

where the second last inequity follows from Markov’s inequality and (3.6), and the last inequality follows by our choice of $k$. This proves the first statement.

Let $v_{1},\cdots,v_{z}$ be the vectors sampled by the algorithm, and $v_{j}$ is picked in iteration $\tau_{j}$, where $1\leqslant j\leqslant z$. We first assume that the algorithm could check the ending condition after adding every single vector. In such case, it holds that

Following the same proof as the first part and noticing that in the final iteration the algorithm chooses at most $O(n)$ extra vectors, we obtain the second statement. ∎

3 Proof of the Main Results

Now we analyze the runtime of the algorithm, and prove the main results. We first analyze the algorithm for sparsifying sums of rank-1 PSD matrices, and prove Theorem 1.2.

By Lemma 3.7, with probability at least $4/5$ the algorithm chooses at most $\frac{10qn}{\varepsilon^{2}}$ vectors, and by Lemma 3.6 the condition number of $A_{k}$ is at most $1+O(\varepsilon)$, implying that the matrix $A_{k}$ is a $(1+O(\varepsilon))$-approximation of $I$. These two results together prove that $A_{k}$ is a linear-sized spectral sparsifier.

For the runtime, Lemma 3.7 proves that the algorithm finishes in $\frac{10qn^{3/q}}{\varepsilon^{2}}$ iterations, and it is easy to see that all the required quantities in each iteration can be approximately computed in $\widetilde{O}(m\cdot n^{\omega-1})$ time using fast matrix multiplication. Therefore, the total runtime of the algorithm is $\widetilde{O}\left(\frac{q\cdot m}{\varepsilon^{2}}\cdot n^{\omega-1+3/q}\right)$. ∎

Next we show how to apply our algorithm in the graph setting, and prove Theorem 1.1. Let $L=\sum_{i=1}^{m}u_{i}u_{i}^{\intercal}$ be the Laplacian matrix of an undirected graph $G$, where $u_{i}u_{i}^{\intercal}$ is the Laplacian matrix of the graph consisting of a single edge $e_{i}$. By setting

for $1\leqslant i\leqslant m$, it is easy to see that constructing a spectral sparsifier of $G$ is equivalent to sparsifing the matrix $\sum_{i=1}^{m}v_{i}v_{i}^{\intercal}$. We will present in the appendix almost-linear time algorithms to approximate the required quantities

in each iteration, and this gives Theorem 1.1.

By applying the same analysis as in the proof of Theorem 1.2, we know that the output matrix $A_{k}$ is a linear-sized spectral sparsifier, and it suffices to analyze the runtime of the algorithm.

By Lemma 3.3 and the Union Bound, with probability at least $9/10$ all the matrices picked in $k=\frac{10qn^{3/q}}{\varepsilon^{2}}$ iterations satisfy

On the other hand, notice that it holds for any $1\leqslant j\leqslant n$ that

Hence, we apply Lemma 4.5 and Lemma 4.6 to compute all required quantities in each iteration up to constant approximation in time

Since by Lemma 3.7 the algorithm finishes in $\frac{10qn^{3/q}}{\varepsilon^{2}}$ iterations with probability at least $4/5$, the total runtime of the algorithm is

Acknowledgment

This work was partially supported by NSF awards 0843915 and 1111109. Part of this work was done while both authors were visiting the Simons Institute for the Theory of Computing, UC Berkeley, and the second author was affiliated with the Max Planck Institute for Informatics, Germany. We thank Zeyuan Allen-Zhu, Zhenyu Liao, and Lorenzo Orecchia for sending us their manuscript of and the inspiring talk Zeyuan Allen-Zhu gave at the Simons Institute for the Theory of Computing. Finally, we thank Michael Cohen for pointing out a gap in a previous version of the paper and his fixes for the gap, as well as Lap-Chi Lau for many insightful comments on improving the presentation of the paper.

References

Omitted Proofs

In this subsection we prove Lemma 3.4. We first list the following two lemmas, which will be used in our proof.

Let $A$ and $B$ be positive definite matrices, and $q\geqslant 1$. Then it holds that

By the assumption of $w^{\intercal}Y^{-1}w\leqslant\frac{\varepsilon}{q}$, we have that

Note that $0\preceq D\preceq\frac{\varepsilon}{q}\cdot I$, and

Now for the second inequality. Let $Z=uI-A$. By the Sherman-Morrison Formula (Lemma 4.1), it holds that

By the assumption of $w^{\intercal}Z^{-1}w\leqslant\frac{\varepsilon}{q}$, it holds that

Combing $E\preceq\frac{\varepsilon}{q}\cdot I$ with the assumption that $q\geqslant 10$ and $\varepsilon\leqslant 1/10$, we have that

2 Implementation of the Algorithm

Let $L$ and $\widetilde{L}$ be the Laplacian matrices of graph $G$ and its subgraph after reweighting. Let $A=L^{-1/2}\widetilde{L}L^{-1/2}$, and assume that

Under Assumption 4.3, the following statements hold:

We can construct a matrix $S_{u}$ such that

and $S_{u}=p(A)$ for a polynomial $p$ of degree $O\left(\frac{\log(1/\varepsilon\eta)}{\eta}\right)$.

since $u^{-1}A\preceq(1-\eta)I$. Notice that $u^{-1/2}I\preceq(uI-A)^{-1/2}$, and therefore

Setting $T=\frac{c\log(1/(\varepsilon\eta))}{\eta}$ for some constant $c$ and defining $S_{u}=u^{-1/2}p_{T}(u^{-1}A)$ gives us that

Using the same analysis as before, we have that

Let $A=\sum_{i=1}^{m}v_{i}v_{i}^{\intercal}$, and suppose that $A$ satisfies Assumption 4.3. Then, we can compute $\{r_{i}\}_{i=1}^{m}$ and $\{t_{i}\}_{i=1}^{m}$ in $\widetilde{O}\left(\frac{m}{\varepsilon^{2}\eta}\right)$ time such that

Define $u_{i}=L^{1/2}v_{i}$ for any $1\leqslant i\leqslant m$. By Lemma 4.4, we have that

We apply a nearly-linear time Laplacian solver to compute $\left\|QBp\left(L^{-1}\widetilde{L}\right)L^{-1}u_{i}\right\|^{2}$ for all $\{u_{i}\}_{i=1}^{m}$ up to $(1\pm\varepsilon/10)$-multiplicative error in time $\widetilde{O}\left(\frac{m}{\varepsilon^{2}\eta}\right)$. This gives the desired $\{r_{i}\}_{i=1}^{m}$.

The computation for $\{t_{i}\}_{i=1}^{m}$ is similar. By Lemma 4.4, it holds for any $1\leqslant i\leqslant m$ that

We invoke the Johnson-Lindenstrauss Lemma and a nearly-linear time Laplacian solver as before to obtain required $\{t_{i}\}_{i=1}^{m}$. The total runtime is $\widetilde{O}\left(\frac{m}{\eta\varepsilon^{2}}\right)$. ∎

Under Assumption 4.3, we can compute values $\alpha,\beta$ in $\widetilde{O}\left(\frac{m}{\eta\varepsilon^{3}}\right)$ time such that

By Lemma 4.4, we have that $S_{u}\approx_{\varepsilon/10}(uI-A)^{-1/2}$. Hence, $\lambda_{\max}(S_{u})^{-2}\approx_{3\varepsilon/10}\lambda_{\min}(uI-A)$, and it suffices to estimate $\lambda_{\max}(S_{u})$. Since

Let $z$ be a polynomial defined by $z(x)=xq^{2}(x)$ and $L^{\prime}=(B^{\prime})^{\intercal}(B^{\prime})$. Then, we have that

Applying the same analysis as before, we can estimate the trace in $\widetilde{O}\left(\frac{m}{\eta\varepsilon^{3}}\right)$ time. ∎