Density Independent Algorithms for Sparsifying $k$-Step Random Walks

TL;DR

This paper introduces faster algorithms for sparsifying transition matrices of k-step random walks on graphs, improving efficiency by leveraging density-independent methods and tighter sampling bounds.

Contribution

It presents novel algorithms that produce sparse approximations of k-step random walk matrices with improved runtime and size bounds, based on better sampling analysis.

Findings

Produces a graph with about n log n edges approximating the k-step walk

Achieves runtime of m + n log^4 n for the sparsification process

Revisits and refines density independent graph sparsification techniques

Abstract

We give faster algorithms for producing sparse approximations of the transition matrices of $k$-step random walks on undirected, weighted graphs. These transition matrices also form graphs, and arise as intermediate objects in a variety of graph algorithms. Our improvements are based on a better understanding of processes that sample such walks, as well as tighter bounds on key weights underlying these sampling processes. On a graph with $n$ vertices and $m$ edges, our algorithm produces a graph with about $n\log{n}$ edges that approximates the $k$-step random walk graph in about $m + n \log^4{n}$ time. In order to obtain this runtime bound, we also revisit "density independent" algorithms for sparsifying graphs whose runtime overhead is expressed only in terms of the number of vertices.

Equations

$$L_{G^k}\stackrel{\mathrm{def}}{=}D-A{\left(D^{-1}A\right)}^{k-1}.$$

$$w_{u_0,\dots ,u_k}\stackrel{\mathrm{def}}{=}\frac{{\prod }_{i=1}^k{w}_{u_{i-1},u_i}}{{\prod }_{i=1}^{k-1}{d}_{u_i}}.$$

$$w_{u_0,u_k}^{G^k}\stackrel{\mathrm{def}}{=}{\left[A{\left(D^{-1}A\right)}^{k-1}\right]}_{u_0{u}_k}=\sum _{u_1,\dots ,u_{k-1}}{w}_{u_0,\dots ,u_k}.$$

$$A⪯B$$

$$A{\approx }_{\kappa }B$$

$$L_H{\approx }_{\kappa }{L}_G,$$

$${\mathcal{R}}_{\mathrm{eff}}^G(u,v)\stackrel{\mathrm{def}}{=}{\mathbf{\chi }}_{uv}^T{L}_G^{†}{\mathbf{\chi }}_{uv},$$

$$\sum _{e\in E}{w}_e{\mathcal{R}}_{\mathrm{eff}}^G\left(u,v\right)=n-1.$$

$${\mathcal{R}}_{\mathrm{eff}}^H(u,v)\le {\mathcal{R}}_{\mathrm{eff}}^G(u,v).$$

$${\mathcal{R}}_{\mathrm{eff}}^{G^k}(u,v)\le 2{\mathcal{R}}_{\mathrm{eff}}^G(u,v)$$

$${\mathcal{R}}_{\mathrm{eff}}^{G^k}(u,v)\le {\mathcal{R}}_{\mathrm{eff}}^{G^2}(u,v).$$

$${\mathcal{R}}_{\mathrm{eff}}^G(u_0,u_k)\le \sum _{0\le i<k}{\mathcal{R}}_{\mathrm{eff}}^G(u_i,u_{i+1}).$$

$${\tilde{\mathbf{\tau }}}_{u_0,u_1,\dots ,u_k}^{(k)}\stackrel{\mathrm{def}}{=}{w}_{u_0,u_1,\dots ,u_k}\cdot \sum _{0\le i<k}{\tilde{r}}_{u_i,u_{i+1}}.$$

$$w_{u_0,\dots ,u_k}\cdot \sum _{0\le i<k}{\tilde{r}}_{u_i,u_{i+1}}$$

$$\sum _{i=0}^{k-1}\frac{1}{k}\cdot \frac{w_{u_i,u_{i+1}}{\tilde{r}}_{u_i,u_{i+1}}}{⟨w,\tilde{r}⟩}\cdot \prod _{j=1}^i\frac{w_{u_{j-1},u_j}}{d_{u_j}}\cdot \prod _{j=i+1}^{k-1}\frac{w_{u_j,u_{j+1}}}{d_{u_j}}=\frac{w_{u_0,\dots ,u_k}}{k⟨w,\tilde{r}⟩}\sum _{i=0}^{k-1}{\tilde{r}}_{u_i,u_{i+1}},$$

$$\sum _{\left(u_0,u_1,\dots ,u_k\right)}\left[w_{u_0,u_1,\dots ,u_k}\cdot \left(\sum _{0\le i<k}{\tilde{r}}_{u_i,u_{i+1}}\right)\right]=k\cdot \left(\sum _{e\in E[G]}{w}_e{\tilde{r}}_e\right).$$

$${\mathcal{R}}_{\mathrm{eff}}^G(u,v)+{\mathcal{R}}_{\mathrm{eff}}^G(v,w)\ge {\mathcal{R}}_{\mathrm{eff}}^{G^2}(u,w).$$

$${\mathcal{R}}_{\mathrm{eff}}^{G^2}(u,w)={\mathcal{R}}_{\mathrm{eff}}^{G\times P_2}(u^{(A)},w^{(A)})\le {\mathcal{R}}_{\mathrm{eff}}^{G\times P_2}(u^{(A)},v^{(B)})+{\mathcal{R}}_{\mathrm{eff}}^{G\times P_2}(v^{(B)},w^{(A)}).$$

$${\mathcal{R}}_{\mathrm{eff}}^{G^2}(u,v)={\mathcal{R}}_{\mathrm{eff}}^{G\times P_2}(u^{\left(A\right)},v^{\left(A\right)}),$$

$$\frac{1}{h\cdot w_{u_0,u_1,\dots ,u_k}\cdot {\sum }_{0\le i<k}{\tilde{r}}_{u_i{u}_{i+1}}}.$$

$$w_{u_0,u_1,\dots u_k}\sum _{0\le i<k}{\mathcal{R}}_{\mathrm{eff}}^G\left(u_i,u_{i+1}\right).$$

$${\mathcal{R}}_{\mathrm{eff}}^{G^k}\left(u,v\right)\le {\mathcal{R}}_{\mathrm{eff}}^{G\times P_2}\left(u^{\left(A\right)},v^{\left(A\right)}\right)={\mathcal{R}}_{\mathrm{eff}}^{G\times P_2}\left(u^{\left(B\right)},v^{\left(B\right)}\right).$$

$${\tilde{r}}_{uv}\ge r^{G\times P_2}\left(u^{(A)},v^{(B)}\right).$$

$$st{r}_T(e)\stackrel{\mathrm{def}}{=}{w}_e\sum _{e^{\prime }\in {\mathcal{P}}_T(e)}\frac{1}{w_{e^{\prime }}}$$

$$st{r}_T(G)\stackrel{\mathrm{def}}{=}\sum _{e\in E}st{r}_T(e)$$

$${\hat{G}}^{(i)}\stackrel{\mathrm{def}}{=}\hat{G}+2^i\cdot T,$$

$$\tilde{r}\stackrel{\mathrm{def}}{=}\kappa \cdot {\mathcal{R}}_{\mathrm{eff}}^T\ge {\mathcal{R}}_{\mathrm{eff}}^G.$$

$$⟨w,\tilde{r}⟩=\kappa \cdot st{r}_T(G)\le \kappa \cdot st{r}_T(G^{\prime })=\hat{O}(\kappa \left|E(G^{\prime })\right|\log n).$$

$$w\left(u_0,u_1,\dots u_k\right)=w\left(u_0,u_1,\dots u_{k-1}\right)\frac{A_{u_{k-1}{u}_k}}{d_{u_{k-1}}}.$$

$$\sum _{\begin{array}{c}(u_0,u_1,\dots u_k)\\ e=u_i{u}_{i+1}\end{array}}w\left(u_0,u_1,\dots u_k\right)$$

Full text

1 Introduction

Random walks in graphs are fundamental objects in both graph algorithms and graph data structures. Problems related to random walks, such as shortest paths and minimum cuts are well studied in both static [Som14] and dynamic settings [HKN14, GHT16]. While some of these problems, such as shortest path, aim to find a single walk, other problems such as flows/cuts [GT14] or triangle densities [BPWZ14, Tso08] aim to capture information related to collections of walks. Algorithms and data structures for such problems often need to store, or can be sped up by, intermediate structures that capture the global properties of multi-step walks [PT12, GP13, ADK*+*16, BHI14]. However, many intermediate structures are inherently dense and therefore expensive to compute explicitly.

Graph sparsification is a technique for efficiently approximating a dense graph by a sparser one, while preserving some key properties such as sizes of graph cuts, distances between vertices, or linear operator properties of matrices associated with the graphs. Spectral sparsifiers provide linear operator approximations that also imply approximation to all graph cuts. Their constructions have some of the simplest interactions with statistical concentration results [BSST13]. Spectral sparsifiers of (possibly dense) intermediate objects have a variety of applications in graph algorithms, such as sampling from graphical models [CCL*+*15], solving linear systems [PS14, KLP*+*16], sampling random spanning trees [DKP*+*16] and maintaining approximate minimum cuts in dynamically changing graphs [ADK*+*16]. In these applications, the optimal performance is achieved by producing a sparsifier of the denser intermediate object directly, instead of generating the larger exact object. Such implicit sparsification routines were first studied for random walk matrices [PS14, CCL*+*15]. These matrices contain the pairwise transition probabilities between vertices under $k$-step walks. Moreover, such matrices are dense even for sparse original graphs with small $k$: the $2$-step walk on the $n$-vertex star has non-zero transition probabilities between any pair of vertices. On the other hand, as the $k$-step random walk can be viewed as a single random process, these vertex-to-vertex transition probabilities correspond to a graph, and therefore have a sparse approximate.

Cheng et al. [CCL*+*15] systematically studied random walk sparsification and its applications. They gave a routine that produces an $\epsilon$-spectral sparsifier (which we will formally define in Subsection 2.2) with $O(\epsilon^{-2}n\log{n})$ edges for a $k$-step walk matrix in $O(\epsilon^{-2}k^{2}m\log^{O(1)}n)$ time. Our main result, which we show in Section 3 is a direct improvement of that routine:

Theorem 1.1**.**

(Sparsifying Laplacian Monomials) Given a graph $G$ and an error $\varepsilon\in(0,1)$, there is an algorithm that outputs an $\varepsilon$-spectral sparsifier of $G^{k}$ with at most $O(\varepsilon^{-2}n\log n)$ edges in $\widehat{O}(m+k^{2}\varepsilon^{-2}n\log^{4}n)$ time. 111 We use $\widehat{O}$ to denote the omission of logarithmic terms lower than the ones shown in the set. In all cases in this paper, we track terms of $\log{n}$ explicitly and such notation hides terms of $\log\log{n}$. In all these cases, this notation hides a term of at most $(\log\log{n})^{2}$.

We term this type of running time with most of the overhead on the number of vertices, $n$, as density independent. Such runtimes arise naturally in many other graph problems [FT87], and was first studied for graph sparsification in an earlier manuscript by a subset of the authors [JK15]. Our results can also be combined with the repeated-squaring technique in [CCL*+*15] to reduce the runtime dependence on $k$ to logarithmic [CC16]. This plus generalizations to general random walk polynomials [CC16] would then supersede all claims from [JK15]. As these steps are much closer to [CCL*+*15], we will focus on the small $k$ case in this paper. Furthermore, as our sparsification algorithm has a much more direct interaction with routines that provide upper bounds of effective resistances, they can likely be combined with tools from [ADK*+*16] to give dynamic algorithms for maintaining $G^{k}$ under insertions/deletions to $G$. However, as there are currently only few applications of such sparsifiers, we believe it may be more fruitful to extend the applications before further developing the tools.

Our algorithms, as with the ones from [JK15, CCL*+*15] are based on implicit sampling of dense graphs by probabilities related to effective resistances. Such approach is the only known efficient approach for even the ‘simpler’ problem of producing cut sparsifiers of $G^{k}$. Our improvements rely on an a key insight from the sparse Gaussian elimination algorithm by Kyng and Sachdeva [KS16]: using triangle inequality between effective resistances to obtain a tighter set of probability upper bounds. This allows us to select the first edge, from which we “grow” a length $k$ walk, via an adaptive sampling process, instead of using uniform sampling as in the previous result [CCL*+*15]. Furthermore, this adaptive process removes any sampling count dependencies on $m$, the number of edges, making a density-independent runtime possible. This type of running time also has analogs in input sparsity time algorithms in randomized numerical linear algebra [W*+*14, CLM*+*15, CMM17]. Obtaining density-independent bounds is critical for graph sparsification algorithms because they are primarily invoked on relatively dense graphs. A graph sparsification routine that produces a sparsifier with $\widehat{O}(n\log^{2}{n})$ edges in $\widehat{O}(m\log^{2}{n})$ time, such as the combinatorial algorithm given in [KPPS17], will only be invoked when $m>n\log^{2}{n}$, which means the running time of the algorithm is actually $\Omega(n\log^{4}{n})$. As a result, we believe that for graph sparsification to work as a primitive for processing large graphs, a running time of $\widehat{O}(m+n\log^{2}{n})$ or better is necessary.

In Section 4, we provide some steps toward this direction by giving a better density-independent spectral sparsification algorithm. We combine ideas from previous density-independent algorithms for sparsifying graphs [KLP15] with recent developments in tree embedding and numerical algorithms to obtain numerical sparsification routines that run in $\widehat{O}(m+n\log^{4}{n})$ time, and combinatorial ones that take $\widehat{O}(m+n\log^{6}{n})$ time. Both of these routines are in turn applicable to the walk sparsification algorithm in Section 3, giving routines for sparsifying $k$-step walks with similar running times: the bound stated in Theorem 1.1 is via the numerical routine. While these results are far from what we think is the best possible, we show a variety of new algorithmic tools for designing algorithms for sparsifying $k$-step random walks.

2 Background

We start with some background information about graphs and matrices corresponding to them. These matrices allow us to define graph approximations, as well as compute key sampling probabilities needed to produce spectral sparsifiers. Due to space constraints, we will only formally define most of the concepts. More intuition on them can be found in notes on spectral graph theory and random walks such as [DS84, Lov93].

2.1 Random Walks and Matrices

Let $G=(V,E,\boldsymbol{\mathit{w}})$ be a weighted undirected graph. We define its adjacency matrix $\boldsymbol{\mathit{A}}$ as $\boldsymbol{\mathit{A}}_{uv}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{\mathit{w}}_{uv}$, and its degree matrix $\boldsymbol{\mathit{D}}$ as $\boldsymbol{\mathit{D}}_{uu}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{v\in V}w_{uv}$ and $\boldsymbol{\mathit{D}}_{uv}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}0$ when $u\neq v$. This leads to the graph Laplacian $\boldsymbol{\mathit{L}}_{G}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}}$.

One step of a random walk can be viewed as distributing the ‘probability mass’ at a vertex evenly among the edges leaving it, and passing them onto its neighbors. In terms of these matrices, it is equivalent to first dividing by $\boldsymbol{\mathit{D}}$, and then multiplying by $\boldsymbol{\mathit{A}}$. Thus, the transition matrix of the $k$th step random walk is given by $(\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{A}})^{k}$. The corresponding Laplacian matrix of the $k$-step random walk is defined by

[TABLE]

The matrices $\boldsymbol{\mathit{A}}(\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{A}})^{k-1}$ can be viewed as a sum over length $k$ walks. This view is particularly useful in our algorithm, as well as the earlier walk sparsification algorithm by Cheng et al. [CCL*+*15] because these walks are a more ‘natural’ unit upon which sparsification by effective resistances is applied. Formally, we can define the weight of a length $k$ walk $u_{0},\dots,u_{k}$ by

[TABLE]

Straightforward checking shows that for any $u_{0},u_{k}\in V$, it holds that

[TABLE]

2.2 Spectral Approximations of Graphs

Our notion of matrix approximations will be through the $\approx$ symbol, which is in turn defined through the Löewner partial ordering of matrices. For two matrices, $\boldsymbol{\mathit{A}}$, and $\boldsymbol{\mathit{B}}$, we say that

[TABLE]

if $\boldsymbol{\mathit{B}}-\boldsymbol{\mathit{A}}$ is positive semidefinite, and

[TABLE]

if there exists bounds $\lambda_{\min}$ and $\lambda_{\max}$ such that $\lambda_{\min}\boldsymbol{\mathit{A}}\preceq\boldsymbol{\mathit{B}}\preceq\lambda_{\max}\boldsymbol{\mathit{A}}$, and $\lambda_{\max}\leq\kappa\lambda_{\min}$. This notation is identical to generalized eigenvalues, and in particular, $\boldsymbol{\mathit{L}}_{G}\approx_{\kappa}\boldsymbol{\mathit{L}}_{H}$ implies that all cuts on them are within a factor of $\kappa$ of each other.

The adjacency matrix of a graph has several of undesirable properties when it comes to operator based approximations: it can have a large number of eigenvalues at [math], which must be exactly preserved under relative error approximations. As a result, graph approximations are defined in terms of graph Laplacians. As we will discuss below, these approximations are often in terms of reducing edges. So formally, we say that a graph $H$ is a $\kappa$-sparsifier of $G$ if

[TABLE]

and our goal is to compute an $\epsilon$-sparsifier of the $k$-step random walk matrix $\boldsymbol{\mathit{L}}_{G^{k}}$.

2.3 Graph Sparsification by Effective Resistances

There are two ways of viewing graph sparsification: either as tossing coins independently on the edges, or sampling a number of them from an overall probability distribution. We take the second view here because it may prove expensive to access all edges in $G^{k}$. The pseudocode of the generic sampling scheme is given in Algorithm 1.

Algorithmically, the sampling step can be implemented by first generating a number uniformly random in $[0,\sum_{e}\widetilde{\boldsymbol{\tau}}_{e}]$, and binary searching among the prefix sums of the $\widetilde{\boldsymbol{\tau}}_{e}$ values until it reaches the edge corresponding to that point. Note also that if we want to generate random numbers with bounded precision, we can also round the $\widetilde{\boldsymbol{\tau}}_{e}$ values of all edges up to the nearest multiple of $1/n$, leading to at most $m/n=O(n)$ extra edges. The guarantees of this routine require defining effective resistances and leverage scores. Effective resistance is a metric on a graph that is defined by:

[TABLE]

where $\boldsymbol{\mathit{L}}_{G}^{{\dagger}}$ denotes the pseudoinverse of $\boldsymbol{\mathit{L}}_{G}$ and $\boldsymbol{\chi}_{uv}$ is the indicator vector with $1$ at $u$ and $-1$ at $v$.

The effective resistances $\mathcal{R}_{\mathrm{eff}}^{G}$ are directly related to the statistical leverage scores $\boldsymbol{\tau}$ by the relation $\boldsymbol{\tau}_{e}=\boldsymbol{\mathit{w}}_{e}\mathcal{R}_{\mathrm{eff}}^{G}(e)$. Moreover, these scores are well defined for general matrices, and have a wide range of applications in randomized linear algebra [W*+*14, CLM*+*15, CMM17]. The guarantees of sampling by weight times effective resistance, or leverage scores, can then be formalized as:

Lemma 2.1**.**

(Sampling by Upper Bounds on Leverage Scores [Tro12]) Suppose $G=(V,E,w)$ is a graph and $\widetilde{\boldsymbol{\tau}}$ is a vector such that $\widetilde{\boldsymbol{\tau}}_{e}\geq\boldsymbol{\mathit{w}}_{e}\mathcal{R}_{\mathrm{eff}}^{G}(e)$ for every edge $e$, then any process that simulates the ideal sampling in Algorithm 1 produces an $\varepsilon$-sparsifier of $G$ with $O(\varepsilon^{-2}\sum_{e}\widetilde{\boldsymbol{\tau}}_{e}\log{n})$ edges in $\widetilde{O}(\varepsilon^{-2}\sum_{e}\widetilde{\boldsymbol{\tau}}_{e}\log^{2}{n})$ time.

The bound on sample count then follows from:

Fact 2.2** (Foster’s Theorem).**

For any undirected graph $G=(V,E,\boldsymbol{\mathit{w}})$, we have:

[TABLE]

Leverage scores are the preferred objects for defining sampling distributions as they are scale invariant: doubling the weights of all edges does not change them. However, we will still make extensive uses of effective resistances because of the need to approximate them across different graphs. Such approximations are difficult to state for leverage scores because spectrally similar graphs may have very different sets of combinatorial edges.

Fact 2.3**.**

If $G$ and $H$ are graphs such that $\boldsymbol{\mathit{L}}_{G}\preceq\boldsymbol{\mathit{L}}_{H}$, then for any vertices $u$ and $v$ we have

[TABLE]

Note that this generalizes Rayleigh’s monotonicity law, which postulates that the effective resistances can only increase as one removes edges from a graph.

3 Random Walk Sparsification via Walk Sampling

In this section we describe our improved algorithm for sparsifying random walk polynomials. The main difficulty faced by such a routine is that the actual walk matrix cannot be constructed. Instead, we need to simulate the ideal sampling routine shown in Algorithm 1 by constructing nearly tight upper bounds of leverages scores in $G^{k}$ that can also be efficiently sampled from, without having explicit access to $G^{k}$.

The first obstacle to obtain such estimates is to get an access to effective resistances in $G^{k}$. To this end, the following lemma from [CCL*+*15] provides a helpful starting point.

Lemma 3.1**.**

[CCL*+*15]** For odd $k$, we have $\frac{1}{2}\boldsymbol{\mathit{L}}_{G}\preceq\boldsymbol{\mathit{L}}_{G^{k}}\preceq k\boldsymbol{\mathit{L}}_{G}$ and for even $k$, we have $\boldsymbol{\mathit{L}}_{G^{2}}\preceq\boldsymbol{\mathit{L}}_{G^{k}}\preceq\frac{k}{2}\boldsymbol{\mathit{L}}_{G^{2}}.$

Furthermore, note that Lemma 3.1 combined with Fact 2.3 implies for odd $k$ that

[TABLE]

and for even $k$ that

[TABLE]

Since $G^{k}$ might be dense, i.e. $E[G^{k}]=\Theta(n^{2})$, it is prohibitive to use (4) and (5) directly. Instead, we upper bound the values with a walk using the triangle inequality of effective resistances [Spi07, Lemma 9.6.1].

Fact 3.2** (Triangle Inequality for Effective Resistances).**

For any graph $G$ and any walk $(u_{0},u_{1},\ldots,u_{k})$, we have

[TABLE]

Now, suppose we have a vector $\boldsymbol{\widetilde{\mathit{r}}}$ that upper bounds the effective resistances, i.e., $\boldsymbol{\widetilde{\mathit{r}}}_{e}\geq\mathcal{R}_{\mathrm{eff}}^{G}(e)$ for all $e$. Then, by Lemma 2.1 and Fact 3.2, to sparsify $G^{k}$, it suffices to sample a length $k$ random walk in $G$ with probability proportional to

[TABLE]

This distribution has the advantage that it is efficiently computable:

Lemma 3.3**.**

For any graph $G=(V,E,\boldsymbol{\mathit{w}})$, and any vector $\boldsymbol{\widetilde{\mathit{r}}}\in\mathbb{R}^{E}$, we can sample length $k$ walks with probability proportional to

[TABLE]

using the following procedure:

Pick uniformly at random an index $i$ in the range $[0,k-1]$. 2. 2.

Choose an edge $(u_{i},u_{i+1})$ with probability proportional to $\boldsymbol{\mathit{w}}_{e}\boldsymbol{\widetilde{\mathit{r}}}_{e}$. 3. 3.

Extend the walk in both directions from $u_{i}$ and $u_{i+1}$ via two random walks.

Proof.

By total law of probability, the procedure samples a fixed walk $(u_{0},\ldots,u_{k})$ with probability equal to

[TABLE]

where the first term is step (1), the second term is step (2) and the third and the fourth are for the two random walks extending the selected edge. ∎

The total number of samples needed by Lemma 2.1 can be extracted from summing over random walks containing a particular edge in a way similar to [CCL*+*15]. For completeness, we present its proof in Appendix A.

Lemma 3.4**.**

For any weighted graph $G=(V,E,\boldsymbol{\mathit{w}})$, any $k\in\mathbb{N}_{+}$, and any vector $\boldsymbol{\widetilde{\mathit{r}}}\in\mathbb{R}^{E}$, it holds

[TABLE]

For every odd $k$, by setting $\boldsymbol{\widetilde{\mathit{r}}}$ to (an approximation of) $\mathcal{R}_{\mathrm{eff}}^{G}$, yields an efficient sampling procedure due to (8) and Lemma 3.3.

However, when $k$ is even Lemma 3.1 gives a bound in terms of $\mathcal{R}_{\mathrm{eff}}^{G^{2}}$ (not $\mathcal{R}_{\mathrm{eff}}^{G}$), i.e. $\mathcal{R}_{\mathrm{eff}}^{G^{k}}(u,v)\leq\mathcal{R}_{\mathrm{eff}}^{G^{2}}(u,v)$. Hence, the distribution in Lemma 3.3 requires an access to the $2$-step random walk matrix $G^{2}$, which might also be dense and therefore expensive to compute.

Moreover, suppose $G$ is a $2$-length path graph $u-v-w$, then $\mathcal{R}_{\mathrm{eff}}^{G^{2}}(e)=+\infty$ for $e\in G$, since $G^{2}$ has only one edge $(u,w)$ (and self-loops). A naive approach to tackle these issues is to substitute $\mathcal{R}_{\mathrm{eff}}^{G^{2}}$ with $\mathcal{R}_{\mathrm{eff}}^{G}$. However, this approach fails shortly since it is not true in general that

[TABLE]

We work around this by using effective resistances from the “double cover” of $G$, instead. The “double cover” $G\times P_{2}$ is the tensor product of $G$ and a path of length $2$. Combinatorially, $G\times P_{2}$ is a bipartite graph with vertex sets $V^{(A)},V^{(B)}$ each a copy of $V$ such that for every edge $(u,v)\in G$ we insert in $G\times P_{2}$ the following two edges: $u^{(A)}v^{(B)}$ and $u^{(B)}v^{(A)}$ with $\boldsymbol{\mathit{w}}_{u^{(A)}v^{(B)}}=\boldsymbol{\mathit{w}}_{u^{(B)}v^{(A)}}=\boldsymbol{\mathit{w}}_{uv}$.

The next lemma (proved in Appendix A) fixes (9) and guarantees for every edge $(u,w)\in G^{2}$ that

[TABLE]

Lemma 3.5**.**

For any vertices $u$ and $v$ in $G$, it holds

[TABLE]

where $u^{(A)}$ and $v^{(A)}$ are the corresponding copies of $u$ and $v$ in $V^{(A)}$, respectively.

Using the preceding results, we design an algorithm with improved sampling count. It takes any procedure that produces effective resistance distribution that dominates the true one, and produces samples that suffice for simulating the ideal sampling algorithm on $G^{k}$ (c.f. Subsection 2.3, Algorithm 1). The pseudocode for this routine is shown in Algorithm 2.

Note that from the perspective of this framework of picking edges with probabilities proportional to $\boldsymbol{\mathit{w}}_{e}\boldsymbol{\widetilde{\mathit{r}}}_{e}$, and extending them into walks, the previous result [CCL*+*15] can be viewed as utilizing a simple $\mathrm{EREstimator}$ that returns $1$ on the effective resistance of every edge.

Theorem 3.6**.**

Given any graph $G$, any values of $k$ and $\varepsilon$, and any effective resistance estimation algorithm $\mathrm{EREstimator}$ that produces w.h.p. effective resistance that sum up to $f(n,m)$, calling $\mathrm{SparsifyG^{k}}(G,k,\varepsilon,\mathrm{EREstimator})$ produces an $\varepsilon$-sparsifier of $G^{k}$ with $O(\varepsilon^{-2}k\log{n}\cdot f(2n,2m))$ edges in time proportion to the cost of one call to $\mathrm{EREstimator}$ on a graph of twice the size, plus an overhead of $\widetilde{O}(\varepsilon^{-2}k^{2}\log^{2}{n}\cdot f(2n,2m))$.

Proof.

By Lemma 2.1, it suffices to show that this algorithm simulates the ideal sampling algorithm given in Algorithm 1. Once again we split into the cases of $k$ being odd or even.

In the case of $k$ being odd, Lemma 3.3 gives that a walk $(u_{0},u_{1},\ldots u_{k})$ is sampled with weight at least

[TABLE]

The quality of the distribution produced, follows from Lemma 3.1 and the triangle inequality in Fact 3.2, which then combined with Lemma 2.1 gives the quality of the output. Also, the total size of the sparsifier, as well as the running time follows from Lemma 3.4.

When $k$ is even, by combining Lemmas 3.1 and Lemma 3.5 we have

[TABLE]

Also, note that because $k$ is even, each $k$ step walk in $G$ also corresponds to a walk in $G\times P_{2}$ that starts/ends on the same side, but alternates sides at each step. Using (10) and the symmetry between $u^{(A)}v^{(B)}$ and $u^{(B)}v^{(A)}$, it suffices to sample length $k$ walks with estimated effective resistances satisfying for every edge $(u,v)\in G$

[TABLE]

The rest of the algorithm follows similarly as in the case of odd $k$.

The extra term $\Theta(k\log n)$ in the overhead’s runtime accounts for performing a random walk of length $k$, i.e. after preprocessing in $O(n)$ time an neighboring edge can be sampled using binary search in $O(\log n)$ time. ∎

This reduces the task of efficiently sampling edges in $G^{k}$ to producing good upper bounds for the effective resistances of an explicitly specified graph, either $G$ or $G\times P_{2}$. In the next section we discuss this routine, with focus on density-independent routines.

4 Faster Density Independent Sparsification of Graphs

Note that the monomial sparsification routine only requires a good distribution that dominates the effective resistances. These effective resistances can in turn be computed w.r.t. an approximate graph in a more efficient manner.

Our approach for obtaining density-independent routines follow the approach given in [KLP15]. Namely, we aggressively make the graph more tree-like, and build sparsifiers backwards, each leveraging access to a sparsifier of a graph that is within a constant factor of it.

The main algorithmic difficulty of designing density-independent schemes is that numerically oriented approaches for estimating effective resistances require $O(m\log{n})$ time.

Instead, a more useful method is to utilize low stretch spanning trees to provide an upper bound on all leverage scores in terms of a tree’s stretch. The advantage of this approach is that the stretch of all edges in $G$ w.r.t. a tree can be computed using lowest common ancestor queries in only $O(m)$ time [HT84]. For a particular tree $T$, we define the stretch of an edge $e$ w.r.t. $T$ as the sum of resistances over the unique path $\mathcal{P}_{T}(e)$ in $T$ connecting $e$’s endpoints:

[TABLE]

Extending this definition, the stretch of a graph $G(V,E,\boldsymbol{\mathit{w}})$ w.r.t. $T$ is given by

[TABLE]

Our analysis relies on the following results:

Lemma 4.1**.**

(Lemma 6.4. in **[KLP15]**) If $G$ and $H$ are two graphs such that $\boldsymbol{\mathit{L}}_{G}\preceq\boldsymbol{\mathit{L}}_{H}$, then for any tree $T$, it holds $str_{T}(G)\leq str_{T}(H)$. 2. 2.

If we have a tree $T\preceq G$, then we can construct an $\varepsilon$-sparsifier of $G$ with $O(\varepsilon^{-2}str_{T}(G)\log{n})$ edges in $O(m)$ time. This is a direct consequence of Lemma 2.1, as the stretch upper bounds leverage scores. 3. 3.

(Theorem 1 in **[AN12]**) For any graph $G$, we can obtain a tree with total stretch $\widehat{O}(m\log{n})$ in $\widehat{O}(m\log{n})$ time in the pointer machine model. 4. 4.

(Lemma 5.9 in **[CMP*+*14]**) For any graph $G$ and any parameter $k$, we can find in $\widehat{O}(m)$ time under the RAM model a tree $T$ and a graph $\widehat{G}$ obtained by removing $O(m/k)$ edges such that $str_{T}(\widehat{G})\leq\widehat{O}(m\log{n})$.

We will use these tools to generate a sequence of graphs based on a single low-stretch subgraph. If we set $k\leftarrow\log^{\Theta(1)}n$ in Part 4, the number of edges omitted in $\widehat{G}$ can be sparsified in $O(m)$ time using any of the sparsification methods [KLP*+*16, KPPS17]. This leads to a scheme that start with $\widehat{G}^{(0)}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\widehat{G}$, and creates a chain of graphs where the scaling factor of the tree increases:

[TABLE]

This sequence quickly leads to a graph whose stretch is small enough that we can obtain an $O(1)$-sparsifier in $O(m)$ time.

Lemma 4.2**.**

For any $i\geq\Omega(\log\log{n})$, an $O(1)$-sparsifier of $\widehat{G}^{(i)}$ with $O(n\log{n})$ edges can be found in $O(m)$ time.

Proof.

Multiplying $T$ by a factor of $2^{i}$ reduces the total stretch by the same factor. From Lemma 4.1 Part 4, we get that the total stretch is bounded by $\widehat{O}(2^{-i}m\log{n})$. This means by Lemma 4.1 Part 2, $\widehat{G}^{(i)}$ has an $O(1)$ sparsifier with $O(m\log^{-\Theta(1)}n)$ edges. Invoking a nearly-linear time graph sparsification algorithm on this graph then gives the result. ∎

We will leverage sparsifiers of $\widehat{G}^{(i+1)}$ to construct iteratively, sparsifiers of $\widehat{G}^{(i)}$ using the subroutine shown in Algorithm 3. We also state its guarantees formally below.

Lemma 4.3**.**

Given a $\kappa$-sparsifier $G^{\prime}$ of $G$ and $\varepsilon>0$, $\mathrm{TreeSparsify}(G,G^{\prime},\kappa,\varepsilon)$ produces an $\varepsilon$-sparsifier of $G$ with at most $\widetilde{O}(\varepsilon^{-2}\log^{2}{n}\cdot\kappa\left|E(G^{\prime})\right|)$ edges in $\widetilde{O}(m+\varepsilon^{-2}\log^{3}{n}\cdot\kappa\left|E(G^{\prime})\right|)$ time.

Proof.

To apply Lemma 2.1, we have to compute a vector $\boldsymbol{\widetilde{\mathit{r}}}\geq\mathcal{R}_{\mathrm{eff}}^{G}$ and give an upper bound on $\langle\boldsymbol{\mathit{w}},\boldsymbol{\widetilde{\mathit{r}}}\rangle$. Since $\boldsymbol{\mathit{L}}_{T}\preceq\boldsymbol{\mathit{L}}_{G^{\prime}}\preceq\kappa\boldsymbol{\mathit{L}}_{G}$ it follows that

[TABLE]

Moreover, by combining (12), $\boldsymbol{\mathit{L}}_{G}\preceq\boldsymbol{\mathit{L}}_{G^{\prime}}$, Lemma 4.1, Part 1 and Part 3 we obtain

[TABLE]

This yields the overall edge count and runtime. ∎

We present two density-independent sparsification algorithms that iteratively construct sparsifiers backwards from $\widehat{G}^{(k)}$ by:

Creating a crude $\epsilon/2$-sparsifier of $\widehat{G}^{(i)}$, $G^{\prime(i)}$ with $\widehat{O}(\epsilon^{-2}n\log^{3}{n})$ edges using TreeSparsify with $\widetilde{G}^{(i+1)}$, the sparsifier of $\widehat{G}^{(i+1)}$ constructed in the previous step as guide. 2. 2.

Further sparsify this crude sparsifier, $G^{\prime(i)}$, with error $\epsilon/2$ to from $\widetilde{G}^{(i)}$, which has the desired edge count, and will be used in the next step.

We refrain from providing pseudocode of these steps because of the subtle differences in the resulting algorithms.

The current fastest sparsification routines compute effective resistances via the Johnson-Lindenstrauss transform [SS11], which in turn requires the use of fast linear system solvers [KLP*+*16].

Lemma 4.4**.**

Given a graph $G$, we can compute $2$-approximations to its effective resistances in $\widehat{O}(m\log{n}+n\log^{2}n)$ time.

This runtime bound can be obtained by letting the depth approach $n$ in the proof of Theorem 1.2 on page 49 of [KLP*+*16]. The effective resistances can in turn be extracted from the call to Sparsify made at $i=0$ in the pseudocode in Figure 11 on page 46. We omit details on these steps in the hope that significantly simpler sparsification routines with similar performances will be developed.

Combining this with the sequence of graphs defined in (11) gives:

Corollary 4.5**.**

There is a routine that takes a weighted undirected graph $G$ with $n$ vertices, $m$ edges, an error $\epsilon>0$, and produces in $\widehat{O}(m+\varepsilon^{-2}\log^{4}{n})$ time an $\epsilon$-sparsifier of $G$ with $O(\epsilon^{-2}n\log{n})$ edges, as well as leverage score upper bounds that sum up to $\widehat{O}(n\log^{2}n)$.

Proof.

Consider the sequence of matrices as defined in (11). The sparsifier for $\widehat{G}^{O(\log\log{n})}$ is given by Lemma 4.2.

Then we can iteratively build $O(n\log{n})$ sized $O(1)$-sparsifiers for $\widehat{G}^{(i)}$ all the way up to $i=1$. The cost of invoking $\mathrm{TreeSparisfy}(\widehat{G}^{(i)},\widetilde{G}^{(i+1)},O(1),2)$ at each step is $\widehat{O}(m+n\log^{4}n)$, while the resulting sparsifier has $\widehat{O}(n\log^{3}n)$ edges. Lemma 4.4 then turns this into an $O(1)$-sparsifier for $\widehat{G}^{(i)}$.

At the last step of $i=0$, we invoke these same routines, but now with error $\epsilon$ to obtain the $\epsilon$-sparsifier. Note that the effective resistance upper bounds are computed during the call to $\mathrm{TreeSparisfy}$. ∎

The guarantees of this routine fits into the requirements of the random walk sampling algorithm from Section 3 and yields the faster, density-independent algorithm for sparsifying $G^{k}$, that is our main result.

Proof.

(Of Theorem 1.1) The leverage scores upper bound obtained by Corollary 4.5, when combined with Theorem 3.6 produces an $\epsilon$-sparsifier for $G^{k}$ with $\widehat{O}(\epsilon^{-2}kn\log^{3}n)$ edges in $\widehat{O}(m+\varepsilon^{-2}k^{2}\log^{4}{n})$ time. Sparsifying this graph once again using Lemma 4.4 then leads to the main result as stated in Theorem 1.1. ∎

There are also purely combinatorial constructions of graph sparsifiers based on spanners [KP12, Kou14, KPPS17]:

Lemma 4.6**.**

(Theorem 4.1. in [KPPS17]) Given a graph $G$ and an error $\varepsilon>0$, we can compute an $\varepsilon$-spectral sparsifier of $G$ with $\widehat{O}(n\log^{2}{n})$ edges in $\widehat{O}(m\log^{2}{n})$ time.

The algorithm in Lemma 4.6, applied to our sparsification scheme gives

Corollary 4.7**.**

There is a combinatorial algorithm that for any graph $G$ on $n$ vertices and $m$ edges, and any error $\epsilon>0$, produces in $\widehat{O}(m+n\log^{6}n)$ time an $\epsilon$-sparsifier with $\widehat{O}(\varepsilon^{-2}n\log^{2}{n})$ edges, as well as leverage score upper bounds that sum up to $\widehat{O}(n\log^{3}{n})$.

Proof.

This is similar to the routine calling numerical sparsifiers outlined in Corollary 4.5. However, each of the $\widetilde{G}^{(i+1)}$ now has $\widehat{O}(n\log^{2}{n})$ edges. Hence, the first crude approximation $\widetilde{G^{\prime}}^{(i)}$ has $\widehat{O}(n\log^{4}{n})$ edges. Sparsifying it down to $\widehat{O}(n\log^{2}{n})$ edges takes $\widehat{O}(n\log^{6}{n})$ time. ∎

Appendix A Deferred Proofs

We provide now some additional details on Lemmas from Section 3 that are direct consequences of steps in previous works. The total summation of the sampling weights follows from a summation identical to the special case of uniform sampling, as presented in [CCL*+*15, Lemma 29]. More precisely, by evaluating the total weights of all random walks that involve a particular edge $e\in G$.

See 3.4

Proof.

We first show by induction that the total weights of all length $k$ walks whose $i$th edge is $e$ is exactly $\boldsymbol{\mathit{w}}_{e}$.

The base case of $k=1$ is trivial as only $e$ is a length $1$ walk between $u_{0}$ and $u_{1}$.

The inductive case of $k>1$ has two cases: $i>0$ or $i<k-1$. We consider the $i>0$ case only, as the other one follows by symmetry. Expanding the weight of a length $k$ walk gives:

[TABLE]

The fact that $i<k-1$ means that $u_{k}$ can be any neighbor of $u_{k-1}$, leading to a sum that cancels the $\boldsymbol{\mathit{d}}_{u_{k-1}}$ term in the denominator. Formally:

[TABLE]

The result then follows from the inductive hypothesis applied to walks of length $k-1$ that have edge $i$ as $e$.

The proof then uses a double counting argument that breaks the summation over the edge $(u_{i},u_{i+1})$, and by noting that the choice over $i$ implies that each edge is picked exactly $k$ times. We can rewrite the original summation as:

[TABLE]

The proof above gives that the term within the bracket is $\boldsymbol{\mathit{w}}_{e}$. So the summation over $i$ is just an extra factor of $k$, by which we obtain the result.

∎

The equivalence of effective resistances in $G^{2}$ and $G\times P_{2}$ requires the definition of Schur complements.

Definition A.1** (Schur Complement).**

The Schur Complement of a symmetric matrix in block form: $\boldsymbol{\mathit{M}}=\begin{pmatrix}\boldsymbol{\mathit{M}}_{[F,F]}&\boldsymbol{\mathit{M}}_{[F,C]}\\ \boldsymbol{\mathit{M}}_{[C,F]}&\boldsymbol{\mathit{M}}_{[C,C]}\end{pmatrix}$ that removes the block $F$ is:

[TABLE]

It can be seen that, if $\boldsymbol{\mathit{M}}$ was the Laplacian of a graph $G$, $Sc(M,F)$ is the Laplacian of the graph $G^{S}$, which can be formed using the following iterative process:

• •

Iteratively for all vertices $u\in F$ (The set of vertices represented by the columns in $F$ of $\boldsymbol{\mathit{M}}$.)

• –

For all pairs of edges $uv_{1}$ and $uv_{2}$ in the current graph (including edges added in previous steps), delete them and add the edge $v_{1}v_{2}$ with weight

[TABLE]

where $\boldsymbol{\mathit{d}}_{u}$ is the weighted degree of $u$ (once again w.r.t. the current graph).

• –

Delete $v$

Lemma A.2**.**

For every vector \boldsymbol{\mathit{z}}=\left(\begin{array}[]{c}\boldsymbol{\mathit{z}}_{1}\\ 0\end{array}\right) it holds that

[TABLE]

By symmetry for any vector \boldsymbol{\mathit{z}}=\left(\begin{array}[]{c}0\\ \boldsymbol{\mathit{z}}_{2}\end{array}\right) it holds that

[TABLE]

In particular, the effective resistances are maintained under Schur complement.

Proof.

Consider the linear system

[TABLE]

Since $\boldsymbol{\mathit{z}}_{2}=0$, we have

[TABLE]

and thus

[TABLE]

∎

See 3.5

Proof.

Notice that

[TABLE]

is the Schur Complement of

[TABLE]

with respect to one half of the vertices, e.g. $V^{(B)}$. So, the lemma follows from Lemma A.2. ∎