Fully Dynamic Spectral Vertex Sparsifiers and Applications

TL;DR

This paper introduces dynamic algorithms for maintaining spectral vertex sparsifiers in graphs, enabling efficient updates and queries for electrical resistances, Laplacian systems, and energy computations, with applications to graph analysis.

Contribution

It presents the first sublinear-time dynamic data structures for spectral sparsifiers and related primitives, without assumptions on graph topology.

Findings

Supports insertions and deletions of edges and terminals in sublinear time.

Provides approximate solutions to Laplacian systems with efficient update/query times.

Maintains all-pairs effective resistance with high probability in dynamic graphs.

Abstract

We study \emph{dynamic} algorithms for maintaining spectral vertex sparsifiers of graphs with respect to a set of terminals $T$ of our choice. Such objects preserve pairwise resistances, solutions to systems of linear equations, and energy of electrical flows between the terminals in $T$. We give a data structure that supports insertions and deletions of edges, and terminal additions, all in sublinear time. Our result is then applied to the following problems. (1) A data structure for maintaining solutions to Laplacian systems $\mathbf{L} \mathbf{x} = \mathbf{b}$, where $\mathbf{L}$ is the Laplacian matrix and $\mathbf{b}$ is a demand vector. For a bounded degree, unweighted graph, we support modifications to both $\mathbf{L}$ and $\mathbf{b}$ while providing access to $\epsilon$-approximations to the energy of routing an electrical flow with demand $\mathbf{b}$, as well as query access to entries of a vector $\tilde{\mathbf{x}}$ such that $\left\lVert \tilde{\mathbf{x}}-\mathbf{L}^{\dagger} \mathbf{b} \right\rVert_{\mathbf{L}} \leq \epsilon \left\lVert \mathbf{L}^{\dagger} \mathbf{b} \right\rVert_{\mathbf{L}}$ in $\tilde{O}(n^{11/12}\epsilon^{-5})$ expected amortized update and query time. (2) A data structure for maintaining All-Pairs Effective Resistance. For an intermixed sequence of edge insertions, deletions, and resistance queries, our data structure returns $(1 \pm \epsilon)$-approximation to all the resistance queries against an oblivious adversary with high probability. Its expected amortized update and query times are $\tilde{O}(\min(m^{3/4},n^{5/6} \epsilon^{-2}) \epsilon^{-4})$ on an unweighted graph, and $\tilde{O}(n^{5/6}\epsilon^{-6})$ on weighted graphs. These results represent the first data structures for maintaining key primitives from the Laplacian paradigm for graph algorithms in sublinear time without assumptions on the underlying graph topologies.

Equations

$$R_{\text{eff}}^G\left(u,v\right):={\mathbf{\chi }}_{u,v}^{\top }{\mathbf{L}}_G^{†}{\mathbf{\chi }}_{u,v}.$$

$$(1-ϵ)R_{\text{eff}}^G(u,v)\le \varphi \le (1+ϵ)R_{\text{eff}}^G(u,v),$$

$$\mathbf{L}:=\left[\begin{array}[]{cc}\mathbf{L}_{\left[F,F\right]}&\mathbf{L}_{\left[F,T\right]}\\ \mathbf{L}_{\left[T,F\right]}&\mathbf{L}_{\left[T,T\right]}\end{array}\right].$$

$$\mathbf{SC}\left(G,T\right)={\mathbf{L}}_{\left[T,T\right]}-{\mathbf{L}}_{\left[T,F\right]}{\mathbf{L}}_{\left[F,F\right]}^{-1}{\mathbf{L}}_{\left[F,T\right]}.$$

$$R_{\text{eff}}^G\left(u,v\right)=R_{\text{eff}}^{\text{textsc}SC\left(G,T\right)}\left(u,v\right).$$

$$(1-\epsilon ){\mathbf{x}}^{\top }{\mathbf{L}}_G\mathbf{x}\le {\mathbf{x}}^{\top }{\mathbf{L}}_H\mathbf{x}\le (1+\epsilon ){\mathbf{x}}^{\top }{\mathbf{L}}_G\mathbf{x}.$$

$$(1-\epsilon )R_{\text{eff}}^G(u,v)\le R_{\text{eff}}^H(u,v)\le (1+\epsilon )R_{\text{eff}}^G(u,v).$$

$$\mathbf{P}({T}):=\left[\begin{array}[]{cc}-\mathbf{L}_{\left[T,F\right]}\mathbf{L}_{\left[F,F\right]}^{-1}&\mathbf{I}_{T}\end{array}\right].$$

$$\left[\mathbf{P}(T)1_u\right](v)={\mathbb{P}}_u\left[t_v<t_{T\setminus v}\right].$$

$${\Vert 1_u-\mathbf{P}(T)1_u\Vert }_{{\mathbf{L}}^{†}}^2=(1_u-\mathbf{P}(T)1_u{)}^{\top }{\mathbf{L}}^{†}(1_u-\mathbf{P}(T)1_u).$$

$$\rho =O\left(\log n{ϵ}^{-2}\right).$$

$$1/\left(\rho \sum _{i=0}^{\ell -1}\left(1/w_{u_i{u}_{i+1}}\right)\right)$$

$$\tilde{O}({\beta }^{-3}+\beta m),$$

$$\ell =c_{BF}\cdot {\delta }^3\cdot {\beta }^{-2}{\log }^{3}n$$

$$\delta \cdot {\beta }^{-1}\cdot \log n$$

$$\begin{array}{rl}{\mathbb{P}}_{T,w}\left[V\left(w\right)\cap T=\varnothing \right]& ={\mathbb{P}}_{w,T}\left[V\left(w\right)\cap T=\varnothing \right]\\ & \le {\mathbb{P}}_w\left[\left|E\left(w\right)\right|\le z\right]\\ & \text{leavevmode}\text{leavevmode}+{\mathbb{P}}_{w:|E\left(w\right)|\ge z,T}\left[V\left(w\right)\cap T=\varnothing \right].\end{array}$$

$$z=\delta \cdot {\beta }^{-1}\cdot \log n,$$

$$\left|E\left(\hat{w}\right)\right|\ge \delta \cdot {\beta }^{-1}\log n.$$

$${\left(1-\beta \right)}^{|E(\hat{w})|}\le {\left(1-\beta \right)}^{\delta \cdot {\beta }^{-1}\log n}\le e^{-\delta \log n}\le n^{-\delta },$$

$$\frac{\tilde{O}(m{\beta }^{-2}{ϵ}^{-4})}{\beta m}=\tilde{O}({\beta }^{-3}{ϵ}^{-4}).$$

$$\tilde{O}\left((\beta m+{\beta }^{-3}){ϵ}^{-4}\right),$$

$$\tilde{O}(m^{3/4}{ϵ}^{-4}).\text{qed}$$

$$s(w):=\sum _{i=1}^t\frac{1}{w_{w_{i-1}{w}_i}},$$

$$p_i(v)$$

$$p^{\mathrm{new}}(i)$$

$$\left[\begin{array}[]{l}\boldsymbol{\mathit{p}}_{i}\\ \mathit{p}^{\mathrm{new}}(i)\end{array}\right]=\boldsymbol{\mathit{W}}\cdot\left[\begin{array}[]{l}\boldsymbol{\mathit{p}}_{i-1}\\ \mathit{p}^{\mathrm{new}}(i-1)\end{array}\right],\quad\forall i\geq 1.$$

$${\hat{p}}_i(v)=p_i(v),\text{leavevmode}\forall v\in U\text{and}{\hat{p}}_i(k+1)=p^{\mathrm{new}}(i).$$

$$(rp-\ell p)={\mathbb{P}}_{X(u,U)}\left[\ell \le X(u,U)\le r\right].$$

$${\mathbb{P}}_{X(u,U)}\left[X(u,U)\le \eta \mid \ell \le X(u,U)\le r\right]$$

$$=\frac{{\mathbb{P}}_{X(u,U)}\left[\ell \le X(u,U)\le \eta \right]}{{\mathbb{P}}_{X(u,U)}\left[\ell \le X(u,U)\le r\right]}$$

$$=\frac{(p^{\mathrm{new}}(\eta )-\ell p)}{(rp-\ell p)}.$$

$$r(v,z)=\frac{q(v)w_{vz}}{R},\text{ where }R:=\sum _{v\in U,z\in V\setminus U}q(v)w_{vz}.$$

Full text

(1) A data structure for dynamically maintaining solutions to Laplacian systems $\mathbf{L}\mathbf{x}=\mathbf{b}$, where $\mathbf{L}$ is the graph Laplacian matrix and $\mathbf{b}$ is a demand vector. For a bounded degree, unweighted graph, we support modifications to both $\mathbf{L}$ and $\mathbf{b}$ while providing access to $\epsilon$-approximations to the energy of routing an electrical flow with demand $\mathbf{b}$, as well as query access to entries of a vector $\tilde{\mathbf{x}}$ such that $\left\lVert\tilde{\mathbf{x}}-\mathbf{L}^{\dagger}\mathbf{b}\right\rVert_{\mathbf{L}}\leq\epsilon\left\lVert\mathbf{L}^{\dagger}\mathbf{b}\right\rVert_{\mathbf{L}}$ in $\tilde{O}(n^{11/12}\epsilon^{-5})$ expected amortized update and query time.

(2) A data structure for maintaining fully dynamic All-Pairs Effective Resistance. For an intermixed sequence of edge insertions, deletions, and resistance queries, our data structure returns $(1\pm\epsilon)$-approximation to all the resistance queries against an oblivious adversary with high probability. Its expected amortized update and query times are $\tilde{O}(\min(m^{3/4},n^{5/6}\epsilon^{-2})\epsilon^{-4})$ on an unweighted graph, and $\tilde{O}(n^{5/6}\epsilon^{-6})$ on weighted graphs.

The key ingredients in these results are (1) the intepretation of Schur complement as a sum of random walks, and (2) a suitable choice of terminals based on the behavior of these random walks to make sure that the majority of walks are local, even when the graph itself is highly connected and (3) maintenance of these local walks and numerical solutions using data structures.

These results together represent the first data structures for maintaining key primitives from the Laplacian paradigm for graph algorithms in sublinear time without assumptions on the underlying graph topologies. The importance of routines such as effective resistance, electrical flows, and Laplacian solvers in the static setting make us optimistic that some of our components can provide new building blocks for dynamic graph algorithms.

1 Introduction

Problems arising from analyzing and understanding graph structures have motivated the development of many powerful tools for storing and compressing graphs and networks. One such tool that has received a considerable amount of attention over the past two decades is graph sparsification [10, 9]. Roughly speaking, a graph sparsifier is a “compressed” version of a large input graph that preserves important properties like distance information [56], cut value [10] or graph spectrum [65]. Graph Sparsifiers fall into two main categories: edge sparsifiers, which are graphs that reduce the number of edges, and vertex sparsifiers, which are graphs that reduce the number of vertices. Both categories have many applications in approximation algorithms [22, 59], machine learning [46, 45, 72], and most recently efficient graph algorithms [66, 47, 63, 36]. While edge sparsifiers have played an instrumental role in obtaining nearly linear time algorithms [9], their practical applicability is somewhat limited due to the fact most of the large networks are already sparse. On the other hand, vertex sparsifiers address the “real” compression of large networks by reducing the number of vertices.

While vertex sparisifers in general are significantly more difficult to generate [51, 13, 49], a notable exception is vertex sparsifiers for quadratic minimization problems, otherwise known as Schur complements. Concretely, given an undirected, weighted graph $G$, a subset of terminal vertices $T$ and its corresponding Laplacian matrix, a graph $H$ with $V(H)=T$ is a vertex resistance sparsifier of $G$ with respect to $T$ if the Laplacian matrix of $H$ is obtained by the Schur complement of the Laplacian of $G$ with respect to $T$. Schur complement is a central concept in physics and linear algebra with a wide range of applications including multi-grid solvers, Markov chains and finite-element analysis [16], and has also recently found extensive applications in graph algorithms [40, 41, 18, 19, 62, 61].

Most of the massive graphs in the real world, such as social networks, the web graph, are subject to frequent changes over time. This dynamic behavior of graphs has been studied for several important graph problems, where the basic idea is to maintain problem solutions as graphs undergo edge insertions and deletions in time faster than recomputing the solution from scratch. Dynamic graph algorithms have also been formulated for many problems that involve edge sparsifiers [32, 70, 35, 26], as well important variants of edge sparsifiers themselves, including minimum spanning trees [33, 52, 73, 53], spanners [8], spectral sparsifiers [2], and low-stretch spanning trees [27]. However, despite the increasing importance of high quality vertex sparsifiers in graph algorithms, very little is known about their maintenance in the dynamic setting.

In this paper we give the first non-trival dynamic algorithms for maintaining Schur complements of general graphs with respect to a set of terminal of our choice. Our data-structure maintains at any point of time a $(1\pm\epsilon)$ approximation to the Schur complement while supporting insertions and deletions of edges, and arbitrary vertex additions to the terminal set. To the best of our knowledge, prior dynamic Schur complement algorithms were only known for minor-free graphs [24, 25].

Lemma 1.1**.**

Given an error parameter $\epsilon>0$, an unweighted undirected multi-graph $G=(V,E)$ with $n$ vertices, $m$ edges, a subset of terminal vertices $T^{\prime}$ and a parameter $\beta\in(0,1)$ such that $|T^{\prime}|=O(\beta m)$, there is a data-structure DynamicSC$(G,T^{\prime},\beta)$ for maintaining a graph $\tilde{H}$ with $\mathbf{L}_{\tilde{H}}\approx_{\epsilon}\mathbf{SC}(G,T)$ for some $T$ with $T^{\prime}\subseteq T$, $|T|=O(\beta m)$, while supporting $O(\beta m)$ operations in the following running times:

• •

Initialize$(G,T^{\prime},\beta)$: Initialize the data-structure, in $\tilde{O}(m\beta^{-2}\epsilon^{-4})$555We use $\tilde{O}(f)$ to denote $O\left(f\cdot\mathrm{poly}\left(\log n\right)\right)$ for a function $f$. expected amortized time.

• •

Insert$(u,v)$: Insert the edge $(u,v)$ to $G$ in $\tilde{O}(1)$ amortized time.

• •

Delete$(u,v)$: Delete the existing edge $(u,v)$ from $G$ in $\tilde{O}(1)$ amortized time.

• •

AddTerminal$(u)$*: Add $u$ to $T^{\prime}$ in $\tilde{O}(1)$ amortized time.

Our guarantees hold against an oblivious adversary.

Our algorithm extends to weighted graphs, albeit with slightly larger running time guarantees. Concretely we give an algorithm that maintains an approximate Schur Complement with $\tilde{O}(m\beta^{-4}\epsilon^{-4})$ expected amortized time for initializing the data-structure, and $\tilde{O}(1)$ amortized time for the remaining operations.

The key algorithmic components behind the result in unweighted graphs are (1) the interpretation of Schur complement as a sum of random walks and (2) randomly picking a terminal vertex subset onto which the vertex resistance sparsifiers is constructed. Specifically, in a novel way we combine random walk based methods for generating resistance vertex sparsifiers [19] with results in combinatorics that bound the speed at which such walks spread among distinct edges [6]. Our result in the weighted case essentially follows the same idea except that the speed at which random walks visit distinct edges in weighted networks could be very slow. To control this, we show how to efficiently sample a distinct edge without the need to simulate every step of the walk. This interacts well with other parts of our data structure and leads to comparable running time guarantees.

We show the applicability of our dynamic Schur complement to two cornerstone problems in graph Laplacian literature, namely dynamic Laplacian solver [66] and dynamic All-Pairs Effective Resistances [67].

Solving linear systems lies at the heart of many problems arising in scientific computing, numerical linear algebra, optimization and computer science. An important subclass of linear systems are Laplacian systems, which arise in many natural contexts, including computation of voltages and currents in electrical network. Solving Laplacian system has received increasing attention over the past years after the breakthrough work of Spielman and Teng [66] who gave the first near-linear time algorithm. Motivated by fast Laplacian solvers in different model of computations [5, 58], we initiate the study of algorithms for dynamically solving Laplacian systems. Concretely, given a graph Laplacian $\mathbf{L}\in\mathbb{R}^{n\times n}$ and a vector $\mathbf{b}\in\mathbb{R}^{n}$ in the range of $\mathbf{L}$, the goal is to maintain an $\mathbf{x}$ such that $\mathbf{L}\mathbf{x}=\mathbf{b}$, while off-diagonals of $\mathbf{L}$ and the entries of $\mathbf{b}$ change over time. To allow for sub-linear query times, here we focus on querying one (or a few) coordinates of $\mathbf{x}$. Formally, given any index $u\in\{1,\ldots,n\}$, the goal is to output $\tilde{\mathbf{x}}(u)$ for some approximation $\tilde{\mathbf{x}}$ of $\mathbf{L}^{\dagger}\mathbf{b}$. Our contribution is the first sub-linear dynamic Laplacian solver in bounded degree graphs.

Theorem 1.2**.**

For any given error threshold $m^{-1}<\epsilon<1$, there is a data-structure for maintaining an unweighted, undirected bounded degree multi-graph $G=(V,E)$ with $n$ vertices, $m$ edges and a vector $\mathbf{b}\in\mathbb{R}^{n}$ that supports the following operations in $\tilde{O}(n^{11/12}\epsilon^{-5})$ expected amortized time:

• •

Insert$(u,v)$: Insert the edge $(u,v)$ with resistance $1$ in $G$.

• •

Delete$(u,v)$: Delete the edge $(u,v)$ from $G$.

• •

Change$(u,\mathbf{b}^{\prime}(u),v,\mathbf{b}^{\prime}(v))$: Change $\mathbf{b}(u)$ to $\mathbf{b}^{\prime}(u)$ and $\mathbf{b}(v)$ to $\mathbf{b}^{\prime}(v)$ while keeping $\mathbf{b}$ in the range of $\mathbf{L}$.

• •

Solve$(u)$: Return $\tilde{\mathbf{x}}(u)$ with $\tilde{\mathbf{x}}$ such that $\left\lVert\tilde{\mathbf{x}}-\mathbf{L}^{{\dagger}}\mathbf{b}\right\rVert_{\mathbf{L}}\leq\epsilon\left\lVert\mathbf{L}^{{\dagger}}\mathbf{b}\right\rVert_{\mathbf{L}}.$

Our guarantees hold against an oblivious adversary.

Note that the $\tilde{\mathbf{x}}$ in the theorem above is not guaranteed to be inside the range of $\mathbf{L}_{G}$ and it only preserves the differences between vertices in the same connected component.

We observe that conditioning on the vector $\mathbf{b}$ having small support, i.e., a small number of non-zero elements, immediately leads to a dynamic solver by just including the corresponding vertices into the Schur complement, and maintaining a dynamic Schur complement onto these vertices augmented with some carefully chosen additional terminals. Upon receipt of a query index, we add the corresponding vertex to the current Schur complement and simply solve a linear system there. However, note that the demand vector may have a large number of non-zero entries, thus preventing us from obtaining a sub-linear time algorithm with this approach. We alleviate this by projecting this vector onto the set of current terminals and showing that such projection can be maintained dynamically while introducing controllable error in the approximation guarantee.

A consequence of dynamic Laplacian solver is that we can maintain the energy of the electrical flow when routing any vector $\mathbf{b}$ in the range of $\mathbf{L}$, which is a generalization of the All-Pairs Effective Resistance problem with $\mathbf{b}$ having exactly two non-zero entries, one of them being $1$ and the other $-1$.

Another application of our technique is dynamic maintainance of effective resistance, a well studied quantity that has direct applications in random walks, spanning trees [48] and graph sparsification [67]. We maintain (approximate) All-Pairs Effective Resistances of a graph $G$ among any pair of query vertices while supporting an intermixed sequence of edge insertions and deletions in $G$. Our study is also motivated in part by the wide usage of commute distances, a random walk-based similarity measure that has been successfully employed in important practical applications such as link predictions [44]. Since commute distance is a scaled version of effective resistance, our dynamic algorithm readily extends to this graph measure while achieving the same approximation and running time guarantees.

Theorem 1.3**.**

For any given error threshold $\epsilon>0$, there is a data-structure for maintaining an unweighted, undirected multi-graph $G=(V,E)$ with up to $m$ edges that supports the following operations in $\tilde{O}(m^{3/4}\epsilon^{-4})$ expected amortized time:

• •

Insert$(u,v)$: Insert the edge $(u,v)$ with resistance $1$ in $G$.

• •

Delete$(u,v)$: Delete the edge $(u,v)$ from $G$.

• •

EffectiveResistance$(s,t)$: Return a $(1\pm\epsilon)$-approximation to the effective resistance between $s$ and $t$ in the current graph $G$.

Our guarantees hold against an oblivious adversary.

Our algorithm can also handle weighted graphs, albeit with a bound of $\tilde{O}(m^{5/6}\epsilon^{-4})$ on the expected amortized update and query time. By running this algorithm on the output of a dynamic spectral sparsifier [2], we obtain a bound of $\tilde{O}(n^{5/6}\epsilon^{-6})$ per operation, which is truly sub-linear irrespective of graph density. We discuss such improvements in Sections 4.3.

We are optimistic that our algorithmic ideas could be useful for dynamically maintaining a wider range of graph properties. Both the results that we give dynamic algorithms for, vertex sparsifiers and Schur complements, have wide ranges of applications in static settings, with the latter being at the core of the ‘Laplacian paradigm’ of graph algorithms [68, 69]. While it’s less clear that solutions across multiple Laplacian solves can be propagated to each other as the input dynamically changes, repeated sparsification on the other hand represents a routine that composes and interacts well with a much wider range of primitives. As a result, we are optimistic that it can be used as a building block in dynamic versions of many existing applications of Laplacian solvers.

1.1 Related Works

The recent data structures for maintaining effective resistances in planar graphs [24, 25] drew direct connections between Schur complements and data structures for maintaining them in dynamic graphs. This connection is due to the preservation of effective resistances under vertex eliminations (Fact 2.2). From this perspective, the Schur complement can be viewed as a vertex sparsifier for preserving resistances among a set of terminal vertices.

The power of vertex or edge graph sparsifiers, which preserve certain properties while reducing problem sizes, has long been studied in data structures [20, 21]. Ideas from these results are central to recent works on offline maintenance for $3$-connectivity [57], generating random spanning trees [18], and new notions of centrality for networks [43]. Our result is the first to maintain such vertex sparsifiers, specifically Schur complements, for general graphs in online settings.

While the ultimate goal is to dynamically maintain (approximate) minimum cuts and maximum flows, effective resistances represent a natural ‘first candidate’ for this direction of work due to them having perfect vertex sparsifiers. That is, for any subset of terminals, there is a sparse graph on them that approximately preserves the effective resistances among all pairs of terminals. This is in contrast to distances, where it’s not known whether such a graph can be made sparse, or in contrast to cuts, where the existence of such a dense graph is not known (assuming that we are not content with large constant or poly-logarithmic approximations).

Dynamic Graph Algorithms.

The maintenance of graph properties in dynamic algorithms has been a major area of ongoing research in data structures. The problems being maintained include $2-$ or $3-$connectivity [21, 32, 33], shortest paths [29, 30, 11, 1, 15], global minimum cut [31, 70, 42, 26], maximum matching [55, 28, 12], and maximal matching [7, 54, 64]. Perhaps most closely related to our work are dynamic algorithms that maintain properties related to paths [23, 21, 32, 35, 73, 53, 52]. In particular, the work of Wulff-Nilsen [73] also utilizes the behavior of random walks under edge deletions to keep track of low-conductance cuts.

Dynamic algorithms for evaluating algebraic functions such as matrix determinant and matrix inverse has also been considered [60]. One application of such algorithms is that they can be used to dynamically maintain single-pair effective resistance. Specifically, using the dynamic matrix inversion algorithm, one can dynamically maintain exact $(s,t)$-effective resistance in $O(n^{1.575})$ update time and $O(n^{0.575})$ query time.

Vertex Sparsifiers.

Vertex sparsifiers have been studied in more general settings for preserving cuts and flows among terminal vertices [51, 13, 39]. Efficient versions of such routines have direct applications in data structures, even when they only work in restricted settings: terminal sparsifiers on quasi-bipartite graphs [4] were core routines in the data structure for maintaining flows in bipartite undirected graphs [2].

Our data structure utilizes vertex sparsifiers, but in even more limited settings as we get to control the set of vertices to sparsify onto. Specifically, the local maintenance of this sparsifier under insertions and deletions hinges upon the choice of a random subset of terminals, while vertex sparsifiers usually need to work for any subset of terminals. Evidence from numerical algorithms [40, 19] suggest this choice can significantly simplify interactions between algorithmic components. We hope this flexibility can motivate further studies of vertex sparsifiers in more restrictive, but still algorithmically useful settings.

Organization.

The paper is organized as follows. We discuss preliminaries in Section 2 and give an overview of the key techniques in Section 3. After that we give a data-structure for dynamic Schur complement on unweighted graphs in Section 4, which can be applied to the dynamic All-Pairs Effective Resistance problem. In Section 4.3, we briefly discuss the extension of our data-structure to weighted graphs. In Section 5, we give a high-level idea behind a data-structure that dynamically maintains the projection of a vector onto a subset of vertices of an unweighted bounded degree graph, which can be combined with dynamic Schur complement to give a dynamic Laplacian solver. In Appendix A, we show that our bound on the maximum load of a vertex from Lemma 5.5 is essentially tight. In Appendix B, we provide details on the graph approximation guarantees and properties of projections that our random walk sampling and other routines rely on. Finally, in Appendix C, we provide an algorithm for approximately sampling the sum of reciprocals of the edge weights of a random walk which allows us to generate long random walks without going through each step.

2 Preliminaries

In our dynamic setting, an undirected, weighted multi-graph undergoes both insertions and deletions of edges. We let $G=(V,E,\boldsymbol{\mathit{w}})$ always refer to the current version of the graph. We will use $n$ and $m$ to denote bounds on the number of vertices and edges at any point, respectively.

For an unweighted, undirected multi-graph $G$, let $\textbf{A}_{G}$ denote its adjacency matrix and let $\textbf{D}_{G}$ its degree diagonal matrix (counting edge multiplicities for both matrices). The graph Laplacian $\mathbf{L}_{G}$ of $G$ is then defined as $\mathbf{L}_{G}=\textbf{D}_{G}-\textbf{A}_{G}$. Let $\mathbf{L}_{G}^{{\dagger}}$ denote the Moore-Penrose pseudo-inverse of $\mathbf{L}_{G}$. Subscripts will be often omitted when the underlying graph is clear from the context. We define the indicator vector $\boldsymbol{1}_{u}\in\mathbb{R}^{n}$ of a vertex $u$ such that $\boldsymbol{1}_{u}(v)=1$ if $v=u$, and $\boldsymbol{1}_{u}(v)=0$ otherwise. Let $\boldsymbol{\mathit{d}}(u)=\sum_{v:(u,v)\in E}\boldsymbol{\mathit{w}}_{u,v}$ be the weighted degree of a vertex $u$.

A walk in $G$ is a sequence of vertices such that consecutive vertices are connected by edges. A random walk in $G$ is a walk that starts at a starting vertex $w_{0}$, and at step $i\geq 1$, the vertex $w_{i}$ is chosen randomly among the neighbors of $w_{i-1}$. If graph $G$ is unweighted, then each of its neighbors becomes $w_{i}$ with equal probability. If $G$ is weighted, the probability $\mathbb{P}_{w}\left[w_{i}=u\mid w_{0},\ldots,w_{i-1}\right]$ is proportional to the edge weight $\boldsymbol{\mathit{w}}_{w_{i-1}u}$.

Effective Resistance.

For our algorithm, it will be useful to define effective resistance using linear algebraic structures. Specifically, given any two vertices $u$ and $v$ in $G$, if $\boldsymbol{\chi}(u,v):=\boldsymbol{1}_{u}-\boldsymbol{1}_{v}$, then the effective resistance between $u$ and $v$ is given by

[TABLE]

Linear systems in graph Laplacian matrices can be solved in nearly-linear time [38]. One prominent application of these solvers is the approximation of effective resistances.

Lemma 2.1**.**

Fix $\epsilon\in(0,1)$ and let $G=(V,E)$ be any graph with two arbitrary distinguished vertices $u$ and $v$. There is an algorithm that computes a value $\phi$ such that

[TABLE]

in $\tilde{O}(m+n/\epsilon^{2})$ time with high probability.

Schur complement.

Given a graph $G=(V,E)$, we can think of the Schur complement as the partially eliminated state of $G$. This relies on some partitioning of $V$ into two disjoint subset of vertices $T$, called terminals and $F:=V\setminus T$, called non-terminals, which in turn partition the Laplacian $\mathbf{L}$ into $4$ blocks:

[TABLE]

The Schur complement onto $T$, denoted by $\mathbf{SC}(G,T)$, is the matrix obtained after eliminating the variables in $F$. Its closed form is given by

[TABLE]

It is well known that $\mathbf{SC}(G,T)$ is a Laplacian matrix of a graph on vertices in $T$. To simplify our exposition, we let $\mathbf{SC}(G,T)$ denote both the Laplacian and its corresponding graph. An important property of Schur complement which we exploit in this work is to view the Schur complement as a collection of random walks. This particular feature will be discussed in more detail in Section 3. The key role of Schur complements in our algorithms stems from the fact that they can be viewed as vertex sparsifiers that preserve pairwise effective resistances (see e.g., [24]).

Fact 2.2** (Vertex Resistance Sparsifier).**

For any graph $G=(V,E)$, any subset of vertices $T$, and any pair of vertices $u,v\in T$,

[TABLE]

Spectral Approximation

Definition 2.3** (Spectral Sparsifier).**

Given a graph $G=(V,E)$ and $\epsilon\in(0,1)$, we say that a graph $H=(V,E^{\prime})$ is a $(1\pm\epsilon)$-spectral sparsifier of $G$ (abbr. $H\approx_{\epsilon}G$) if $E^{\prime}\subseteq E$, and for all $\mathbf{x}\in\mathbb{R}^{n}$

[TABLE]

In the dynamic setting, Abraham et al. [2] recently showed that $(1\pm\epsilon)$-spectral sparsifiers of a dynamic graph $G$ can be maintained efficiently. This algorithm will be invoked in several places throughout this paper.

Lemma 2.4** ([2], Theorem 4.1).**

*Given a graph $G$ with polynomially bounded edge weights, with high probability, we can dynamically maintain a $(1\pm\epsilon)$-spectral sparsifier of size $\tilde{O}(n\epsilon^{-2})$ of $G$ in $O(\log^{9}n\epsilon^{-2})$ expected amortized time per edge insertion or deletion. The running time guarantees hold against an oblivious adversary. *

The above result is useful because matrix approximations also preserve approximations of their quadratic forms. As a consequence of this fact, we get the following lemma.

Lemma 2.5**.**

If $H$ is a $(1\pm\epsilon)$-spectral sparsifier of $G$, then for any pair of vertices $u$ and $v$

[TABLE]

2.1 Projection matrix and its properties

We next define a matrix that naturally appears when performing Gaussian elimination on the non-terminal vertices. Concretely, given a graph $G=(V,E)$ and terminals $T\subseteq V$, the matrix-projection of the non-terminals $F=V\setminus T$ onto $T$ is given by

[TABLE]

We now review some useful properties about the matrix projection $\mathbf{P}({T})$. Consider the Laplacian system $\mathbf{L}\mathbf{x}=\mathbf{b}$, where $\mathbf{L}$ is partitioned into block-matrices as in Equation (1). This in turn partitions the solution vector into non-terminals and terminals, i.e., $\mathbf{x}=\begin{bmatrix}\mathbf{x}_{F}\;\mathbf{x}_{T}\end{bmatrix}^{\top}$.

Lemma 2.6**.**

Let $\mathbf{x}_{T}$ be a solution vector such that $\mathbf{SC}(G,T)\mathbf{x}_{T}=\mathbf{P}({T})\mathbf{b}$. Then there exists an extension $\mathbf{x}$ of $\mathbf{x}_{T}$ such that $\mathbf{L}\mathbf{x}=\mathbf{b}$.

The following lemma draws a connection between the projection matrix and certain random walk probabilities which will allow us to take a combinatorial view on several problems we deal with.

Lemma 2.7**.**

Consider a graph $G=(V,E)$. For any subset of vertices $T\subseteq V$, a vertex $v\in T$, and a vertex $u\in F=V\setminus T$, let $\mathbb{P}_{u}\left[t_{v}<t_{T\setminus v}\right]$ be the probability that the random walk that starts at $u$ hits $v$ before hitting any other vertex from $T\setminus v$. Then we have that

[TABLE]

In fact, $\{\mathbf{P}({T})\boldsymbol{1}_{u}\}_{v\in T}$ is a probability distribution for any fixed vertex $v\in F$.

Given a demand vector $\mathbf{b}\in\mathbb{R}^{n}$, we say that $\mathbf{P}({T})\cdot\mathbf{b}$ is the projection of $\mathbf{b}$ onto $T$. In general, the projection of $\mathbf{b}$ is shorter than the original vector $\mathbf{b}$. However, for the sake of exposition, often we consider $\mathbf{P}({T})\cdot\mathbf{b}$ to be an $n$-dimensional vector by assuming that all coordinates in $F$ are [math].

Lemma 2.8**.**

Consider a graph $G=(V,E)$. Let $T\subseteq V$ be a subset of vertices, and let $u\in F=V\setminus T$. Consider the demand vector $\boldsymbol{1}_{u}-\mathbf{P}({T})\boldsymbol{1}_{u}$ that requests to send one unit of flow from $u$ to $T$ according to the probability distribution $\{\mathbf{P}({T})\boldsymbol{1}_{u}\}_{v\in T}$. Then the minimum energy needed to route this demand is given by

[TABLE]

3 Overview

The core building block of our algorithms is a fast routine that generates and maintains an approximate Schur complement onto a set of terminals $T$ of our choice under insertion and deletions of edges as well as terminal additions, with all of these operations being supported in sub-linear time. One of the key ideas is to view to the Schur complement as a sum of random walks, and then observe that sampling exactly one walk per edge in the original graph already yields the desired object. Concretely, we build upon ideas introduced in sparsifying random walk polynomials [14], and Schur complements [40, 19] to show that it suffices to keep a union of these walks. The following result is implicit in these works, and we will review its proof in the full version of this paper.

Theorem 3.1**.**

Let $G=(V,E,w)$ be an undirected, weighted multi-graph with a subset of vertices $T$. Furthermore, let $\epsilon\in(0,1)$, and let $\rho$ be some parameter related to the concentration of sampling given by

[TABLE]

Let $H$ be an initially empty graph, and for every edge $e=(u,v)$ of repeat $\rho$ times the following procedure:

Simulate a random walk starting from $u$ until it first hits $T$ at vertex $t_{1}$, 2. 2.

Simulate a random walk starting from $v$ until it first hits $T$ at vertex $t_{2}$, 3. 3.

Combine these two walks (including $e$) to get a walk $u=(t_{1}=u_{0},\ldots,u_{\ell}=t_{2})$, where $\ell$ is the length of the combined walk. 4. 4.

Add the edge $(t_{1},t_{2})$ to $H$ with weight

[TABLE]

The resulting graph $H$ satisfies $\mathbf{L}_{H}\approx_{\epsilon}\mathbf{SC}(G,T)$ with high probability.

The output approximate Schur complement of $H$ onto $T$ has up to $\tilde{O}(m\epsilon^{-2})$ edges, and thus is very dense to be leveraged as a sparsifier for our applications. Fortunately, there already exist efficient dynamic spectral sparsifiers, and we can always afford to keep a sparsifier $\tilde{H}$ of $H$ whose size is only $\tilde{O}(|T|\epsilon^{-2})$.

The performance of our data structure depends on how fast we can generate the random walks used to create $H$. Note that even on the length $n$ path with terminals $T$ concentrated on one end, the lengths of these walks may be as long as $\Omega(n^{2})$. To overcome this we shorten the walks by augmenting $T$ with roughly $O(\beta m)$ random vertices from a carefully chosen distribution. This random augmentation of $T$ ensures that any vertex $v$ in $G$ is roughly $O(\beta^{-1})$ apart from a vertex in $T$, and then our problem reduces to understanding the rate at which a random walk spreads among distinct edges. Concretely, our goal is to efficiently generate the first $k$ distinct edges visited by a walk in $G$. We distinguish the following cases.

For unweighted graphs, we utilize a result by Barnes and Feige [6] that shows that with high probability a walk reaches $k$ distinct edges in about $k^{2}$ steps. 2. 2.

For weighted graphs, we employ an event driven simulation of walks. Specifically, by computing the exit probability on the current set of edges visited so far, we sample the first $k$ new edges reached by the walk in $\mathrm{poly}(k)$ time. Then, because we know the order that each edges is first reached, the first among them that belongs to $T$ gives the intersection of the walk with $T$.

Following Point (1), our dynamic Schur complement data-structure $H$ with respect to a randomly augmented $T$ is initialized by generating for each edge $e\in E$, $\rho$ random walk of length roughly $\beta^{-2}$. This operation costs roughly $O(m\beta^{-2})$. We then make the observation that the ability to add terminals into $T$ means we only need to consider insertions/deletions between vertices in $T$. Specifically, for each affected edge we append its endpoints to $T$. A further advantage of this approach is that additions to $T$ only shorten random walks in $H$, and the cost of shortening or truncating these random walks in $H$ can be charged to the cost of constructing them during the initialization. Thus, it follows that we can support terminal additions, and thus insert or delete edges in $O(1)$ amortized time. Maintaining a sparsifier $\tilde{H}$ of $H$ introduces only polylogarithmic overheads, so this step does not affect much our running times. We next discuss the applicability of this result.

The data-structure we presented readily gives a sub-linear dynamic Laplacian solver for the case where $\mathbf{b}$ has small support, namely fewer than $\beta m$ vertices of $\mathbf{b}$ are non-zero. This can be accomplished by simply appending the entries of $\mathbf{b}$ (more precisely, their corresponding vertices) to the Schur complement $H$, and solving the system on $H$ upon receipt of an index query. The resulting solution vector can then be lifted back to the original Laplacian using Lemma 2.6. However, note that our data-structure can only support up to $O(|T|)=O(\beta m)$ operations if we want to keep the the size of $H$ small. Thus, to limit the growth in $|T|$ we periodically rebuild the entire data structure (i.e., we resample the set of new terminals completely) after $\beta m$ operations, which in turn gives an amortized update time of $O(m\beta^{-2}/(\beta m))=O(\beta^{-3})$. Combining this with the bound of $O(\beta m)$ on the query time we obtain the following trade-off

[TABLE]

which is minimized when $\beta=m^{-1/4}$, thus giving an update and query time of $O(m^{3/4})$.

So it remains to address the case where $\mathbf{b}$ has a large number of non-zero entries. We overcome this difficulty by projecting this vector onto the current set of terminals $T$ using the matrix $\mathbf{P}({T})$ and analyzing the error incurred by this projection. Our main observation is that the standard notion of error in Laplacian solvers, namely the $\mathbf{L}$-norm, corresponds to energies of electrical flows. This allows us to incur error in some of the $\mathbf{b}(u)$ values and then bound the energy of fixing them. To find such flows, we once again consider our problem from a random walk perspective, namely we view the projection of $\mathbf{b}$ onto $T$ being equivalent to moving $\mathbf{b}$ around via random walks (Lemma 2.7). As such walks are short on unweighted graphs, we can relate their energies to the length of the walks times $\mathbf{b}(u)^{2}$ (Lemma 2.8).

One final obstacle is that if we move some vertex $u$ from outside of $T$ into $T$, the walks affected may be from multiple $\mathbf{b}(u)$s. To address this, we bound the ‘load’ of a vertex, defined as the number of walks that go through it, by the total length of the walks. The latter follows from the uniform distribution of random walks being stationary. Thus, as long as we picked $T$ so that all the entries in $V\setminus T$ have small magnitudes, each move of some $u$ into $T$ incurs some small error. Bounding the accumulation of such errors, and rebuilding appropriately gives the overall dynamic solver result.

One application of the dynamic Laplacian solver is that we can maintain the energy of electrical flow for routing $\mathbf{b}$. This can also be viewed as an extension of our dynamic effective resistances data-structure, which can only maintain the energy of electrical flows for $\mathbf{b}$ with two non-zeros. Some further extensions in this direction that we believe would be useful are providing implicit access to the dual electrical flows, as well as finding the $k$ largest entries either in the flow edges or the solution vector $\mathbf{x}$. However, such extensions will likely require a better understanding of the graph sparsifier component [2], which is treated as a black box in this paper.

For dynamically maintaining effective resistance in unweighted graphs, we essentially follow the same approach as with the dynamic solver for small support demand vectors, and thus obtain a running time of $O(m^{3/4})$ on both update and query time. For weighted graphs, we employ the weighted dynamic Schur complement algorithm (following Point(2)) which gives slightly weaker guarantees, namely a bound of $\tilde{O}(m^{5/6})$ on the update and query time. Interestingly, this weighted version has another immediate advantage; by running the data-structure on the output of a dynamic spectral sparsifier (Lemma 2.4), we obtained a bound of $\tilde{O}(n^{5/6})$ per operation, which is truly sub-linear irrespective of graph density.

4 Dynamic Schur Complement

In this section we show how to dynamically maintain approximate Schur complements. We first restrict our attention to unweighted graphs (i.e., prove Lemma 1.1), and then discuss how this result extends to the weighted case. We also present one of the applications of our data-structure, namely dynamic maintenance of effective resistance on unweighted (Theorem 1.3).

4.1 Dynamic Schur Complement on Unweighted Graphs

In this section we design a data-structure for maintaining approximte Schur complements under the assumption that the dynamic graph remains unweighted throughout the updates. Specifically, we have the following lemma.

Lemma 4.1** (Restatement of Lemma 1.1).**

Given an error threshold $\epsilon>0$, an unweighted undirected multi-graph $G=(V,E)$ with $n$ vertices, $m$ edges, a subset of terminal vertices $T^{\prime}$ and a parameter $\beta\in(0,1)$ such that $|T^{\prime}|=O(\beta m)$, there is a data-structure DynamicSC$(G,T^{\prime},\beta)$ for maintaining a graph $\tilde{H}$ with $\mathbf{L}_{\tilde{H}}\approx_{\epsilon}\mathbf{SC}(G,T)$ for some $T$ with $T^{\prime}\subseteq T$, $|T|=O(\beta m)$, while supporting $O(\beta m)$ operations in the following running times:

• •

Initialize$(G,T^{\prime},\beta)$*: Initialize the data-structure, in $\tilde{O}(m\beta^{-2}\epsilon^{-4})$ expected amortized time. *

• •

Insert$(u,v)$: Insert the edge $(u,v)$ to $G$ in $\tilde{O}(1)$ amort. time.

• •

Delete$(u,v)$: Delete the existing edge $(u,v)$ from $G$ in $\tilde{O}(1)$ amortized time.

• •

AddTerminal$(u)$*: Add $u$ to $T^{\prime}$ in $\tilde{O}(1)$ amortized time. *

Our guarantees hold against an oblivious adversary.

To prove the lemma above, we first review the interpretation of Schur Complements using random walks, and then discuss how to generate and maintain these walks under edge updates and addition of terminal vertices.

Given a graph $G=(V,E)$ and a subset of terminals $T$ recall that $\mathbf{SC}(G,T)$ was defined using an algebraic expression that involved the Laplcian of $G$. However, since it is still unclear how to exploit this expression in the dynamic setting we instead take a different, more ‘combinatorial’, view on $\mathbf{SC}(G,T)$. Concretely, we will interpret $\mathbf{SC}(G,T)$ as a collection of random walks, each starting at an edge of $G$ and terminating in $T$, as described in Theorem 3.1.

Let $H$ be the output graph from the construction in Theorem 3.1. Recall that $H$ is an approximate Schur Complement onto $T$ that has up to $\rho m=\tilde{O}(m\epsilon^{-2})$ edges (that is, $\rho$ for each edge in $G$, where $\rho=O(\log n\epsilon^{-2})$ is the sampling parameter). As we will next show, $H$ does not change too much (in amortized sense) upon inserting or deleting an edge in $G$. We will be able to maintain $H$ such that the distribution of $H$ is the same as $H(G)$ of the current graph $G$. Therefore, we can maintain these changes using a dynamic spectral sparsifier $\tilde{H}$ of $H$ (Lemma 2.4), and whenever a query comes, we answer it on $\tilde{H}$ in $\tilde{O}(\left|T\right|\epsilon^{-2})=\tilde{O}(\beta m\epsilon^{-2})$ time.

While it is widely known how to generate random walks efficiently, we note that the length of the walks generated in Theorem 3.1 could be prohibitively large if $T$ is being picked arbitrarily. To see this, recall our example where we considered a path of length $n$ with terminals $T$ being places in one end. The length of such random walks may be as long as $\Omega(n^{2})$. To shorten these random walks, we augment $T$ with a random subset of vertices. Coming back to the path example, $\beta n$ uniformly random vertices will be roughly $\beta^{-1}$ apart, and random walks will reach one of these $\beta n$ vertices in about $\beta^{-2}$ steps. Because $G$ could be a multi-graph, and we want to support queries involving any vertex, we pick $T$ as the end points of a uniform subset of edges. A case that illustrates the necessity of this choice is a path except one edge has $n$ parallel edges. In this case it takes $\Theta(n)$ steps in expectation for a random walk to move away from the end points of that edge. This choice of $T$ completes the definition of our data structure, which we summarize in Algorithm 1.

The performance of our data structures hinge upon the properties of the generated random walks. We start by formalizing such a structure involving the set of the augmented terminals described above while parameterizing it with a more general probability $\beta$ for including the endpoints of the edges.

Definition 4.2** ($\beta$-shorted walks).**

Let $G$ be an weighted, undirected multi-graph and $\beta\in(0,1)$ a parameter. A collection of $\beta$-shorted walks $W$ on $G$ is a set of random walks created as follows:

Choose a subset of terminal vertices $T$, obtained by including the endpoints of each edge independently with probability at least $\beta$. 2. 2.

For each edge $e\in E$, generate $\rho=O(\log n\epsilon^{-2})$ walks from its endpoints either until $\Omega(\beta^{-1}\log{n})$ different edges have been hit, they reach $T$, or they visited each edge that is in the same connected component as $e$.

As we will shortly seee, the main property of the collection $W$ is that its random walks are short. Moreover, we will also prove that all walks in $W$ will reach $T$ with high probability. These guarantees are summarized in the following theorem.

Theorem 4.3**.**

Let $G=(V,E)$ be any undirected multi-graph, and $\beta\in(0,1)$ a parameter. Any set of $\beta$-shorted walks $W$, as described in Definition 4.2, satisfies the following:

• •

With high probability, any random walk in $W$ starting in a connected component containing a vertex from $T$ terminates at a vertex in $T$.

Note that Theorem 4.3 is conditioned upon the connected component having a vertex in $T$: this is necessary because walks stay inside a connected component. However, this does not affect our queries: our data-structure has an operation for making any vertex $u$ a terminal, which we call during each query to ensure both $s$ and $t$ are terminal vertices. Such an operation interacts well with Theorem 4.3 because it can only increase the probability of an edge’s endpoints being chosen.

Proving the theorem requires to determine the rate at which a random walk visits at least $\beta^{-1}\log n$ edges. It turns out that a random walk of length $\tilde{O}(\beta^{-2})$ is highly likely to achieve this. For formally showing this, we need the following result by Barnes and Feige [6].

Theorem 4.4** ([6], Theorem 1.2).**

There is an absolute constant $c_{BF}$ such that for any undirected unweighted connected multi-graph $G$ with $n$ vertices and $m$ edges, any vertex $u$ and any value ${\hat{{m}}}\leq m$, the expected time for a random walk starting from $u$ to visit at least ${\hat{{m}}}$ distinct edges is at most $c_{BF}{\hat{{m}}}^{2}$.

The above theorem can be amplified into a with high probability bound by repeating the walk $O(\log{n})$ times.

Corollary 4.5**.**

In any undirected unweighted connected multi-graph $G$ with $m$ edges, for any starting vertex $u$, any length $\ell\leq m^{2}$, and a parameter $\delta\geq 1$, a walk of length $c_{BF}\cdot\delta\cdot\ell\log n$ from $u$ visits at least $\ell^{1/2}$ distinct edges with probability at least $1-n^{-\delta}$.

Proof.

We can view each such walk as a concatenation of $\delta\log n$ sub-walks, each of length $2\cdot c_{BF}\cdot\ell$.

We call a sub-walk good if the number of distinct edges that it visits is at least $\ell^{1/2}$. Applying Markov’s inequality to the result of Theorem 4.4, a walk takes more than $2\cdot c_{BF}\cdot\ell$ steps to visit $\ell^{1/2}$ distinct edges with probability at most $1/2$.

This means that each subwalk fails to be good with probability at most $1/2$. Thus, the probability that all subwalks fail to be good is at most $2^{-\delta\log n}=n^{-\delta}$. The result then follows from an union bound over all starting vertices $u\in V$. ∎

We now have all the tools to prove Theorem 4.3.

Proof of Theorem 4.3.

For any walk $w$, we define $V(w)$ (respectively, $E(w)$) to be the set of distinct vertices (respectively, edges) that a walk $w$ visits. Consider a random walk $w$ that starts at $u$ of length

[TABLE]

where $\delta$ is a constant related to the success probability.

If the connected component containing the walk has fewer than

[TABLE]

edges, then Corollary 4.5 gives that we have covered this entire component with high probability, and the guarantee follows from the assumption that this component contains a vertex of $T$.

Otherwise, we will show that $w$ reached enough edges for one of their endpoints to be picked to be picked into $T$ with high probability. The key observation is that because $w$ is generated independently from $T$, we can bound the probability of this walk not hitting $T$ by first generating $w$, and then $T$. Specifically, for any size threshold $z$, we have

[TABLE]

By Corollary 4.5 and the choice of $\ell$, if we set

[TABLE]

then the first term in Equation (2) is bounded by $n^{-\delta}$. For bounding the second term, we can now focus on a particular walk $\widehat{w}$ that visits at least $\delta\cdot\beta^{-1}\cdot\log{n}$ distinct edges, i.e.,

[TABLE]

Recall that we independently added the end points of each of these edges into $T$ with probability $\beta$. If any of them is selected, we have a vertex that is both in $V(\widehat{w})$ and $T$. Thus, the probability that $T$ contains no vertices from $V(\widehat{w})$ is at most

[TABLE]

which completes the proof. ∎

Corollary 4.5 together with Theorem 4.3 yield the following lemma.

Lemma 4.6**.**

Algorithm 1 runs in $\tilde{O}(m\beta^{-2}\epsilon^{-2})$ time and outputs a graph $H$ with $\mathbf{L}_{H}\approx_{\epsilon}\mathbf{SC}(G,T)$, with high probability.

Proof.

By Corollary 4.5, the length of each walk generated in Algorithm 1 is bounded by $O(\beta^{-2}\log^{3}n)$. In addition, note that each step in a random walk can be simulated in $O(1)$ time. This is due to the fact that we can sample an integer in $[0,n-1]$ by drawing $x\in[0,1]$ uniformly and taking $\lfloor xn\rfloor$. Combining these with the fact that the algorithm generates $\rho m=\tilde{O}(m\epsilon^{-2})$ walks, it follows that the running time of the algorithm is dominated by $\tilde{O}(m\beta^{-2}\epsilon^{-2})$.

Note that the collection of generated walks form the set $W$ of $\beta$-shorted walks. By Theorem 4.3, with high probability, each of the walks that starts at a component containing a vertex in $T$ hits $T$. Conditioning on the latter, Theorem 3.1 gives that with high probability, $\mathbf{L}_{H}\approx_{\epsilon}\mathbf{SC}(G,T)$. ∎

Handling edge updates and terminal additions.

We start by observing that there is always a one-to-one correspondence between the collection of $\beta$-shorted walks $W$ and our approximate Schur complement $H$. Accordingly, our primary concern will be supporting the Insert, Delete, and AddTerminal operations in the collection $W$. However, as $W$ undergoes changes, we need to efficiently update the sparsifier $H$. To handle these updates, we would ideally have efficient access to walks in $W$ that are affected by the corresponding updates.

To achieve this, we index into walks that utilize a vertex or an edge, and thus set up a reverse data structure pointing from vertices and edges to the walks that contain them. The following lemma says that we can modify this representation with minimal cost.

Lemma 4.7**.**

For the collection of $\beta$-shorted walks $W$, let $W_{v}$ and $W_{e}$ be the specific walks of $W$ that contain vertex $v$ and edge $e$, respectively. We can maintain a data structure for $W$ such that for any vertex $v$ or edge $e$ it reports

• •

*all walks in $W_{v}$ or $W_{e}$ in $O(|W_{v}|)$ or $O(|W_{e}|)$ time, respectively, *

with an additional $O(\log{n})$ overhead for any changes made to $W$.

Proof.

For every vertex (respectively, edge), we can maintain a balanced binary search tree consisting of all the walks that use it in time proportional to the number of vertices (respectively, edges) in the walks. Supporting rank and select operations on such trees then gives the claimed bound. ∎

As a result, any update made to the collection of walks can be updated in the approximate Schur complement $H$ generated from these walks in $O(\log n)$ time. We now have all the necessary ingredients to prove Lemma 1.1.

Proof of Lemma 1.1.

We give a two-level data-structure for dynamically maintaining Schur complements on unweighted graphs, which keeps at any time a terminal set $T$ of size $\Theta(m\beta)$. This entails maintaining

an approximate Schur complement $H$ of $G$ with respect to $T$ (Theorem 3.1), 2. 2.

a dynamic spectral sparsifier $\tilde{H}$ of $H$ (Lemma 2.4).

We implement the procedure Initialize by running Algorithm 1, which produces a graph $H$ and then computing a spectral sparsifier $\tilde{H}$ of $H$ using Lemma 2.4. Note that by construction of our data-structure, every update in $H$ will be handled by the black-box dynamic sparsifier $\tilde{H}$.

As we will shortly see, operations Insert and Delete will be reduced to adding terminals to the set $T$. Thus, the bulk of our effort is devoted to implementing the procedure AddTerminal. Let $u$ be a non-terminal vertex that we want to append to $T$. We set $T\leftarrow T\cup\{u\}$, and then shorten all the walks at the first location they meet $u$. This shortening of walks induces in turn edge insertions and deletions to $H$, which are then processed by $\tilde{H}$. The pseudocode for this operation is summarized in Algorithm 2. To quickly locate the first appearances of $u$ in the random walks from $W$, we make use of the data-structure from Lemma 4.7. Let us first describe the construction of such data-structure during the preprocessing phase. Let $W_{u}$ be the balanced binary search tree consisting of all the walks that use the vertex $u$ in $W$. Fix $w\in W_{u}$. For any $t\geq 0$, if $w$ visits $u$ after $t$ steps, we check whether $W_{u}$ contains $w$ or not. If the latter holds, we know that $u$ has appeared before in $w$ and we do not need to add $w$ to $W_{u}$. Otherwise, we add $w$ to $W_{u}$ as this is the first time the walk $w$ visits $u$. After locating the first appearances of $u$, we cut the walks in these locations, delete the corresponding affected walks (together with their weight from $H$), and insert the new shorter walks to $H$. Note that we can simply use arrays to represent each random walk in $W$.

For implementing operations Insert and Delete we proceed as follows. Specifically, upon insertion or deletion of an edge $e=(u,v)$ in $G$, we append both $u$ and $v$ to the terminal set $T$. Now, all the walks that pass through $u$ or $v$ in $W$ must be shorten at the first location they meet $u$ or $v$. For inserting an edge $(u,v)$ with weight $\boldsymbol{\mathit{w}}_{uv}$ in $G$ (in fact, Lemma 1.1 is restricted to $\boldsymbol{\mathit{w}}_{u,v}=1$), we simply add $\rho$ trivial random walks (i.e., the edge $(u,v)$) of weight $\frac{\boldsymbol{\mathit{w}}_{uv}}{\rho}$ to $H$ (which sum up to the edge $(u,v)$ itself). For deleting the edge $(u,v)$ with weight $\boldsymbol{\mathit{w}}_{uv}$ from $G$, simply delete these $\rho$ random walks between $u$ and $v$ in $H$ (which exist since we guaranteed that $u$ and $v$ are added as terminals to $H$).

We next analyze the performance of our data-structure. Let us start with the pre-processing time. First, by Lemma 4.6 we get that the cost for constructing $H$ on a graph with $m$ edges is bounded by $\tilde{O}(m\beta^{-2}\epsilon^{-2})$. Next, since $H$ has $\tilde{O}(m\epsilon^{-2})$ edges, constructing $\tilde{H}$ takes $\tilde{O}(m\epsilon^{-4})$ time. Thus, the amortized time of Initialize operation is bounded by $\tilde{O}(m\beta^{-2}\epsilon^{-4})$.

Following the above discussion, for analyzing the cost of the update operations it suffices to bound the time for adding a vertex to $T$, which in turn (asymptotically) bounds the update time for edge insertions and deletions. The main observation we make is that adding a vertex to $T$ only shortens the existing walks, and Lemma 4.7 allows us to find such walks in time proportional to the amount of edges deleted from the walk. Since the walk needed to be generated in the Initialize operation, the deletion of these edges takes equivalent time to generating them. Moreover, we note that (1) handling the updates in $\tilde{H}$ induced by $H$ introduces additional $O(\mathrm{poly}(\log n)\epsilon^{-2})$ overheads, and (2) adding or deleting $\rho$ edges until the next rebuild costs $\tilde{O}(\beta m\epsilon^{-2})$, since we process only up to $\beta m$ operations. These together imply that the amortized cost for adding a terminal can be charged against the pre-processing time, which is bounded by $\tilde{O}(m\beta^{-2}\epsilon^{-4})$, up to poly-logarithmic factors. Thus, it follows that the operations AddTerminal, Insert and Delete can be implemented in $\tilde{O}(1)$ amortized update time. ∎

4.2 Dynamic All-Pairs Effective Resistance on Unweighted Graphs

In this section we present one of the applications of our dynamic Schur complement data structure for unweighted graphs. Concretely, we design a dynamic algorithm that supports an intermixed sequence of edge insertions, deletion and pair-wise resistance queries, and returns a $(1\pm\epsilon)$-approximation to all the resistance queries.

We start by reviewing two natural attempts for solving this problem.

• •

First, since spectral sparsifiers preserve effective resistances (Lemma 2.5), we could dynamically maintain a spectral sparsifier (Lemma 2.4), and then compute the $(s,t)$ effective resistance on this sparsifier. This leads to a data structure with $\mathrm{poly}(\log n,\epsilon^{-1})$ update time and $\tilde{O}(n\epsilon^{-2})$ query time.

• •

Second, by the preservation of effective resistances under Schur complements (Fact 2.2), we could also utilize Schur complements to obtain a faster query time among a set of $\beta m$ terminals, $T$, for some reduction factor $\beta\in(0,1)$, at the expense of a slower update time. Specifically, after each edge update, we recompute an approximate Schur complement of the sparsifier onto $T$ in time $\tilde{O}(m\epsilon^{-2})$ [18], after which each query takes $\tilde{O}(\beta m\epsilon^{-2})$ time.

The first approach obtains sublinear update time, while the second approach gives sublinear query time. Our algorithm stems from combining these two methods, with the key additional observation being that adding more vertices to $T$ makes the Schur complement algorithm more local. Specifically, using Lemma 1.1 leads to a data-structure for dynamically maintaining all-pairs effective resistances.

Proof of Theorem 1.3.

Let $\mathcal{D}(\tilde{H})$ denote the data structure that maintains a dynamic (sparse) Schur complement $\tilde{H}$ of $G$ (Lemma 1.1). Since $\mathcal{D}(\tilde{H})$ supports only up to $\beta m$ operations, we rebuild $\mathcal{D}(\tilde{H})$ on the current graph $G$ after such many operations. Note that the operations Insert and Delete on $G$ are simply passed to $\mathcal{D}(\tilde{H})$. For processing the query operation $\textsc{EffectiveResistance}(s,t)$, we declare $s$ and $t$ terminals (using the operation AddTerminal of $\mathcal{D}(\tilde{H})$), which ensures that they are both now contained in $\tilde{H}$. Finally, we compute the (approximate) effective resistance between $s$ and $t$ in $\tilde{H}$ using Lemma 2.1.

We now analyze the performance of our data-structure. Recall that the insertion or deletion of an edge in $G$ can be supported in $\tilde{O}(1)$ expected amortized time by $\mathcal{D}(\tilde{H})$. Since our data-structure is rebuilt every $\beta m$ operations, and rebuilding $\mathcal{D}(\tilde{H})$ can be implemented in $\tilde{O}(m\beta^{-2}\epsilon^{-4})$, it follows that the amortized cost per edge insertion or deletion is

[TABLE]

The cost of any $(s,t)$ query is dominated by (1) the cost of declaring $s$ and $t$ terminals and (2) the cost of computing the $(s,t)$ effective resistance to $\epsilon$ accuracy on the graph $\tilde{H}$. Since (1) can be performed in $\tilde{O}(1)$ time, we only need to analyze (2). We do so by first giving a bound on the size of $T$. To this end, note that each of the $m$ edges in the current graph adds two vertices to $T$ with probability $\beta$ independently. By a Chernoff bound, the number of random augmentations added to $T$ is at most $2\beta m$ with high probability. In addition, since $\mathcal{D}(\tilde{H})$ is rebuilt every $\beta m$ operations, the size of $T$ never exceeds $4\beta m$ with high probability. The latter also bounds the size of $\tilde{H}$ by $\tilde{O}(\beta m\epsilon^{-2})$ and gives that the query cost is $\tilde{O}(\beta m\epsilon^{-4})$.

Combining the above bounds on the update and query time, we obtain the following trade-off

[TABLE]

which is minimized when $\beta=m^{-1/4}$, thus giving an expected amortized update and query time of

[TABLE]

4.3 Dynamic Schur Complement on Weighted Graphs

In this section we present an extension of Lemma 1.1 to weighted graphs while slightly increasing the running time guarantees. Concretely, we prove the following lemma.

Lemma 4.8**.**

Given an error threshold $\epsilon>0$, a weighted, undirected multi-graph $G=(V,E,\boldsymbol{\mathit{w}})$ with $n$ vertices, $m$ edges, a subset of terminal vertices $T^{\prime}$ and a parameter $\beta\in(0,1)$ such that $|T^{\prime}|=O(\beta m)$, there is a data-structure WeightedDynamicSC$(G,T^{\prime},\beta)$ for maintaining a graph $\tilde{H}$ with $\tilde{H}\approx_{\epsilon}\mathbf{SC}(G,T)$ for some $T$ with $T^{\prime}\subseteq T$, $|T|=O(\beta m)$, while supporting $O(\beta m)$ operations in the following running times:

• •

Initialize$(G,T^{\prime},\beta)$: Initialize the data-structure in $\tilde{O}(m\beta^{-4}\epsilon^{-4})$ expected amortized time.

• •

Insert$(u,v,w)$: Insert the edge $(u,v)$ with weight $w$ to $G$ in $\tilde{O}(1)$ amortized time.

• •

Delete$(u,v)$: Delete the existing edge $(uf,v)$ from $G$ in $\tilde{O}(1)$ amortized time.

• •

AddTerminal$(u)$: Add $u$ to $T^{\prime}$ in $\tilde{O}(1)$ amortized time.

While the extension of our data-structure to weighted graphs builds upon the ideas we used in the unweighted case, there are a few obstacles that force us to introduce new components in our algorithm in order to make such an extension feasible. To illustrate, consider the weighted graph shown in Figure 1 and recall that the running time our data-structure depends on the speed at which random walks visit distinct edges in a graph. Due to the structure of the edge weights, a random walk in this graph is expected to take $\Theta(n^{10})$ steps before hitting four different edges. This shows that the naive generation of random walks in weighted graphs may be computationally prohibitive for our purposes.

Fast Generation of Random Walks in Weighted Graphs.

To rectify the above issue, we make the important observation that it is not necessary to keep information for every single step of a random walk. Instead, it would suffice if we could efficiently determine the step at which the walk meets a new vertex along with the corresponding weight associated with the walk, which defines the edge weight that is added to the sparsifier. This high-level idea allows us to generate random walks much faster, and we next make this more precise.

Following the notation we used in the unweighted case, for an arbitrary vertex $v\in V$, a set of terminals $T\subseteq V$ and a parameter $\beta\in(0,1)$, a $\beta$-shorted walk with respect to $v$ and $T$ is a random walk that starts at a given vertex $v\in V$ and halts whenever $\Omega(\beta^{-1}\log n)$ different vertices have been hit, it reaches a vertex in $T$, or it has hit every edge in the connected component containing $v$. The main contribution of this section is summarized in the following lemma.

Lemma 4.9**.**

Let $G=(V,E,\boldsymbol{\mathit{w}})$ be an undirected, weighted graph with $\boldsymbol{\mathit{w}}_{e}=[1,n^{c}]$ for each $e\in E$, where $c$ is a positive constant. Let $T\subseteq V$ be a set of terminals and $v\in V$ be an arbitrary vertex. Then there is an algorithm that generates a $\beta$-shorted random walk with respect to $v$ and $T$ and approximates its corresponding weight up to a $\pm\epsilon$ relative error in $\tilde{O}(\beta^{-4}\epsilon^{-2})$ time.

We first give an intuition behind the algorithm in the above lemma and also describe how this algorithm interacts with other parts of our dynamic data-structure. Let $w=(w_{0},\ldots,w_{t})$ be a random walk that starts at an endpoint of an edge, and define

[TABLE]

to be its corresponding weight. Recall that before adding the walk $w$ to $H$, we must scale it proportionally to $1/s(w)$ (Theorem 3.1). Observe that throughout our dynamic algorithm, the only modification we might do to $w$ is to truncate it at the first location it meets a new vertex $u$ that is being declared a terminal. Moreover, after this modification, note that the old value of $s(w)$ is no longer valid and we need to extract $s(w)$ that corresponds to the new walk. To allow efficient access to such information, we can view the walk $w$ as being split into sub-walk segments by the first locations $w$ meets new vertices and store the weights of each such sub-walks. As we will next see, this bookkeeping alone allows us to proceed with the same algorithm as in the unweighted case.

We next give the three main components for implementing the algorithm stated in Lemma 4.9.

1. (A)

Sample the number of steps needed for a random walk $w$ to visit a new vertex. 2. (B)

Sample a new distinct vertex that $w$ hits, and its corresponding edge. 3. (C)

Sample the (approximate) weight of a random walk between two given vertices.

After describing each of them, we will see that their combination naturally leads to our desired result.

Let us first discuss A. For any $t\geq 0$, consider a $t$ step random walk $w$ and let $U=\{w_{0},\ldots,w_{t}\}$ be the set of distinct vertices that $w$ has visited up to step $t$. Define $u:=w_{t}\in U$ to be the current vertex of the walk $w$. Our goal to to efficiently sample the number of steps the walk $w$ needs to visit a vertex not in $U$. To this end, we start by introducing some useful notation. For any $i\geq 0$, let $\mathit{p}^{\mathrm{new}}(i)$ be the probability that $w$ meets a new vertex that is not in $U$ in $w_{t+1},\ldots,w_{t+i}$. For $v\in U$, let $\boldsymbol{\mathit{p}}_{i}(v)$ be the probability that $w^{(t+i)}=v$, conditioned on $w$ not having met any new vertex in $w_{t+1},\ldots,w_{t+i-1}$. Then it can be easily verified that both $\mathit{p}^{\mathrm{new}}(i)$ and $\boldsymbol{\mathit{p}}_{i}(v)$ are just linear combinations of $\mathit{p}^{\mathrm{new}}(i)$ and $\boldsymbol{\mathit{p}}_{i-1}(v)$

[TABLE]

Next, using the linearity of the recurrences in (5) and (4) we can find a matrix $\boldsymbol{\mathit{W}}$ of dimension $(k+1)\times(k+1)$, where $k=|U|$, satisfying the following equality

[TABLE]

The main advantage introducing such a matrix is that it allows us to efficiently compute $\mathit{p}^{\mathrm{new}}(i)$ and $\boldsymbol{\mathit{p}}_{i}$ using fast exponentiation via repeated squaring. Specifically, let $\boldsymbol{\mathit{p}}_{0}$ be a unit vector of dimension $k$, where for the current vertex $u$ of the walk $w$ we have that $\boldsymbol{\mathit{p}}_{0}(u)=1$, and [math] otherwise. Let $\hat{\boldsymbol{\mathit{p}}}_{0}=\left[\boldsymbol{\mathit{p}}_{0}\quad\mathit{p}^{\mathrm{new}}(0)\right]^{\top}$ be the extended $k+1$ dimension vector, where $\mathit{p}^{\mathrm{new}}(0)=0$. For any $i\geq 1$, repeatedly applying Equation (6) and letting $\hat{\boldsymbol{\mathit{p}}}_{i}:=\boldsymbol{\mathit{W}}^{i}\hat{\boldsymbol{\mathit{p}}}_{0}$ yields

[TABLE]

Using the above relation, we can use fast exponentiation via repeated squaring to compute $\mathit{p}^{\mathrm{new}}(i)$ in $O(k^{3}\log(i))$ time. This follows directly from the following well-known lemma, which we will exploit in a few other places throughout this work.

Lemma 4.10**.**

Let $A$ be a matrix of dimension $n\times n$, and $A^{i}$ denote the $i$-th power of $A$, for any $i\geq 1$. Then there is an algorithm that computes $A^{i}$ in $O(n^{3}\log(i))$ time.

We now have all the tools to describe the sampling procedure for computing the number of steps that the walk needs to visit a vertex that is distinct from the vertices in $U$. We accomplish this using a “binary search”-inspired subroutine, which works as follows. As an input, our algorithm is given a $(k+1)\times(k+1)$ matrix $\boldsymbol{\mathit{W}}$ (as defined in Equation (6)), the vector $\hat{\boldsymbol{\mathit{p}}}_{0}$, and an integer $M$, which is an upper-bound on the cover time of $G$. The algorithm also maintains variables $\ell,r,\ell p,rp$ with the following initialization $\ell\leftarrow 0$, $r\leftarrow M$, $\ell p\leftarrow 0$ and $rp\leftarrow 1$. As long as $(\ell\neq r)$, it defines the average $\eta=\lfloor(\ell+r)/2\rfloor$ and then proceeds to compute $\hat{\boldsymbol{\mathit{p}}}_{\eta}=W^{\eta}\hat{\boldsymbol{\mathit{p}}}_{0}$ using Lemma 4.10. Note that $\mathit{p}^{\mathrm{new}}(\eta)=\hat{\boldsymbol{\mathit{p}}}_{\eta}(k+1)$ by Equation (7). Finally, the algorithm uses $\mathit{p}^{\mathrm{new}}(\eta)$ to randomly decide whether $w$ meets a new vertex in the next $\eta$ steps or not. In other words, it updates the maintained variables using the rule below:

with probability $(\mathit{p}^{\mathrm{new}}(\eta)-\ell p)/(rp-\ell p)$, set $r\leftarrow\eta$, and $rp=\mathit{p}^{\mathrm{new}}(\eta)$, 2. 2.

otherwise, with probability $(rp-\mathit{p}^{\mathrm{new}}(\eta))/(rp-\ell p)$, set $\ell\leftarrow\eta+1$, $\ell p=\mathit{p}^{\mathrm{new}}(\eta)$.

If $(\ell=r)$, then the algorithm returns $\ell$. This procedure is summarized in Algorithm 3.

We next show the correctness of the above procedure. To do so, we first need the following notation. For a $t$-step random walk $w$ and a current vertex $u=w_{t}\in U$, let $X(u,U)$ be the smallest number of steps of steps needed for $w$ to visit a vertex not in $U$, i.e., $X(u,U)=\min\{i\mid i\geq 1,\leavevmode\nobreak\ w_{t+i}\not\in U\}$. Note that $X(u,U)$ is a random variable, and $X(u,U)\leq M$ by definition of $M$.

Lemma 4.11**.**

Let $w$ be a $t$-step random walk and $U$, $u=w_{t}\in U$, and $k=|U|$ be the number of distinct vertices $w$ has visited so far. For $\boldsymbol{\mathit{W}}$, $\hat{\boldsymbol{\mathit{p}}}_{0}$, and $M$ defined as above, BinarySearch$(\boldsymbol{\mathit{W}},\hat{\boldsymbol{\mathit{p}}}_{0},M)$ correctly samples $X(u,U)$, i.e., the number of steps $w$ needs to visit a vertex not in $U$, in $O(k^{3}\log^{2}M)$ time.

Proof.

By Equation (7) and Line 3 in Algorithm 3, note that $\mathit{p}^{\mathrm{new}}(\eta)$ is the probability that $w$ meets a new vertex in the fist $\eta$ steps. The correctness of BinarySearch can be proven using an inductive argument on the number of iterations of the while loop. Here, we just show the crucial parts for being able to apply such an argument. First, observe that right after Line 3 in the while loop, we have that

[TABLE]

The latter holds because $\ell p=\mathbb{P}_{X(u,U)}\left[X(u,U)\leq\ell\right]$ and $rp=\mathbb{P}_{X(u,U)}\left[X(u,U)\leq r\right]$, which in turn can be verified for each assignment of $\ell$ and $r$. Next, we prove that conditioning on $\ell\leq X(u,U)\leq r$ right after Line 3 in the while loop, Algorithm 3 samples $X(u,U)$ from the correct distribution. This is true when $(\ell=r)$, since the condition of the while loop is no longer satisfied and the algorithm returns $\ell$. If, however $(\ell\neq r)$, then we need to compute the following probabilities: (1) $\mathbb{P}_{X(u,U)}\left[X(u,U)\leq\eta\mid\ell\leq X(u,U)\leq r\right]$ and (2) $\mathbb{P}_{X(u,U)}\left[X(u,U)>\eta\mid\ell\leq X(u,U)\leq r\right]$. To determine (1), we get that

[TABLE]

The probability from case (2) can be shown similarly. Since Line 3 in Algorithm 3 updates the search boundaries $\ell$ and $r$ and their corresponding values $\ell p$ and $rp$ using probabilities (1) and (2), the correctness of the algorithm follows.

For the running time, observe that the number of iterations until the condition of the while loop is no longer satisfied is bounded by $O(\log M)$. Moreover, the running time of one iteration is dominated by the time needed to compute $\boldsymbol{\mathit{W}}^{\eta}\hat{\boldsymbol{\mathit{p}}}_{0}$. Since $\boldsymbol{\mathit{W}}$ is a $(k+1)\times(k+1)$ dimensional matrix and $\eta\leq M$, Lemma 4.10 implies that the matrix powering step can be computed in $O(k^{3}\log M)$. Thus, it follows that Algorithm 3 can be implemented in $O(k^{3}\log^{2}M)$ time. ∎

We next explain how to sample a new distinct vertex, and its corresponding edge of a $t$-step random walk $w$, i.e., we discuss component B. Let $X(u,U)$ be the index computed by BinarySearch routine. We first compute the probability dsitribution $\boldsymbol{\mathit{q}}$ over vertices in $U$ after performing the next $(X(u,U)-1)$ steps of the random walk $w$, conditioning on $w$ not leaving $U$. Next, we proceed to computing the probability distribution $\boldsymbol{\mathit{r}}$ over the edges in in $(U,V\setminus U)$ conditioning on $w_{0},\ldots,w_{t}$ and $w_{t+X(u,U)}$ being the first vertex not in $U$. Formally, for $v\in U$, $z\in V\setminus U$, we have

[TABLE]

Finally, we sample $(w_{t+X(u,U)-1},w_{t+X(u,U)})$ according to $\boldsymbol{\mathit{r}}$, where $w_{t+X(u,U)}$ is the first vertex not in $U$. The lemma below shows that we can efficiently sample from $\boldsymbol{\mathit{r}}$.

Lemma 4.12**.**

Let $w$ be a $t$-step random walk and let $U$ with $k=|U|$ be the set of distinct vertices $w$ has visited so far. Given the number of steps $X(u,U)$ needed for $w$ to visit a vertex not in $U$, there exists an algorithm that samples an edge leaving $U$, and the first vertex not in $U$ in $O(k^{3}\log M)$ time.

Proof.

We start by showing how to compute the distribution $\boldsymbol{\mathit{q}}$. To this end, recall that $\boldsymbol{\mathit{p}}_{i}(v)$ is the probability that $w_{t+i}=v$, conditioned on $w$ not having met any vertex different from $U$ in $w_{t+1},\ldots,w_{t+i-1}$. Thus, by Equation (6), we can use the fast exponentiation routine (Lemma 4.10) to compute the vector $\hat{\boldsymbol{\mathit{p}}}_{X(u,U)-1}=\boldsymbol{\mathit{W}}^{X(u,U)-1}\hat{\boldsymbol{\mathit{p}}}_{0}$. Since by Equation (7) we have that $\hat{\boldsymbol{\mathit{p}}}_{X(u,U)-1}(v)=\boldsymbol{\mathit{p}}_{X(u,U)-1}(v)$ for each $v\in U$, it follows that $\boldsymbol{\mathit{q}}(v)$ is exactly $\boldsymbol{\mathit{p}}_{X(u,U)-1}(v)$. Note that the running time for implementing this step is $O(k^{3}\log M)$ as $X(u,U)\leq M$.

We next describe how to efficiently sample from the distribution $\boldsymbol{\mathit{r}}$. First, it will be helpful to to sample a vertex $v\in U$ conditioning on the on $w_{t+X(u,U)}$ being the first vertex not in $U$. Specifically, we are interested in sampling a vertex $v\in U$ with probability

[TABLE]

For being able to efficiently sample from this distribution, we need to compute $\boldsymbol{\mathit{w}}(v,V\setminus U)$, which in turn may require examining up to $\Omega(n)$ edges incident to $v$. However, this is not sufficient for our purposes as our ultimate goal is to sample from $\boldsymbol{\mathit{r}}$ in time only proportional to $k$. To alleviate this, observe that $\boldsymbol{\mathit{w}}(v,V\setminus U)=(\boldsymbol{\mathit{d}}(v)-\sum_{z\in U}w_{v,z})$. Thus, maintaining weighted degree $\boldsymbol{\mathit{d}}(v)$ for each $v\in V$, allows us to compute $\boldsymbol{\mathit{w}}(v,V\setminus U)$ in $O(k)$ time. Similarly, rearranging the sums in the definition of $R$ we get

[TABLE]

which in turn implies that $R$ can be computed in $O(k^{2})$ time. The latter gives that the distribution defined in Euqation (9) can be computed in $O(k^{2})$ time. For sampling a vertex $v\in U$ from this distribution we simply generate a uniformly-random value $x\in[0,1]$, and then perform binary search on the prefix sum array of the probability distribution. Since computing the prefix sum array and performing binary search can be done in $O(k)$ and $O(\log n)$ time, respectively, we get that sampling $v\in U$ according to distribution defined in Equation (9) can be performed in $O(k^{2})$ time.

We next explain how to sample an edge $(v,z)$, where $z\in V\setminus U$ and $v\in U$ is the vertex we sampled from above. The probability distribution from which $(v,z)$ is sampled is as follows

[TABLE]

To see the idea behind this choice, note that Equation (10) combined with Equation (9) yields the distribution $\boldsymbol{\mathit{r}}$ as defined in Equation (8), which ensures that the edge is sampled correctly. However, one complication we face with is that $v$ may be incident to $\Omega(n)$ edges. Remember that for sampling an edge one needs access to the prefix sum array, which is expensive for our purposes. A natural attempt is to compute such an array during preprocessing. Nevertheless, this alone does not suffice as the set $U$ will change over the course of our algorithm. Instead, for every vertex $v\in V$, we maintain an augmented Balanced Binary Tree (BBT) on the edge weights incident to $v$. Augmented BBT is a data-structure that supports operations such as (1) computing prefix sums and (2) updating the edge weights incident to $v$, both in $O(\log n)$ time.

We employ the augmented BBT data-structure as follows. First, for each vertex $U$ and the sampled vertex $v\in U$, we update the weights of the edges from $v$ to $U$ to [math] in the augmented BBT of $v$. We then sample a uniformly-random value $x\in[0,W]$, and use the prefix sums computation in the tree to determine the range in which $x$ lies together with the corresponding edge $(w_{t+X(u,U)-1},w_{t+X(u,U)})$, where $w_{t+X(u,U)}$ is the first vertex not in $U$. After having sampled the edge, we undo all the changes we performed in the augmented BBT of $v$. It follows that sampling an edge according to Equation (10) can be implemented in $O(k\log n)$. Putting together the above running times, we conclude that sampling an edge leaving $U$ as well as the first vertex not in $U$ can be implemented in $O(k^{3}\log M)$ time. ∎

The last ingredient we need is an efficient way to sample the sum of weights in the random walk starting at $w_{t}$ and ending at $w_{t+X(u,U)}$, where $X(u,U)$ is the number of steps needed for the walk to leave the vertex set $U$ (Component C). In other words, we need to sample the following sum

[TABLE]

We accomplish this task by employing a doubling technique. To illustrate, for any pair of vertices $u,v\in V$ and $s(w)$ as defined in Equation (3), let

[TABLE]

be the probability mass function of $s(w)$ conditioning on (1) $w$ being a random walk that starts at $u$ and ends at $v$, i.e., $w\sim w_{u,v}$ and (2) length of the walk $\ell(w)$ is $\ell$ in $G$. Then it can be shown that

[TABLE]

where $*$ denotes the convolution between two probability mass functions. Equivalently, the convolution is the probability mass function of the sum of the two corresponding random variables. The above relation suggests that if (1) we have some approximate representation of the probability mass functions $f_{s(w),\ell/2}^{u,v}$ for all $u,v\in V$, and (2) we are able to compute the convolution of the two mass functions under such representation, we can produce approximations for $f_{s(w),\ell}^{u,v}$, where $u,v\in V$. This idea is formalized in the following lemma.

Lemma 4.13**.**

Let $G=(V,E,\boldsymbol{\mathit{w}})$ be a undirected, weighted graph with $\boldsymbol{\mathit{w}}_{e}=[1,n^{c}]$ for each $e\in E$, where $c$ is a positive constant. For any finite random walk $w$ of length $\ell$ with $\ell\leq n^{d}$, where $d$ is a positive constant, let $s(w)$ be the sum of the inverse of its edge weights, i.e.,

[TABLE]

Moreover, for any $u,v\in V$, let

[TABLE]

be the probability mass function of $s(w)$ conditioning on (1) $w$ being a random walk that starts at $u$ and ends at $v$, and (2) length of the walk $\ell(w)$ is $\ell$ in $G$. Then, for any pair $u,v\in V$, there exists an algorithm that that samples from $f_{s(w),\ell}^{u,v}$ and outputs a sampled $s(w)$ up to $\pm\epsilon$ relative error in $\tilde{O}(n^{3}\epsilon^{-2})$ time.

We finally describe a procedure that generates a $\beta$-shorted walk with respect to some vertex $v$ and set of terminals $T$. Concretely, the algorithm maintains (1) a set $U$, initialized to $\{v\}$, of the distinct vertices visited so far by a random walk $w$ starting at $v$, (2) the number of steps $t$ the walk $w$ has performed so far and (3) two lists $L_{w}$ and $L_{s}$, initially set to empty, containing the first occurrences of distinct vertices of $w$ and the weightes of the sub-walks induced by the distinct vertices, respectively. Next, as long as $w$ does not hit a vertex in $T$ or there are vertices in the component containing $v$ that are still not visited by $w$, for the next $\Theta(\beta^{-1}\log n)$ steps, the algorithm repeatedly generates a new vertex not in the current $U$ by using components A, B and C. In each iteration, the maintained quantities $U$, $t$, $L_{w}$ and $L_{s}$ are updated accordingly. Note that this procedure indeed outputs all necessary information we need from a $\beta$-shorted walk. A detailed implementation of the algorithm is summarized in Figure 4.

We now have all the necessary tools to prove Lemma 4.9.

Proof of Lemma 4.9.

We first show correctness. By Lemma 4.11 it follows that BinarySearch correctly samples the number of steps before a walk meets a new vertex. Next, Lemma 4.12 implies that the we can sample the new distinct vertex and its corresponding edge. Finally, by Lemma 4.13 we know that the weight of each sub-walk of a $\beta$-shorted walk is approximated within a $(1+\epsilon)$ relative error. Bringing these approximation together we get that the weight of the $\beta$-shorted walk itself is approximated within the same relative error.

We now analyse the running time of procedure GenerateSingleWalk. We start by bounding the cover time of $G$, which in turn bounds the number of steps for a random walk to meet a new vertex. To this end, note that it takes expected $O(m^{2})$ time to meet a vertex in the same component ([3]). Thus, if we perform a random walk of length $O(m^{3})$ we are guaranteed that it covers ever vertex in the component containing the starting vertex, with high probability.

Next, we analyze the running time for the steps executed within one iteration of the for loop. Observe that $k=|U|=O(\beta^{-1}\log n)$ at any point of time throughout our algorithm. The latter together with Lemma 4.11 give that it takes $O(k^{3}\log^{2}M)=\tilde{O}(\beta^{-3})$ time to sample the minimum number of steps for a random walk to visit a vertex not in $U$, where $M=O(m^{3})$ by the discussion above. Furthermore, by Lemma 4.12 we can sample the new vertex not in $U$, and its corresponding edge in $\tilde{O}(\beta^{-3})$ time. Finally, Lemma 4.13 implies that the weight $s(w)$ of the random sub-walk between the current vertex and the new generated vertex can be approximately sampled in $\tilde{O}(\beta^{-3}\epsilon^{-2})$ time. The latter holds because Lemma 4.13 is invoked on top of the graph $G[U]$ for which $|V(G[U])|=O(\beta^{-1}\log n)$. Combining the above running times, we get that one iteration can be implemented in $\tilde{O}(\beta^{-3}\epsilon^{-2})$ time. Since there are $O(\beta^{-1}\log n)$ iterations, we conclude that the overall running time of our procedure is $\tilde{O}(\beta^{-4}\epsilon^{-2})$.∎

We now present the procedure for generating a Schur complement on weighted graphs. The idea behind this algorithm is the same as in the unweighted setting, except that now we use GenerateSingleWalk to extract the information needed to simulate $\beta$-shorted walks. For the sake of completeness we summarize the details of this modified procedure in Algorithm 5.

Lemma 4.14**.**

Algorithm 5 runs in $\tilde{O}(m\beta^{-4}\epsilon^{-4})$ time and outputs a graph $H$ satisfying $\mathbf{L}_{H}\approx_{\epsilon}\mathbf{SC}(G,T)$, with high probability.

Proof.

We first bound the running time of Algorithm 5. By Lemma 4.9, the time needed to generate a $\beta$-shorted walk is $\tilde{O}(\beta^{-4}\epsilon^{-2})$. Combining this with the fact that the algorithm generates $\rho m=\tilde{O}(m\epsilon^{-2})$ walks, it follows that the running time of the algorithms is dominated by $\tilde{O}(m\beta^{-4}\epsilon^{-4})$.

We next show the correctness of our procedure. First, note that procedure GenerateSingleWalk generates a valid $\beta$-shorted walk with its weight being approximated up to a $(1+\epsilon)$ relative error (Lemma 4.9). Assume for now that there is an oracle that fixes this approximate weight of a walk to its original exact weight. Then the collection of generated walks from Algorithm 5 forms the set $W$ of $\beta$-shorted walks, and let $\hat{H}$ be the corresponding output graph. By Theorem 4.3, with high probability, each of the walks that starts at a component containing a vertex in $T$ hits $T$. Conditioning on the latter, Theorem 3.1 gives that with high probability, $\mathbf{L}_{\hat{H}}\approx_{\epsilon}\mathbf{SC}(G,T)$.

Finally, let $H$ be the graph where the edge weights are correct up to a $(\pm\epsilon)$ relative error. In other words, the weight of each edge $e$ in $H$ differs from the corresponding weight $\boldsymbol{\mathit{w}}_{\hat{H}}(e)$ in $\hat{H}$ by $\epsilon\boldsymbol{\mathit{w}}_{\hat{H}}(e)$. Summing over all the edges we get that $\mathbf{L}_{H}\approx_{\epsilon}\mathbf{L}_{\hat{H}}$. Since $\mathbf{L}_{\hat{H}}\approx_{\varepsilon}\mathbf{SC}(G,T)$ by the discussion above, we get that $\mathbf{L}_{H}\approx_{O(\epsilon)}\mathbf{SC}(G,T)$. Scaling $\epsilon$ appropriately completes the correctness. ∎

We now have all the necessary tools to present our dynamic algorithm for maintaining the collection of walks $W$ (equivalently, the approximate Schur complement $H$), on weighted graphs.

Proof of Lemma 4.8.

Similarly to the unweighted case, we give a two-level data-structure for dynamically maintaining Schur complements on weighted graphs. Specifically, we keep the terminal set $T$ of size $\Theta(m\beta)$. This entails maintaining

an approximate Schur complement $H$ of $G$ with respect to $T$ (Theorem 3.1), 2. 2.

a dynamic spectral sparsifier $\tilde{H}$ of $H$ (Lemma 2.4).

We implement the procedure Initialize by running Algorithm 5, which produces a graph $H$ and then computing a spectral sparsifier $\tilde{H}$ of $H$ using Lemma 2.4. Note that by construction of our data-structure, every update in $H$ will be handled by the black-box dynamic sparsifier $\tilde{H}$.

Similarly to the unweighted case, operations Insert and Delete are reduced to adding terminals to the set $T$ and we refer the reader to the previous section for details on this reduction. Thus, the bulk of our effort is devoted to implementing the procedure AddTerminal. Let $u$ be a non-terminal vertex that we want to append to $T$. We set $T\leftarrow T\cup\{u\}$, and then shorten all the walks at the first location they meet $u$. This shortening of walks induces in turn edge insertions and deletions to $H$, which are then processed by $\tilde{H}$. To quickly locate the first appearances of $u$ in the random walks from $W$, we maintain a linked list $W_{u}$ for each $u\in V$. This linked list contains the first appearances of $w$ in the collections of random walks $W$. Note that constructing such lists can be performed at no additional costs during preprocessing phase, since Algorithm 4 directly gives the first appearances of vertices in every walk belonging to $W$. After locating the first appearances of $u$, we cut the walks in these locations, delete the corresponding affected walks (together with their weight from $H$), and insert the new shorter walks to $H$. Note that we can simply use arrays to represent each random walk in $W$.

We next analyze the performance of our data-structure. Let us start with the preprocessing time. First, by Lemma 4.14 we get that the cost for constructing $H$ on a graph with $m$ edges is bounded by $\tilde{O}(m\beta^{-4}\epsilon^{-4})$. Next, since $H$ has $\tilde{O}(m\epsilon^{-2})$ edges, constructing $\tilde{H}$ takes $\tilde{O}(m\epsilon^{-4})$ time. Thus, the amortized time of Initialize operation is bounded by $\tilde{O}(m\beta^{-4}\epsilon^{-4})$.

We now analyze the update operations. By the above discussion, note that it suffices to bound the time for adding a vertex to $T$, which in turn (asymptotically) bounds the update time for edge insertions and deletions. The main observation we make is that adding a vertex to $T$ only shortens the existing walks, and by the above discussion we can find such walks in time proportional to the amount of edges deleted from the walk. Since the walk needed to be generated in the Initialize operation, the deletion of these edges take equivalent time to generating them. Moreover, we note that (1) handling the updates in $\tilde{H}$ induced by $H$ introduces additional $O(\mathrm{poly}(\log n)\epsilon^{-2})$ overheads, and (2) adding or deleting $\rho$ edges until the next rebuild costs $\tilde{O}(\beta m\epsilon^{-2})$, since we process only up to $\beta m$ operations. These together imply that the amortized cost for adding a terminal can be charged against the preprocessing time, which is bounded by $\tilde{O}(m\beta^{-4}\epsilon^{-4})$, up to poly-logarithmic factors. Thus it follows that the operations AddTerminal, Insert and Delete can be implemented in $\tilde{O}(1)$ amortized update time. ∎

4.4 Dynamic All-Pair Effective Resistance on Weighted Graphs

Following exactly the same arguments as in the proof of Theorem 1.3, we can use the above data-structure to efficiently maintain effective resistances on weighted, undirected dynamic graphs.

Theorem 4.15**.**

For any given error threshold $\epsilon>0$, there is a data-structure for maintaining an weighted, undirected multi-graph $G=(V,E,\boldsymbol{\mathit{w}})$ with up to $m$ edges that supports the following operations in $\tilde{O}(m^{5/6}\epsilon^{-4})$ expected amortized time:

• •

Insert$(u,v,w)$: Insert the edge $(u,v)$ with resistance $1/w$ in $G$.

• •

Delete$(u,v)$: Delete the edge $(u,v)$ from $G$.

• •

EffectiveResistance$(s,t)$: Return a $(1\pm\epsilon)$-approximation to the effective resistance between $s$ and $t$ in the current graph $G$.

Proof.

Let $\mathcal{D}(\tilde{H})$ denote the data structure that maintains a dynamic (sparse) Schur complement $\tilde{H}$ of $G$ (Lemma 4.8). Since $\mathcal{D}(\tilde{H})$ supports only up to $\beta m$ operations, we rebuild $\mathcal{D}(\tilde{H})$ on the current graph $G$ after such many operations. Note that the operations Insert and Delete on $G$ are simply passed to $\mathcal{D}(\tilde{H})$. For processing the query operation $\textsc{EffectiveResistance}(s,t)$, we declare $s$ and $t$ terminals (using the operation AddTerminal of $\mathcal{D}(\tilde{H})$), which ensures that they are both now contained in $\tilde{H}$. Finally, we compute the (approximate) effective resistance between $s$ and $t$ in $\tilde{H}$ using Lemma 2.1.

We now analyze the performance of our data-structure. Recall that the insertion or deletion of an edge in $G$ can be supported in $\tilde{O}(1)$ expected amortized time by $\mathcal{D}(\tilde{H})$. Since our data-structure is rebuilt every $\beta m$ operations, and rebuilding $\mathcal{D}(\tilde{H})$ can be implemented in $\tilde{O}(m\beta^{-4}\epsilon^{-4})$, it follows that the amortized cost per edge insertion or deletion is

[TABLE]

The cost of any $(s,t)$ query is dominated by (1) the cost of declaring $s$ and $t$ terminals and (2) the cost of computing the $(s,t)$ effective resistance to $\epsilon$ accuracy on the graph $\tilde{H}$. Since (1) can be performed in $\tilde{O}(1)$ time, we only need to analyze (2). We do so by first giving a bound on the size of $T$. To this end, note that each of the $m$ edges in the current graph adds two vertices to $T$ with probability $\beta$ independently. By a Chernoff bound, the number of random augmentations added to $T$ is at most $2\beta m$ with high probability. In addition, since $\mathcal{D}(\tilde{H})$ is rebuilt every $\beta m$ operations, the size of $T$ never exceeds $4\beta m$ with high probability. The latter also bounds the size of $\tilde{H}$ by $\tilde{O}(\beta m\epsilon^{-2})$ and gives that the query cost is $\tilde{O}(\beta m\epsilon^{-2})$.

Combining the above bounds on the update and query time, we obtain the following trade-off

[TABLE]

which is minimized when $\beta=m^{-1/6}$, thus giving an expected amortized update and query time of

[TABLE]

5 Dynamic Laplacian Solver in Sub-linear Time

In this section we extend our dynamic approximate Schur complement algorithm to obtain a dynamic Laplacian solver for unweighted, bounded degree graphs. Specifically, as described in Theorem 1.2, our goal is to design a data-structure that maintains a solution to the Laplacian system $\mathbf{L}\mathbf{x}=\mathbf{b}$ under updates to both the underlying graph and the demand vector vector $\mathbf{b}$ while being able to query a few entries of the solution vector. For the sake of exposition, in what follows we assume that the underlying graph is always connected.

Consider the dynamic Schur complement data-structure provided by Lemma 1.1. If the demand vector $\mathbf{b}$ has up to $O(\beta m)$ non-zero entries, for some parameter $\beta\in(0,1)$, we can simply incorporate the vertices corresponding to these entries in the terminal set $T$ using operation AddTerminal of the dynamic Schur complement data-structure (Lemma 1.1). Upon receipt of a query index, we add the corresponding vertex to the Schur complement and (approximately) solve a Laplacian system on the maintained Schur complement. The obtained solution vector can then be lifted back to the Laplacian matrix using the following lemma, which we introduced in the preliminaries.

See 2.6 However, the demand vector $\mathbf{b}$ may have a large number of non-zero entries, thus preventing us from obtaining a sub-linear algorithm with this approach. We alleviate this by projection this demand vector onto the current set of terminals and showing that such a projection can be maintained dynamically while introducing controllable error in the approximation guarantee. At a high level, our solver can be viewed as an one layer version of sparsified block-Cholesky algorithms [40].

We next discuss specific implementation details. Recall that $\mathbf{P}({T})$ is the matrix projection of non-terminal vertices $F$ onto $T$. By Lemma 2.6, it is sufficient to maintain a solution $\mathbf{x}_{T}=\mathbf{SC}(G,T)^{\dagger}\mathbf{P}({T})\mathbf{b}$ dynamically. Since Lemma 1.1 already allows us to maintain a dynamic Schur complement, we need to devise a routine that maintains the projection $\mathbf{P}({T})\mathbf{b}$ of $\mathbf{b}$ under vertex additions to the terminal set.

To this end, we describe an algorithm that maintains such a projection which in turn allows us to again achieve sub-linear running times. The algorithm itself can be viewed as a numerically minded generalization of the approach for the small-support case. Concretely, let $S$ denote the current set of terminals that the algorithm maintains ($S$ and $T$ will always be equal, and we differentiate between them only for the sake of presentation). We initialize $S$ with $O(\beta m)$ vertices from the corresponding entries in $\mathbf{b}$ that have the largest value. Our key structural observation is that if the entries of $\mathbf{b}$ are small, adding vertices to $S$ does not change the projection significantly. To measure the error incurred by declaring some vertex a terminal, we exploit the fact that the projection $\mathbf{P}({S})\mathbf{b}$ itself is tightly connected to specific random walks in the underlying graph. In Subsection 5.3, we show that one can reuse earlier projections, even when new terminals are added to $S$, while paying an error corresponding to the lengths of these random walks and the magnitude of entries in $\mathbf{b}$. We then analyze how to control the accumulation of these errors over a sequence of terminal additions, and also describe an initialization procedure that involves solving a Laplacian system for computing the starting (approximate) projection vector. These together lead to the main lemma of this section, whose implementation details and analysis are deferred to Section 5.1.

Lemma 5.1**.**

Given an error parameter $\epsilon>0$, an unweighted unweighted bounded-degree $G=(V,E)$ with $n$ vertices, a vector $\mathbf{b}\in\mathbb{R}^{n}$ in the image of $\mathbf{L}$, a subset of terminal vertices $S^{\prime}$ and a parameter $\beta\in(0,1)$ such that $|S^{\prime}|=O(\beta m)$, there is a data-structure DynamicProj$(G,S^{\prime},\beta)$ for maintaining a vector $\tilde{\mathbf{b}}$ with $\left\lVert\tilde{\mathbf{b}}-\mathbf{P}({S})\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}\leq\epsilon\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}$ for some $S$ with $S^{\prime}\subseteq S$, $|S|=O(\beta m)$, while supporting at most $\beta^{3}m^{1/2}\epsilon(\mathrm{poly}\log n)^{-1}$ operations in the following running times:

• •

Initialize$(G,S^{\prime},\beta)$: Initialize the data-structure $\tilde{O}(m)$ time.

• •

Insert$(u,v)$: Insert the edge $(u,v)$ to $G$ in $O(1)$ time while keeping $G$ bounded-degree.

• •

Delete$(u,v)$: Delete the edge $(u,v)$ from $G$ in $O(1)$ time.

• •

Change$(u,\mathbf{b}^{\prime}(u),v,\mathbf{b}^{\prime}(v))$: Change $\mathbf{b}(u)$ to $\mathbf{b}^{\prime}(u)$ and changes $\mathbf{b}(v)$ to $\mathbf{b}^{\prime}(v)$ while keeping $\mathbf{b}$ in the range of $\mathbf{L}$ in $O(1)$ time.

• •

AddTerminal$(u)$: Add $u$ to $S$ in $O(1)$ time.

• •

Query$()$: Output the maintained $\tilde{\mathbf{b}}$ in $O(\beta m)$ time.

The following lemma, whose proof will be shortly provided, allows us to combine the approximation guarantees of the data-structures (1) dynamic Schur complement and (2) dynamic Projection.

Lemma 5.2**.**

Let $0<\epsilon\leq\frac{1}{2}$. Let $k$ be a positive number such that $\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}\leq k$. Suppose $\tilde{\mathbf{L}}\approx_{\epsilon}\mathbf{L}$, $\left\lVert\tilde{\mathbf{b}}-\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}\leq\epsilon k$ and $\left\lVert\tilde{\mathbf{x}}-\tilde{\mathbf{L}}^{\dagger}\tilde{\mathbf{b}}\right\rVert_{\tilde{\mathbf{L}}}\leq\epsilon\left\lVert\tilde{\mathbf{L}}^{\dagger}\tilde{\mathbf{b}}\right\rVert_{\tilde{\mathbf{L}}}$. Then $\left\lVert\tilde{\mathbf{x}}-\mathbf{L}^{\dagger}\mathbf{b}\right\rVert_{\mathbf{L}}\leq 10\epsilon k$.

We now have all the necessary tools to present the data-structure for solving Laplacian systems in bounded-degree graphs, which essentially entails combining Lemma 1.1 and Lemma 5.1.

Proof of Theorem 1.2.

Let $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$ denote the data-structure that maintains a dynamic (sparse) Schur complement $\tilde{H}$ of $G$ and an approximate dynamic Projection $\tilde{\mathbf{b}}$ of $\mathbf{P}({S})\mathbf{b}$, respectively. Set $\epsilon\leftarrow(\epsilon/10)$ for both data-structures. Our dynamic solver simultaneously maintains $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$. Since $\mathcal{D}(\tilde{\mathbf{b}})$ supports only up to $\beta^{3}m^{1/2}\epsilon(\mathrm{poly}\log n)^{-1}$, we rebuild both data-structures after such many operations.

We now describe the implementation of the operations. First, we find the first $\beta m$ entries with maximum value in $\mathbf{b}$. We then take the corresponding vertices and initialize $S^{\prime}$ and $T^{\prime}$ to be these $\beta m$ vertices. The implementation of these data-structures involves including the endpoints of each edge with probability $\beta$ to $S$ and $T$, respectively. We couple these algorithms such that $S=T$, and this property will be maintained throughout the algorithm. The operations Insert and Delete on $G$ are simply passed to $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$. The operation Change$(u,\mathbf{b}^{\prime}(u),v,\mathbf{b}^{\prime}(v))$ is passed to $\mathcal{D}(\tilde{\mathbf{b}})$. Upon receipt of a query $\mathbf{x}(u)$, for some vertex $u\in V$, i.e., operation $\textsc{Solve}(u)$, we declare $u$ a terminal (using the operation $\textsc{AddTerminal}(u)$ of both $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$). We then proceed by extracting an approximate Schur complement $\tilde{H}$ of $G$ from $\mathcal{D}(\tilde{H})$ and an approximate projection vector $\tilde{\mathbf{b}}$ from $\mathcal{D}(\tilde{\mathbf{b}})$. Finally, using a black-box Laplacian solver [38], we compute a solution vector $\tilde{\mathbf{x}}_{T}$ to the system $\mathbf{L}_{\tilde{H}}\tilde{\mathbf{x}}_{T}=\tilde{\mathbf{b}}$ and output $\tilde{\mathbf{x}}_{T}(u)$ (this is possible since $u$ was added to $T$).

We next show the correctness of the operation $\textsc{Solve}(u)$. The Laplacian solver guarantees that the vector $\tilde{\mathbf{x}}_{T}$ satisfies

[TABLE]

Data-structure $\mathcal{D}(\tilde{\mathbf{b}})$ guarantees that

[TABLE]

By Lemma 5.8, we know $\left\lVert\mathbf{P}({T})\mathbf{b}\right\rVert_{\mathbf{SC}(G,T)^{\dagger}}\leq\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}$. Bringing together Equations (12) and (13) and applying Lemma 5.2 with $k=\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}$, $\mathbf{L}:=\mathbf{SC}(G,T)$, $\mathbf{b}:=\mathbf{P}({T})\mathbf{b}$, $\mathbf{L}_{\tilde{H}}$ and $\tilde{\mathbf{b}}$ yield

[TABLE]

Using Lemma 2.6 we can lift the vector $\tilde{\mathbf{x}}_{T}$ to a solution $\tilde{\mathbf{x}}$ such that

[TABLE]

Finally, we bound the running time of our dynamic solver. Changes in the demand vector $\mathbf{b}$ can be performed in $O(1)$ times, thus having negligible affect in our running times. The insertion or deletion of an edge in $G$ can be supported in $\tilde{O}(1)$ expected amortized time by both $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$. Since we build our data-structures every $\beta^{3}m^{1/2}\epsilon(\mathrm{poly}\log n)^{-1}$ operations, and the total rebuild cost is dominated by $\tilde{O}(m\beta^{-2}\epsilon^{-4})$, it follows that the amortized cost per edge insertion or deletion is

[TABLE]

The cost of any query is dominated by (1) the cost of declaring the queried vertex $u$ a terminal and (2) the cost of extracting $\tilde{H}$ and $\tilde{\mathbf{b}}$. Since (1) can be performed in $\tilde{O}(1)$ amortized time, we only need to analyze (2). Size of the terminal set $S=T$, which can be easily shown to be $O(\beta m)$ with high probability, immediately implies that the running time for (2) is dominated by $\tilde{O}(\beta m\epsilon^{-2})=\tilde{O}(\beta m\epsilon^{-5})$, which also bounds the query cost.

Combining the above bounds on the query and update time, we obtain the following trade-off

[TABLE]

which is minimized when $\beta=m^{-1/12}$, thus giving an expected amortized update and query time of

[TABLE]

We can replace $m$ by $n$ in the above running time guarantee since by assumption $G$ has bounded-degree throughout the algorithm. ∎

We next prove Lemma 5.2.

Proof of Lemma 5.2.

We will use triangle inequality to decompose the error as:

[TABLE]

and bound each of them separately.

The first term can be bounded by first invoking the similarity of $\mathbf{L}$ and $\tilde{\mathbf{L}}$ to change the norm to $\tilde{\mathbf{L}}$, and applying the guarantees of the solve involving $\tilde{\mathbf{L}}$:

[TABLE]

This norm can in turn be transferred back to $\mathbf{L}$, and the discrepancy between $\mathbf{b}$ and $\tilde{\mathbf{b}}$ absorbed using triangle inequality:

[TABLE] 2. 2.

The second term follows from combining the norms in $\tilde{\mathbf{L}}$ and $\mathbf{L}$ using the approximations between these matrices:

[TABLE]

and once again converting the norm back from $\tilde{\mathbf{L}}$ to $\mathbf{L}$:

[TABLE] 3. 3.

The third term can first be written in terms of the norm of $\mathbf{b}$ against a matrix involving the difference between $\mathbf{L}$ and $\tilde{\mathbf{L}}$:

[TABLE]

where because $\mathbf{L}^{1/2}\left(\tilde{\mathbf{L}}^{\dagger}-\mathbf{L}^{\dagger}\right)\mathbf{L}^{1/2}$ is a symmetric matrix, we have by the definition of operator norm:

[TABLE]

Composing both sides of $\tilde{\mathbf{L}}\approx_{\epsilon}\mathbf{L}$ by $\mathbf{L}^{1/2}$ gives $\mathbf{L}^{1/2}\tilde{\mathbf{L}}\mathbf{L}^{1/2}\approx_{\epsilon}\mathbf{L}^{1/2}\mathbf{L}\mathbf{L}^{1/2}$, or upon rearranging:

[TABLE]

or $\left\|\mathbf{L}^{1/2}(\tilde{\mathbf{L}}^{\dagger}-\mathbf{L}^{\dagger})\mathbf{L}^{1/2}\right\|_{2}^{2}\leq\epsilon$. Substituting this bound into Equation 15 above then gives the result.

Summing up these three cases as in Equation 14 then gives the overall result

[TABLE]

∎

5.1 Dynamic Projection

We next discuss the main ideas behind the dynamic algorithm that maintains an approximate projection in Lemma 5.1 and then formally describe the implementation of this data-structure together with its running time guarantees. To this end, suppose we are given an approximate projection $\tilde{\mathbf{b}}$ of $\mathbf{P}({S})\mathbf{b}$ satisfying the following inequality

[TABLE]

The crucial idea is to exploit the fact that the right hand side of the above inequality $\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}$ corresponds to the square root of the energy need by the electrical flow to route $\mathbf{b}$ (see Lemma 2.1 in [50]). Since we assume that our dynamic graph $G$ has bounded-degree, this energy is lower-bounded by

[TABLE]

Let $S^{\prime}$ be the set of $\beta m$ vertices such that their corresponding coordinates in $\mathbf{b}$ have the largest values. Without loss of generality, scale all the entries in the vector $\mathbf{b}$ such that

[TABLE]

By definition of $S^{\prime}$, after up to $(\beta m)/2$ operations in our data-structure, we know that at least half of the vertices in $S^{\prime}$ will keep their corresponding $\mathbf{b}$ values. Thus the allowable error from right hand side of Equation (16), $\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{{\dagger}}}$, is lower bounded by $\Omega(\sqrt{\beta m})$. Our goal is to control the error between the maintained approximate projection $\tilde{\mathbf{b}}$ and the true projection $\mathbf{P}({S})\mathbf{b}$. Our algorithm has two main components. First, it shows how to use a Laplacian solver that computes an approximate projection $\tilde{\mathbf{b}}$ of $\mathbf{P}({S})\mathbf{b}$ satisfying Equation (16) in nearly-linear time. Second, it gives a way to control the error of the projection $\mathbf{P}({S})\mathbf{b}$ under terminal additions to $S$ with respect to the $\left\lVert\cdot\right\rVert_{\mathbf{L}^{{\dagger}}}$ norm.

The initialization lemma, whose proof is deferred to Subsection 5.2 is given below.

Lemma 5.3**.**

Given an unweighted graph $G=(V,E)$ with $n$ vertices and $m$ edges, a demand vector $\mathbf{b}\in\mathbb{R}^{n}$, set of vertices $S\subseteq V$ and an error parameter $\epsilon>0$, there is an $\tilde{O}(m)$ time algorithm that computes a vector $\tilde{\mathbf{b}}$ such that

[TABLE]

To elaborate on the second component of the algorithm, consider the error induced on $\mathbf{P}({S})\mathbf{b}$ when we add a vertex $u$ to some terminal set $S$

[TABLE]

In other words, the above expression gives the error when we simply keep the same vector $\mathbf{b}$ under a terminal addition to the set $S$. We will show that over a certain number of such addition we can bound the compounded error by $O(\sqrt{\beta m})$. Since the latter is a lower bound on $\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}$, it follows that the maintained projection still provides a good approximation guarantee. The following lemma, whose proof is deferred to Subsection 5.3, bounds the error after one terminal addition.

Lemma 5.4**.**

Consider an unweighted undirected bounded-degree graph $G=(V,E)$, a demand vector $\mathbf{b}\in\mathbb{R}^{n}$ and a parameter $\beta\in(0,1)$. Let $S\subseteq V$ with $|S|=O(\beta m)$ such that $\left|\mathbf{b}(u)\right|\geq 1$ for all $u\in S$, and $\left|\mathbf{b}(u)\right|\leq 1$ for all $u\in V\setminus S$. For each edge in $G$, include its endpoints to $S$ independently, with probability at least $\beta$. Then, for any vertex $u\in V\setminus S$, with high probability

[TABLE]

We now have all the necessary tools to give a dynamic data-structure that maintains an approximate projection, i.e., prove Lemma 5.1.

Proof of Lemma 5.1.

Given the input demand vector $\mathbf{b}$, let $S^{\prime}$ be the set of $\beta m$ vertices such that their corresponding coordinates in $\mathbf{b}$ have the largest values. Without loss of generality, scale $\mathbf{b}$ according to Equation (17). For each edge in $G$ include its endpoints to $S^{\prime}$ independently, with probability at least $\beta$.

We next describe the implementation of the operations. For implementing procedure $\textsc{Initialize}(G,S^{\prime},\beta)$, we invoke Lemma 5.3 with $\epsilon/2$. Let $\tilde{\mathbf{b}}$ be the output approximate projection satisfying Equation (16) with error parameter $\epsilon/2$ and set $S=S^{\prime}$. As we will shortly see, operations Insert and Delete will be reduced to adding terminals to the set $S$. Thus we first discuss the implementation of the operation AddTerminal. To this end, let $u$ be a non-terminal vertex that we want to append to $S$. We set $S=S\cup\{u\}$ and simply add an entry $\tilde{\mathbf{b}}(u)=0$ to $\tilde{\mathbf{b}}$ while keeping the rest of the entries unaffected. To insert or delete an edge from the current graph, we simply run AddTerminal procedure for the edge endpoints.

Consider the operation $\textsc{Change}(u,\mathbf{b}(u)^{\prime},v,\mathbf{b}(v)^{\prime})$. We first invoke AddTerminal on both $u$ and $v$ and then add $\mathbf{b}(u)^{\prime}-\mathbf{b}(u)$ to $\tilde{\mathbf{b}}(u)$ and $\mathbf{b}(v)^{\prime}-\mathbf{b}(v)$ to $\tilde{\mathbf{b}}(v)$. Finally, to implement Query we simply return the approximate projection $\tilde{\mathbf{b}}$.

We next analyze the correctness of our data-structure which solely depends on the correctness of AddTerminal and Change operations. We will show that after $k$ many such operations, our maintained approximate projection $\tilde{\mathbf{b}}$ satisfies

[TABLE]

where $S$ denotes the set of terminals after $k$ operations. Note that when $(k=0)$, the above inequality holds by Lemma 5.3 that implements the initialization. Let us analyze the error when a single terminal is added to $S$, i.e., $(k=1)$. Then Lemma 5.4 implies that

[TABLE]

Combining these two guarantees and applying triangle inequality, we get that the error after one terminal addition is

[TABLE]

Next, to analyze the changes in the values of $\mathbf{b}$, let the updated $\mathbf{b}^{\prime}$ be the updated vector $\mathbf{b}$. Let $\tilde{\mathbf{b}}^{\prime}$ be the updated $\tilde{\mathbf{b}}$ and let $\mathbf{P}({S})\mathbf{b}^{\prime}$ be the updated $\mathbf{P}({S})\mathbf{b}$. Using the fact that $u$ and $v$ are added to $S$, we get that

[TABLE]

which in turn implies that

[TABLE]

and thus the error vector does not change.

We showed that after each operation, either the correct vector moves by at most $\tilde{O}(\beta^{-5/2})$ with respect to its $\left\lVert\cdot\right\rVert_{\mathbf{L}^{{\dagger}}}$ norm, or $\tilde{\mathbf{b}}-\mathbf{P}({S})$ does not changes. Thus repeating the above argument $k$ times yields Equation (18). Setting $k=c_{\textrm{EN}}\cdot m^{1/2}\beta^{3}\epsilon(\mathrm{poly}\log n)^{-1})$ such that $\left\lVert\mathbf{b}\right\rVert_{\mathbf{L}^{\dagger}}=c_{\textrm{EN}}\cdot\sqrt{\beta m}$, we get that

[TABLE]

For the running time, Lemma 5.3 implies that the initialization cost is bounded by $\tilde{O}(m)$. Since the size of the maintained vector $\tilde{\mathbf{b}}$ is bounded by $|S|$, it follows that the query cost is $O(\beta m)$. All the remaining operations can be implemented in $O(1)$ time. ∎

5.2 Initialization of Approximate Projection Vector

In this subsection we show how to compute an initial approximate projection vector of $\mathbf{P}({S})\mathbf{b}$, i.e., we prove Lemma 5.3.

Proof of Lemma 5.3.

Define $F=V\setminus S$ and let $G^{\prime}$ be an $n^{\prime}$-vertex graph obtained from $G$ by contracting all vertices in $S$ within $G$ into a single vertex $s$ and keeping parallel edges. Let $\mathbf{L}^{\prime}$ denote the corresponding Laplacian matrix of $G^{\prime}$ and consider the induced vertex mapping $\gamma:V\rightarrow V(G^{\prime})$ with $\gamma(u)=u$ for $u\in F$ and $\gamma(u)=s$ for $u\in S$. Let $\mathbf{b}^{\prime}\in\mathbb{R}^{n^{\prime}}$ be the corresponding demand vector in $G^{\prime}$ such that for $u\in V$, $\mathbf{b}^{\prime}(\gamma(u))=\mathbf{b}(u)$ if $\gamma(u)=u$ and $\mathbf{b}^{\prime}(\gamma(u))=\sum_{v\in S}\mathbf{b}(v)$ otherwise. For the given error parameter $\epsilon>0$, we can invoke a black-box Laplacian solver to compute an approximate solution vector $\tilde{\mathbf{v}}^{\prime}$ to $\mathbf{v}^{\prime}=\mathbf{L}^{\prime{\dagger}}\mathbf{b}^{\prime}$ such that

[TABLE]

Now, to lift back the vector $\tilde{\mathbf{v}}^{\prime}$ to $G$ we define new vectors $\tilde{\mathbf{v}}$ and $\mathbf{v}$ such that for all $u\in V$

[TABLE]

Observe that for any edge $e=(u,v)$ in $G$, we have that

[TABLE]

The above relations imply that the approximation guarantee from Equation (19) can be written as follows

[TABLE]

It is well known that if we interpret $G$ as a resistor network, $\mathbf{v}$ represents the voltage vector on the vertices induced by the electrical flow that routes a certain demand in the network (see e.g., [17]). Thus, by linearity of electrical flows and our construction, we can view $\mathbf{v}$ as being the sum of the voltage vectors corresponding to the electrical flows that route $\mathbf{b}(u)$ amount of flow to $S$, where sum is taken over all $u\in F$. By Lemma 2.8, for each $u\in F$, the demand corresponding to the electrical flow that send $\mathbf{b}(u)$ units of flow to $S$ is given by

[TABLE]

Summing over all $u\in F$ we get the demand vector corresponding to $\mathbf{v}$

[TABLE]

where $\mathbf{b}|_{U}$ is the restriction of $\mathbf{b}$ on the subset $U$ with $\mathbf{b}|_{U}(u)=\mathbf{b}(u)$ if $u\in U$, and $\mathbf{b}|_{U}(u)=0$ otherwise. Since we determined the demand vector corresponding to $\mathbf{v}$, we get that

[TABLE]

Define the approximate project vector $\tilde{\mathbf{b}}$ that our algorithm outputs using the following relation

[TABLE]

where $\tilde{\mathbf{v}}$ is the extended voltage vector we defined above. To complete the proof of the lemma, it remains to bound the difference between $\tilde{\mathbf{b}}$ and $\mathbf{P}({S})\mathbf{b}$ with respect to the $\mathbf{L}^{\dagger}$ norm. To this end, using Equations (21) and (22) we have

[TABLE]

Using the approximate guarantee in Equation (20) we have that

[TABLE]

where the last inequality follows from the fact that the minimum energy needed to route $\mathbf{b}$ becomes smaller when contracting vertices. ∎

5.3 Stability of Projected Vectors

In this subsection we prove our core structural observation, namely that the projection vectors remain stable under the addition of a new terminal vertex, as stated in Lemma 5.4.

We start by considering the projection vector $\mathbf{P}({S})\boldsymbol{1}_{u}$, where $u\in F=V\setminus S$. Recall that for $s\in S$, Lemma 2.7 gives that $[\mathbf{P}({S})\boldsymbol{1}_{u}](s)$ is the probability that the random walk that starts at $u$ hits the set $S$ at the vertex $s$. Equivalently, we can view the probability of this walk as routing a fraction of $\boldsymbol{1}_{u}$ from $u$ to $s$. Now, consider the operation of adding a non-terminal $u\in F$ to $S$, i.e., $\tilde{S}=S\cup\{u\}$. We observe that the fraction of $\boldsymbol{1}_{u}$ that we routed to some vertex $v$ in $S$ might have used the vertex $u\in F$. This indicates that this this fraction should have stopped at $u$, instead of going to other vertices in $S$, which in turn implies that the old projection vector $\mathbf{P}({S})\boldsymbol{1}_{u}$ is not valid anymore. We will later show that this change is tightly related to the load that random walks from other vertices in $F$ put on the new terminal vertex $u$. In the following we focus on showing a provable bound on this load, which in turn will allow us to control the error for the maintained projection vector.

Concretely, for each vertex $u\in F$, we want to bound the load incurred by the random walks of the other vertices $v\in F\setminus u$ to the set $S$. For the purposes of our proof, it will be useful to introduce some random variable. For $v\in F$, let $Z_{v}(S)$ be the set of vertices visited in a random walk starting at $v$ and ending at some vertex in $S$. For $t\geq 0$, let $X_{v}(t)$ be the set of vertices visited in a random walk starting at $v\in F$ after $t$ steps. For a demand vector $\mathbf{b}$ and any two vertices $u,v\in F$, the contribution of $v$ to the load of $u$, denoted by $Y_{v}(u)$, is defined as follows

[TABLE]

The load of a vertex $u\in F$, denoted by $N_{u}$, is obtained by summing the contributions over all vertices in $F$, i.e.,

[TABLE]

The following lemma gives a bound on the expected load of every non-terminal vertex.

Lemma 5.5**.**

For a parameter $\beta\in(0,1)$ and every vertex $u\in F$ we have that ${\mathbb{E}}\left[N_{u}\right]=\tilde{O}(\beta^{-2})$.

For proving the above lemma it will be useful to rewrite the load quantity. To this end, recall that in the proof of Theorem 4.3 we have shown that any random walk that start at a vertex $v$ of length $\ell=\tilde{O}(\beta^{-2})$ hits a vertex in the terminal set $T$ with probability at least $1-1/n^{c}$, for some large constant $c$. Note that by construction of $S$ in Lemma 5.4, the exact same argument applies to the set $S$. Thus, instead of terminating the random walks once they hit $S$, we can run all the walks from the vertices in $F$ up to $\ell$ steps. The latter together with the assumption $\mathbf{b}(v)\leq 1$ for all $v\in F$ (provided by Lemma 5.4) give that

[TABLE]

It turns out that that the term contained in the brackets of Equation (5.3) equals $\deg(u)$. Formally, we have the following lemma.

Lemma 5.6**.**

Let $G$ be an undirected unweighted graph. For any vertex $u\in V$ and any length $t\geq 0$, we have

[TABLE]

To prove this, we use the reversibility of random walks, along with the fact that the total probability over all edges of a walk starting at $e$ is $1$ at any time. Below we verify this fact in a more principled manner.

Proof of Lemma 5.6.

The proof is by induction on the length of the walks $t$. When $t=0$, we have

[TABLE]

which gives a total of $\deg(u)$.

For the inductive case, assume the result is true for $t-1$. The probability of a walk reaching $u$ after $t$ steps can then be written in terms of its location at time $t-1$, the neighbor $x$ of $u$, as well as the probability of reaching there:

[TABLE]

Substituting this into the summation to get

[TABLE]

which upon rearranging of the two summations gives:

[TABLE]

By the inductive hypothesis, the term contained in the bracket is precisely $\deg(x)$, which cancels with the division, and leaves us with $\deg(u)$. Thus the inductive hypothesis holds for $t$ as well. ∎

Plugging Lemma 5.6 in Equation (5.3), along with the fact that by assumption $G$ has bounded degree we get that

[TABLE]

thus proving Lemma 5.5.

We now have all the tools to prove Lemma 5.4.

Proof of Lemma 5.4.

Recall that $\tilde{S}=S\cup\{u\}$, where $u$ is vertex in $F=V\setminus S$. We want to obtain a bound on the difference $(\mathbf{P}({S})\mathbf{b}-\mathbf{P}({\tilde{S}})\mathbf{b})$ with respect to the $\mathbf{L}^{{\dagger}}$ norm. We distinguish the following types of entries of the difference vector: (1) newly added terminal $u$, (2) the old terminals $S$ and (3) the remaining non-terminal vertices $F\setminus\{u\}$. Note that $\mathbf{P}({S})\mathbf{b}$ and $\mathbf{P}({\tilde{S}})\mathbf{b}$ are not $n$-dimensional vectors, so we assume that all missing entries are appended with zeros. This also makes it possible to compute the $\mathbf{L}^{{\dagger}}$ norm.

In what follows, we will repeatedly make use of the following relation by Lemma 2.7 for vertices $u\in F$ and $v\in S$

[TABLE]

For the type (1) entry, i.e., newly added terminal $u$, using the definition of the load $N_{u}$, we get:

[TABLE]

Note that for type (3) entries, i.e., the remaining non-terminals $v\in F\setminus\{u\}$, we have that

[TABLE]

Finally, for type (2) entries, i.e., old terminals $v\in S$, similarly to the type (1) entries we get

[TABLE]

Bringing together Equations (5.3), (25) and (5.3) we get that

[TABLE]

The right-hand side of the equation can be interpreted as routing ${\mathbb{E}}\left[N_{u}\right]$ unit of flows from $u$ to $S$. Thus, to measure the error, we simply need to upper-bound the square root of the energy need to route ${\mathbb{E}}\left[N_{u}\right]$ amount of flow from $u$ to $S$ (Lemma 2.8),i.e.,

[TABLE]

By the simplifying assumption that $G$ is connected and the fact that each endpoint of an edge in $E$ is added to $S$ independently, with probability at least $\beta$, it is easy to show that with high probability, there exists a path $p(v,S)$ from $u$ to $S$ that uses at most $O(\beta^{-1}\log n)$ edges. Hence, if we route ${\mathbb{E}}\left[N_{u}\right]$ units of flow from $u$ to $S$ along the path $p(v,S)$, the energy of such a flow is upper-bounded by

[TABLE]

Using the latter we get that

[TABLE]

where the last inequality uses the fact that ${\mathbb{E}}\left[N_{u}\right]=\tilde{O}(\beta^{-2})$ by Lemma 5.5. This completes the proof the lemma. ∎

5.4 Maintaining Energy of the Electrical Flow

In this subsection we show how our dynamic solver can be extended to maintain the energy of electrical flows for routing an arbitrary demand vector $\mathbf{b}$. This extension can be thought as a generalization of the All-Pair Effective Resistance problem with $\mathbf{b}=(\boldsymbol{1}_{u}-\boldsymbol{1}_{v})$.

Theorem 5.7**.**

For any given error threshold $m^{-1}<\epsilon<1$, there is a data-structure for maintaining an unweighted, undirected bounded degree multi-graph $G=(V,E)$ with $n$ vertices, $m$ edges and a vector $\mathbf{b}\in\mathbb{R}^{n}$ that supports the following operations in $\tilde{O}(\epsilon^{-25/6}m^{5/6})$ expected amortized time:

• •

Insert$(u,v)$: Insert the edge $(u,v)$ with resistance $1$ in $G$.

• •

Delete$(u,v)$: Delete the edge $(u,v)$ from $G$.

• •

Change$(u,\mathbf{b}^{\prime}_{u},v,\mathbf{b}^{\prime}_{v})$: Change $\mathbf{b}_{u}$ to $\mathbf{b}^{\prime}_{u}$ and $\mathbf{b}_{v}$ to $\mathbf{b}^{\prime}_{v}$ while keeping $\mathbf{b}$ in the range of $\mathbf{L}$.

• •

Energy$()$: Return the energy of the electrical flow routing $\mathbf{b}$ on $G$, up to an $\pm\epsilon$ relative error.

The main idea behind the proof of the above theorem is to express the energy of the electrical flow routing $\mathbf{b}$ as sum of the energy to route $\mathbf{b}_{F}$ to $S$ and the energy to route $\mathbf{P}({S})\mathbf{b}$ inside some arbitrary vertex set $S\subseteq V$, where $F=V\setminus S$ and $\mathbf{P}({S})$ is the projection vector. We formalize this approach in the lemma below.

Lemma 5.8**.**

For any undirected graph $G=(V,E)$, demand vector $\mathbf{b}$ in the image of $\mathbf{L}$ and vertex set $S\subseteq V$,

[TABLE]

Proof.

Since it is easy to verify that

[TABLE]

it follows that

[TABLE]

Multiplying the above equation with $\mathbf{b}^{\top}$ from the left and $\mathbf{b}$ from the right yields

[TABLE]

∎

Let us have a closer look at the equation of Lemma 5.8. In particular, consider the sum on the left-hand side of this equation and its two terms. Observe that the second term can be approximated using a black-box Laplacian solver since we already know how to dynamically maintain a sparsifier $\tilde{H}$ of $\mathbf{SC}(G,S)$ as well as an approximation projection $\tilde{\mathbf{b}}$ of $\mathbf{P}({S})\mathbf{b}$. For the first term, similar to the algorithm for maintaining projections, we employ lazy updates. Concretely, we show that that this term does not change too much when moving a vertex $u$ to $S$, which is equivalent to deleting it from $F$. The latter is formalized in the lemma below.

Lemma 5.9**.**

Consider an unweighted, undirected, bounded-degree graph $G=(V,E)$, a demand vector $\mathbf{b}\in\mathbb{R}^{n}$ and a parameter $\beta\in(0,1)$. Let $S\subseteq V$ with $|S|=O(\beta m)$ such that $\left|\mathbf{b}(u)\right|\geq 1$ for all $u\in S$ and $\mathbf{b}(u)\leq 1$ for all $u\in V\setminus S$. For each edge in $G$, include its endpoints to $S$ with probability at least $\beta$. For any $u\in(F=V\setminus S)$, define $\tilde{F}=F\setminus\{u\}$. Then the following two equations hold with high probability,

[TABLE]

and

[TABLE]

Proof.

We first show Equation (28). To this end, let us start by showing that

[TABLE]

Indeed, the following chain of equations give that

[TABLE]

Using this equality and letting $\tilde{S}=S\cup\{u\}$ we get that

[TABLE]

Using the fact that $\mathbf{P}({S})\mathbf{b}_{S}=\mathbf{b}_{S}$ for any $\mathbf{b}$ and $S\subseteq V$ yields

[TABLE]

where the last equality follows from Lemma 5.4, thus proving Equation (28).

Equation (27) follows from the following chain of inequalities

[TABLE]

Proof of Theorem 5.7.

For the sake of exposition, we rename the energy terms described in Lemma 5.8. Specifically, let $\mathcal{E}(\mathbf{b}):=\left\lVert\mathbf{b}\right\rVert^{2}_{\mathbf{L}^{{\dagger}}}$, $\mathcal{E}(\mathbf{b}_{F}):=\left\lVert\mathbf{b}_{F}\right\rVert^{2}_{\mathbf{L}^{{\dagger}}_{F,F}}$ and $\mathcal{E}(\mathbf{P}({S})\mathbf{b}):=\left\lVert\mathbf{P}({S})\mathbf{b}\right\rVert^{2}_{\mathbf{SC}(G,S)^{{\dagger}}}$.

Let $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$ denote the data-structure that maintains a dynamic (sparse) Schur complement $\tilde{H}$ of $G$ and an approximate dynamic projection $\tilde{b}$ of $\mathbf{P}({S})\mathbf{b}$, respectively. Set $\epsilon\leftarrow(\epsilon/7)$. Similar to the dynamic Laplacain solver, our algorithm simultaneously maintains $\mathcal{D}(\tilde{H})$ and $\mathcal{D}(\tilde{\mathbf{b}})$ and it rebuilds after $\beta^{3}m^{1/2}\epsilon(\mathrm{poly}\log n)^{-1}$ operations.

We now describe the implementation of the operations. The operations Inserte and Delete are implemented similarly to the dynamic solver so we only focus on the differences here. In the initialization step, we calculate an approximation $\tilde{\mathcal{E}}(\mathbf{b}):=\left\lVert\tilde{\mathbf{x}}-\mathbf{L}_{G}^{\dagger}\mathbf{b}\right\rVert_{\mathbf{L}}^{2}$ to the energy $\mathcal{E}(\mathbf{b})$ by using a black-box approximate Laplacian solver. Moreover, using the approximate projection $\tilde{\mathbf{b}}$ and the approximate Schur Complement $\tilde{H}$, we calculate an approximation $\tilde{\mathcal{E}}(\mathbf{P}({S})\mathbf{b})$ to $\mathcal{E}(\mathbf{P}({S})\mathbf{b})$ using black-box Laplacian solver. Finally, we define

[TABLE]

and we will shortly show that this a good approximation to $\tilde{\mathcal{E}}(\mathbf{b}_{F})$. We note that the approximation $\tilde{\mathcal{E}}(\mathbf{b}_{F})$ is remains unchanged until we rebuild the data-structure from scratch.

To handle the Query operation, i.e., output the (approximate) energy needed to route $\mathbf{b}$ in $G$, we calculate $\tilde{\mathcal{E}}(\mathbf{P}({S})\mathbf{b})$ using black-box Laplacian solver as in the initialization step, and output

[TABLE]

We next show the correctness of our algorithm. Consider the errors introduced when approximating $\mathcal{E}(\mathbf{b})$ and $\mathcal{E}(\mathbf{P}({S})\mathbf{b})$ during initialization. For the former we have that

[TABLE]

For the latter, we can similarly show that $\tilde{\mathcal{E}}(\mathbf{P}({S})\mathbf{b})$ approximates $\mathcal{E}(\mathbf{P}({S})\mathbf{b})$ within a $\pm 2\epsilon\mathcal{E}(\mathbf{b})$ absolute error. Combining these two errors gives that $\tilde{\mathcal{E}}(\mathbf{b}_{F})$ approximates $\mathcal{E}(\mathbf{b}_{F})$ within a $\pm 4\epsilon\mathcal{E}(\mathbf{b})$ absolute error.

We next study the correctness of the Query operation. Let us first consider the error incurred by keeping the same $\tilde{\mathcal{E}}(\mathbf{b}_{F})$ until the next rebuild. By Lemma 5.9, Equation 27, the change of $\mathcal{E}(\mathbf{b}_{F})$ is bounded by

[TABLE]

where $\mathcal{E}^{\prime}(\mathbf{b}_{F})$ and $\mathbf{b}^{\prime}$ are the corresponding values before the update and

[TABLE]

Observe that $\mathcal{E}(\mathbf{b})$ is $O(\beta m)$. If $\mathcal{E}_{\max}(\mathbf{b}_{F},\mathbf{b}_{\tilde{F}})$ is $O(\beta m)$, then the above error is bounded by $\tilde{O}((\beta^{-5/2}/\sqrt{\beta m})\cdot\mathcal{E}(\mathbf{b}))$. Otherwise, $\mathcal{E}_{\max}(\mathbf{b}_{F},\mathbf{b}_{\tilde{F}})$ is at least $c\sqrt{\beta m}$ for some constant $c$. Since we only process up to $\beta^{3}m^{1/2}\epsilon(\mathrm{poly}\log n)^{-1}$ operations and $\sqrt{\mathcal{E}(\mathbf{b})}$ can change by at most $\tilde{O}(\beta^{-5/2})$ per update (Lemma 5.9 Equation 28), the current $\sqrt{\mathcal{E}(\mathbf{b})}$ is at least $O(\sqrt{\beta m}-\epsilon\sqrt{\beta m})=\Omega(\sqrt{\beta m})$. So again the above error is bounded by $\tilde{O}((\beta^{-5/2}/\sqrt{\beta m})\cdot\mathcal{E}(\mathbf{b}))$.

As the number of updates until next rebuild is bounded by $\beta^{3}m^{1/2}\epsilon(\mathrm{poly}\log n)^{-1}$, we get the total error incurred by keeping the same $\tilde{\mathcal{E}}(\mathbf{b}_{F})$ is $\tilde{O}(\frac{\beta^{-5/2}}{\sqrt{\beta m}}\mathcal{E}(\mathbf{b})\cdot\epsilon\beta^{3}m^{1/2}(\mathrm{poly}\log n)^{-1}))=\epsilon\mathcal{E}(\mathbf{b})$.

Recall that upon a query our algorithm outputs $\tilde{\mathcal{E}}(\mathbf{P}({S})\mathbf{b})+\tilde{\mathcal{E}}(\mathbf{b}_{F})$. We next show that this value approximates $\mathcal{E}(\mathbf{b})$ within a $\pm\epsilon\mathcal{E}(\mathbf{b})$ absolute error, thus concluding the correctness proof. Indeed, summing over the error of $\tilde{\mathcal{E}}(\mathbf{P}({S})\mathbf{b})$, the initial error of $\tilde{\mathcal{E}}(\mathbf{b}_{F})$ and the error incurred during the updates for $\tilde{\mathcal{E}}(\mathbf{b}_{F})$ yields an approximation with $\pm 7\epsilon\mathcal{E}(\mathbf{b})=\pm\epsilon\mathcal{E}(\mathbf{b})$ absolute error, where the last equality follows by our choice of $\epsilon$.

The time complexities of our update and query operations in the data-structure are similar to the dynamic solver and can be shown to be $\tilde{O}(n^{11/12}\epsilon^{-5})$. ∎

Acknowledgements

We thank Daniel D. Sleator for helpful comments on an earlier version of this paper, and Xiaorui Sun for showing us the tight example for the load of random walks as discussed in Appendix A.

Appendix A Lower Bound on Load of Edge

We show that our bound on the maximum load of a vertex from Lemma 5.5 is close to tight. This bound is crucial for bounding the effect of a single insertion/deletion on the random walks since all walks that interact with that vertex go into the error term.

Claim A.1**.**

For any value $k$, there exists a bounded degree graph such that if we take a random walk starting at every vertex until it reaches $k$ distinct vertex, the maximum load of an edge is $\Omega(k^{2}/\log{k})$.

Note that the bounded degree assumption allow us to ignore the distinction between edge and vertex loads.

We build our hard instance by having a very well connected core graph, and paths leading from it. The main observation is that vertices from these paths need to walk on this core piece in order to reach $k$ distinct vertices. An example of such a graph is shown in Figure 2 and a formal definition is given below.

Definition A.2**.**

We say that a graph $G=(V,E)$ is a path augmented expander if it is defined as follows

Let $G$ be a bounded-degree graph of constant expansion on $k$ vertices, denoted as the core. 2. 2.

Extend $G$ by adding for each core vertex $u$, a path of length $k/10$, denoted as the ‘ray’ of $u$.

For our overall bound on the number of steps that a random walk takes in the core component, we need the following bound on the expected number of steps of a walk within each ray.

Lemma A.3**.**

The probability that a random walk starting from one end of a path reaches $n$ vertices to the right before returning to the end of the path is at most $1/n$.

Proof.

The resistance between the second vertex and the vertex $n$ is $(n-1)$, while the resistance between the vertices $1$ and $2$ is $1$. Since the long term behavior of random walks is akin to electrical flow, we get that such a random walk reaches vertex $n$ before vertex $1$ with probability $1/(1+n-1)=1/n$. ∎

We now have all the tools to prove our lower bound.

Proof of Claim A.1.

Consider the graph as given in Definition A.2. Specifically, consider the random walk starting at each vertex from some ray. Since each ray has at most $k/10$ vertices, at least $k/2$ distinct vertices must come from outside of this ray. We show that with probability at least $1/2$, such a walk must reach $\Omega(k/\log{k})$ core vertices before visiting $k$ distinct vertices.

At each core vertex, Lemma A.3 gives that the probability that the walk hits at least $i$ vertices is

[TABLE]

Thus the expected number of distinct vertices visited after $k$ steps is at most most $k/10$. The later gives that with probability at least $1/2$, such a walk takes at least $\Omega(k/\log{k})$ steps among core vertices. As the core graph is an expander, the number of distinct core vertices visited is at least $\Omega(k/\log{k})$ as well.

Since the random walks from each vertex are independent, applying Chernoff bound gives that with high probability at least $\Omega(k^{2})$ walks spend at least $\Omega(k^{2}/\log{k})$ steps in the core. As there are $k$ vertices in the core, one of them must have load $\Omega(k^{2}/\log{k})$. ∎

Appendix B Proofs about Schur Complement

See 2.6

Proof of Lemma 2.6.

We assume without loss of generality that the underlying graph $G$ is connected. Consider the following extended linear system

[TABLE]

Using the definitions of Schur complement and projection matrix, we can rewrite the above equation as follows:

[TABLE]

Multiplying both sides from the left with

[TABLE]

we get that

[TABLE]

what we wanted to show.

∎

See 2.7

Proof of Lemma 2.7.

First, note that if there is no path from vertices in $T$ to $F=V\setminus T$, then the lemma holds trivially. Thus assume $T$ and $F$ are connected by paths. Next, let

[TABLE]

where $\textbf{D}_{F}$ is the diagonal of $\mathbf{L}_{[F,F]}$ and $\textbf{A}_{F}$ is the negation of the off-diagonal entries, and then expand $\mathbf{L}_{[F,F]}^{-1}$ using the Jacobi series:

[TABLE]

The above series converges due to the fact that $\mathbf{L}_{[F,F]}$ is strictly diagonally dominant. Concretely, the latter implies $(\textbf{A}_{F}\textbf{D}_{F}^{-1})^{\ell}$ tends to zero as $\ell$ tends to infinity. Substituting this in the definition of $\mathbf{P}({T})$ and letting $\boldsymbol{1}_{u}=\begin{bmatrix}\boldsymbol{1}_{u}^{F}&\boldsymbol{1}_{u}^{T}\end{bmatrix}^{\top}$ we get that

[TABLE]

In particular, it follows that for any $v\in T$

[TABLE]

See 2.8

Proof of Lemma 2.8.

Given a valid demand vector $\mathbf{b}$ with $\mathbf{b}^{\top}\boldsymbol{1}=0$, Lemma 2.1 due to Miller and Peng [50] shows that the minimum energy for routing $\mathbf{b}$ is given by $\mathbf{b}^{\top}\mathbf{L}^{\dagger}\mathbf{b}$. Since by construction we have that $[\boldsymbol{1}_{u}-\mathbf{P}({K})\boldsymbol{1}_{u}]^{\top}\boldsymbol{1}=0$, substituting this demand vector in place of $\mathbf{b}$ gives the lemma. ∎

Now we prove Theorem 3.1, which states that sampling random walks generates sparsifiers of Schur complements: See 3.1

Note that this rescaling by $1/\rho\left(\sum_{1\leq i\leq\ell}\frac{1}{\boldsymbol{\mathit{w}}_{w_{i-1}w_{i}}}\right)$ is quite natural: Consider the degenerate case where $T=V$. This routine generates $\rho$ copies of each edge, which then need to be rescaled by $1/\rho$ to ensure approximation to the original graph.

Similar to other randomized graph sparsification algorithms [67, 37, 2, 19, 34], our sampling scheme directly interacts with Chernoff bounds. Our random matrices are ‘groups’ of edges related to random walks starting from the edge $e$. We will utilize Theorem 1.1 due to [71], which we paraphrase in our notion of approximations.

Theorem B.1**.**

Let $\mathbf{X}_{1},\mathbf{X}_{2}\ldots\mathbf{X}_{k}$ be a set of random matrices satisfying the following properties:

Their expected sum is a projection operator onto some subspace, i.e.,

[TABLE] 2. 2.

For each $\mathbf{X}_{i}$, its entire support satisfies:

[TABLE]

Then, with high probability, we have

[TABLE]

Re-normalizations of these bounds similar to the work of [67] give the following graph theoretic interpretation of the theorem above.

Corollary B.2**.**

Let $E_{1}\ldots E_{k}$ be distributions over random edges satisfying the following properties:

Their expectation sums to the graph $G$, i.e.,

[TABLE] 2. 2.

For each $E_{i}$, any edge in its support has low leverage score in $G$, i.e.,

[TABLE]

Then, with high probability, we have

[TABLE]

To fit the sampling scheme outlined in Theorem 3.1 into the requirements of Corollary B.2, we need (1) a specific interpretation of Schur complements in terms of walks, and (2) a bound on the effective resistances between two vertices at a given distance.

Given a walk $w=u_{0},\ldots,u_{\ell}$ of length $\ell$ in $G$ with a subset a vertices $T$, we say that $w$ is a terminal-free walk iff $u_{0},u_{\ell}\in T$ and $u_{1},\ldots,u_{\ell-1}\in V\setminus T$.

Fact B.3** ([19], Lemma 5.4).**

For any undirected, weighted graph $G$ and any subset of vertices $T\subseteq V$, the Schur complement $\mathbf{SC}(G,T)$ is given as an union over all multi-edges corresponding to terminal-free walks $u_{0},\ldots,u_{\ell}$ with weight

[TABLE]

The fact below follows by repeatedly applying the triangle inequality of the effective resistances between two vertices.

Fact B.4**.**

In an weighted undirected graph $G$, the effective resistance between two vertices that are connected by a path $p=(p_{0},\ldots,p_{\ell})$ is at most $\sum_{i=0}^{\ell-1}1/\boldsymbol{\mathit{w}}_{p_{i}p_{i+1}}$.

Combining the above results gives the guarantees of our sparsification routine.

Proof of Theorem 3.1.

For every edge $e\in E$, let $W_{e}$ be the random graph corresponding the the terminal-free random walk that started at edge $e$. Define $H=\rho\cdot\sum_{e}W_{e}$ to be the output graph by our sparsification routine, where $\rho=O(\log n\epsilon^{-2})$ is the sampling overhead. To prove that $\mathbf{L}_{H}\approx_{\epsilon}\mathbf{SC}(G,T)$ with high probability, we need to show that (1) ${\mathbb{E}}\left[H\right]=\mathbf{SC}(G,T)$ and (2) for any edge $f$ in $W_{e}$, its leverage score $\boldsymbol{\mathit{w}}_{f}R_{\textrm{eff}}^{W_{e}}(f)$ is at most $\leq\epsilon^{2}/\log n$ (by Corollary B.2). Note that (2) immediately follows from the effective resistance bound of Fact B.4 and the choice of $\rho=O(\log n/\epsilon^{2})$. We next show (1).

To this end, we start by describing the decomposition of $\mathbf{SC}(G,T)$ into random multi-edges, which correspond to random terminal-free walks in Fact B.3. The main idea is to sub-divide each walk $u_{0}\ldots u_{\ell}$ of length $\ell$ in $G$ into $\ell$ walks of the same length, each starting at one of the $\ell$ edges on the walk, and each having weight

[TABLE]

By construction of our sparsification routine, note that every random graph $W_{e}$ is a distribution over walks $u_{0}\ldots u_{\ell}$, each picked with probability

[TABLE]

Thus, to retain expectation, when such a walk is picked, our routine correctly adds it to $H$ with weight $1/(\rho\sum_{i=0}^{\ell-1}1/\boldsymbol{\mathit{w}}_{u_{i}u_{i+1}}).$

Formally, we get the following chain of equalities

[TABLE]

∎

Appendix C Sampling Weights of a Random Walk

In this section, we show that given a random walk $w$ of length $\ell$ in a weighted $G$ with polynomially bounded weights, we can efficiently sample an approximation to $s(w)=\sum_{i=1}^{\ell}(1/\boldsymbol{\mathit{w}}_{w_{i-1}w_{i}})$. Concretely, we prove the following lemma from Section 4.3.

See 4.13

To prove the above lemma, we employ a doubling technique. Specifically, for any pair of vertices $u,v\in V$, and a random walk $w$ of length $\ell$ that starts at $u$ and ends at $v$, it is easy to see that

[TABLE]

where $*$ denotes the convolution between two probability mass functions. However, one challenge here is that we cannot afford dealing with exact representations of probability mass functions as this would be computationally expensive. Instead, we introduce an approximate representation of such functions, and then give an algorithm that allows computing the convolution between such approximate representations. Before proceeding further, note that we can scale down the edge weights so that $\boldsymbol{\mathit{w}}_{e}\leq 1$, and thus $1/\boldsymbol{\mathit{w}}_{e}\geq 1$ for every $e\in E$. In addition, we remark that $\boldsymbol{\mathit{w}}$ does not need to be integral.

Let us introduce a compact way to represent any given probability mass function approximately $f$. The main idea is to ‘move’ each number in the support of $f$ by $(1+\epsilon)$, which in turn results in a $(1+\epsilon)$ approximation of the sampled value for $f$. Formally, let $f$ be a probability mass function such that $f(x)=0$, for each $x\not\in\{0,\ldots,n^{c}\}$, where $c$ is a positive constant. For $j\geq 1$, let $I_{k}^{j}$ be the interval $[(1+\epsilon)^{k},(1+\epsilon)^{k+j})$ for $k\in\{0,\ldots,L\}$ where $L=O\left((c+d)\epsilon^{-1}\log n\right)$. Note that the upper bound $L$ is chosen in such a way that $\cup_{k}I_{k}^{1}$ covers the range of $f_{s(w),\ell}^{u,\ell}$ for every possible triplet $(u,v,\ell)$. For $j\geq 1$ and $\epsilon>0$, we say that $g$ is an $(\epsilon,j)$-approximation of a probability mass function $f$ iff there exists a matrix $\boldsymbol{\mathit{H}}$ satisfying the following properties:

1. (a)

$\sum_{k=0}^{L}\boldsymbol{\mathit{H}}_{x,k}=f(x),\leavevmode\nobreak\ \forall x\in\{0,\ldots,n^{c}\}$, 2. (b)

$\sum_{x=0}^{n^{c}}\boldsymbol{\mathit{H}}_{x,k}=g(k),\leavevmode\nobreak\ \forall k\in\{0,\ldots,L\}$, 3. (c)

$\boldsymbol{\mathit{H}}_{x,k}=0,\leavevmode\nobreak\ \forall x\not\in I_{k}^{j}$.

Note that an $(\epsilon,j)$-approximation of $f$ is also an $\left(\epsilon,(j+1)\right)$-approximation of $f$. Moreover, observe that the intervals $\{I_{k}^{1}\}_{k\in\{0\}\cup L}$ are disjoint for different $k$ but $I_{k}^{j}$ overlaps with $I_{k^{\prime}}^{j}$ whenever $j\geq 2$ and $|k-k^{\prime}|<j$.

Next we show how to compute the convolution of two probability mass functions under their approximate representations. Let and $g^{(1)}$ and $g^{(2)}$ be ($\epsilon,j$)-approximations of probability mass functions $f^{(1)}$ and $f^{(2)}$, respectively. Now consider two intervals $I_{k_{1}}^{j}$ and $I_{k_{2}}^{j}$. Without loss of generality, assume that $k_{1}\leq k_{2}$. If $x\in I_{k_{1}}^{j}$ and $y\in I_{k_{2}}^{j}$, then

[TABLE]

Furthermore, let $I_{k_{3}}^{1}$ be an interval such that $\textrm{le}\in I_{k_{3}}^{1}$. The latter implies that $(1+\epsilon)^{k_{3}}\leq\textrm{le}<(1+\epsilon)^{k_{3}+1}$. Since $\textrm{ri}=\textrm{le}\cdot(1+\epsilon)^{j}$, it follows that $\textrm{ri}<(1+\epsilon)^{k_{3}+j+1}$. Bringing together the above bounds we get that $(1+\epsilon)^{k_{3}}\leq\text{le}<(1+\epsilon)^{k_{3}+j+1}$, i.e., $I^{\prime}\subseteq I_{k_{3}}^{j+1}$. Since $k_{3}$ depends on $k_{1},k_{2},$ and $j$ we sometimes write $k_{3}(k_{1},k_{2},j)$ instead of $k_{3}$.

Since the above approach gives us a way to combine two different intervals, it is now straightforward to compute the convolution between two probability mass functions. This task is performed in the standard way and we review its implementation details in Algorithm 6 for the sake of completeness.

Lemma C.1**.**

Let $j\geq 1$ and $\epsilon>0$ by two parameters. Given any two $(\epsilon,j)$-approximations $g^{(1)}$ and $g^{(2)}$ of probability mass functions $f^{(1)}$ and $f^{(2)},$ Convolute$(g^{(1)},g^{(2)},\epsilon,j)$ (Algorithm 6) computes in $\tilde{O}(\epsilon^{2})$ time an $\left(\epsilon,(j+1)\right)$-approximation $g:=g^{(1)}*g^{(2)}$ of the convolution $f:=f^{(1)}*f^{(2)}$.

Proof.

We first show the correctness. Since $g^{(1)}$ and $g^{(2)}$ are $(\epsilon,j)$-approximations to $f^{(1)}$ and $f^{(2)}$ by assumption of the lemma, we know that there exists matrices $\boldsymbol{\mathit{H}}^{(1)}$ and $\boldsymbol{\mathit{H}}^{(2)}$ satisfying properties (a), (b) and (c). To show that the output $g$ is correct we need to construct a matrix $\boldsymbol{\mathit{H}}$ that satisfies each of these properties. By construction of the algorithm, the new matrix $\boldsymbol{\mathit{H}}$ is defined as follows:

[TABLE]

We start by showing property (a) for $\boldsymbol{\mathit{H}}$. Concretely, for any $z\in\{0,\ldots,n^{c}\}$ we get that

[TABLE]

Next, $\boldsymbol{\mathit{H}}$ satisfies property (b) since for any $k_{3}\in\{0,\ldots,L\}$ we get that

[TABLE]

where the penultimate equality follows by Algorithm 6.

Finally, for every $x\not\in I_{k}^{j}$, we have that $\boldsymbol{\mathit{H}}_{x,k}=0$, i.e., property (c) holds for $\boldsymbol{\mathit{H}}$. The latter holds since $x\in I_{k_{1}}^{j}$ and $y\in I_{k_{2}}^{j}$ gives that $x+y\in I_{k_{3}(k_{1},k_{2},j)}^{j+1}$. Thus, by definition of approximate probability mass function, it follows that $g=g^{(1)}*g^{(2)}$ is an $\left(\epsilon,(j+1)\right)$-approximation of $f=f^{(1)}*f^{(2)}$.

For the running time first recall that $L=O\left((c+d)\epsilon^{-1}\log n\right)=\tilde{O}(\epsilon^{-1})$. Since the cost for implementing Convolute is bounded by $\tilde{O}(L^{2})$, it follows that we can implement this procedure in $\tilde{O}(\epsilon^{-2})$ time. ∎

The last ingredient we need is to show that given a family of probability mass functions, and their corresponding approximations, choosing one of these functions according to some probability distribution yields a random approximation in the natural way. Specifically, for an index set $Q$, let $\{f^{(q)}\}_{q\in Q}$ be a set of probability mass functions. Let $\hat{q}$ be be a random variable (independent from $\{f^{(q)}\}_{q\in Q}$) such that for every $q\in Q$, $\displaystyle\mbox{Pr}[\hat{q}=q]=\boldsymbol{\mathit{p}}(q)$, and $\sum_{q\in Q}\boldsymbol{\mathit{p}}(q)=1$. Furthermore, define

[TABLE]

.

Lemma C.2**.**

Suppose $g^{(q)}$ is an $(\epsilon,j)$-approximation of the probability mass function $f^{q}$, for all $q\in Q$. Let $f$ be the probability mass function as defined above. Then

[TABLE]

is an $(\epsilon,j)$-approximation of $f$.

Proof.

By definition of an $(\epsilon,j)$-approximation, we know that there exist matrices $\{\boldsymbol{\mathit{H}}^{(q)}\}_{q\in Q}$ for $\{g^{(q)}\}_{q\in Q}$ satisfying properties (a), (b) and  (c). We need to show that for $g$ as defined in the lemma, there exist a suitable matrix $\boldsymbol{\mathit{H}}$ that satisifes each of these properties. To this end, define $\boldsymbol{\mathit{H}}$ as follows

[TABLE]

We start by showing property (a). Concretely, for any $z\in\{0,\ldots,n^{c}\}$ we get that

[TABLE]

Next, $\boldsymbol{\mathit{H}}$ satisfies property (b) since for any $k\in\{0,\ldots,L\}$ we get that

[TABLE]

Finally, for every $x\not\in I_{k}^{j}$, we have that $\boldsymbol{\mathit{H}}_{x,k}=0$, i.e., property (c) is satisfied for $\boldsymbol{\mathit{H}}$. The latter holds since for all $x\not\in I_{k}^{j}$ we have that $\boldsymbol{\mathit{H}}^{(q)}_{x,k}=0$ and thus $\boldsymbol{\mathit{H}}_{x,k}=\sum_{q\in Q}\boldsymbol{\mathit{p}}(q)\boldsymbol{\mathit{H}}^{(q)}_{x,k}=0$. As a result we conclude that $g$ is an $(\epsilon,j)$-approximation of $f$ with matrix $\boldsymbol{\mathit{H}}$ satisfying all the required properties. ∎

We now describe how to compute a probability distribution that will in turn allow us to sample approximately from $f_{s(w),\ell}^{u,v}$. At a high level we accomplish this task by employing the “doubling technique” together with the approximate representations and their convolution. As an input, the algorithm receives a weighted graph $G$ with polynomially bounded weights, a length parameter $\ell\geq 1$, an error parameter $\epsilon>0$ and two vertices $u,v\in V$. The procedure computes and outputs a vector $\left(j_{u,v},g_{\ell}^{u,v},p_{\ell}^{u,v}\right)$, where $j_{u,v}\geq 1$ is a precision parameter, $g_{\ell}^{u,v}$ is an $(\epsilon,j_{u,v})$-approximation of $f_{s(w),\ell}^{u,v}$, and $p_{\ell}^{u,v}=\mathbb{P}_{u}\left[w_{\ell}=v\right]$ is the probability that the random walk $w$ that originates at $u$ hits $v$ after $\ell$ steps.

If $(\ell=1)$, then there are two possibilities depending on whether $(u,v)\in E$ or not. If the former holds, then the algorithms simply returns $(1,\vec{0},0)$ as it is not possible to reach $v$ after performing one step of the random walk from $v$. Otherwise, we simply return $(1,g_{1}^{u,v},p_{\ell}^{u,v})$, where $g_{1}^{u,v}(\frac{1}{\boldsymbol{\mathit{w}}_{u,v}})=1$ and $p_{\ell}^{u,v}=\frac{\boldsymbol{\mathit{w}}_{u,v}}{\boldsymbol{\mathit{d}}_{G}(v)}$.

However, if $(\ell>1)$, then it first halves $\ell$ into two parts $\ell^{\prime}=\lfloor\ell/2\rfloor$ and $\ell^{\prime\prime}=\lceil\ell/2\rceil$. Next, for each $y\in V$ it recursively calls itself with input parameters $(G,u,y,\ell^{\prime},\epsilon)$ and $(G,y,v,\ell^{\prime\prime},\epsilon)$. The outputs from these two calls are then combined using the convolution manipulations described above to produce the final output. Exact details for implementing this procedure are summarized in Algorithm 7. The following lemma proves the correctness and the running time of the algorithm.

Lemma C.3**.**

Given a weighted graph $G=(V,E,\boldsymbol{\mathit{w}})$ with $\boldsymbol{\mathit{w}}_{e}\in[1,n^{c}]$ for each $e\in E$ and $c>0$, two vertices $u,v\in V$, a length parameter $\ell\in[1,n^{d}]$ and an error parameter $\epsilon>0$, ComputeDistrib$(G,u,v,\ell,\epsilon)$ (Algorithm 7) correctly computes a vector $(j_{u,v},g_{\ell}^{u,v},p_{\ell}^{u,v})$ in $\tilde{O}(n^{3}\epsilon^{-2})$ time, where $g_{\ell}^{u,v}$ is an $(\epsilon,j_{u,v})$-approximation to $f_{s(w),\ell}^{u,v}$ and $p_{\ell}^{u,v}$ is the probability that the random walk $w$ the starts at $u$ hits $v$ after $\ell$ steps. Moreover, the output $j_{u,v}$ cannot exceed $O(\log n)$.

Proof.

We first prove that the third coordinate of the output vector equals $\mathbb{P}_{u}\left[w_{\ell}=v\right]$. We proceed by induction on the length of the walk $\ell$. If $(\ell=1)$, it is easy to check that the condition holds by construction of the algorithm. Next assume $(\ell\geq 2)$ and note that $(\ell^{\prime}<\ell)$ and $(\ell^{\prime\prime}<\ell)$. Applying induction hypothesis on each recursion call, we know $p^{u,y}_{\ell^{\prime}}$ is $\mathbb{P}_{u}\left[w_{\ell}=y\right]$ and $p^{y,v}_{\ell^{\prime\prime}}$ is $\mathbb{P}_{y}\left[w_{\ell^{\prime\prime}}=u\right]$. The latter along with the fact that $(\ell^{\prime}+\ell^{\prime\prime}=\ell)$ imply

[TABLE]

We next prove that the second coordinate $g_{\ell}^{u,v}$ is an $(\epsilon,j)$-approximation of $f_{s(w),\ell}^{u,v}$. First, since $j_{u,v}=(\max_{y}\max(j_{v,y},j_{y,u}))+1$, Lemma C.1 implies that $\left(g_{\ell^{\prime}}^{u,y}*g_{\ell^{\prime\prime}}^{y,v}\right)$ is an $(\epsilon,j_{u,v})$-approximation of $f^{u,y,v}_{s(w),\ell^{\prime},\ell^{\prime}}$ where we define $f^{u,y,v}_{s(w),\ell^{\prime},\ell^{\prime}}$ to be the probability mass function of $s(w)$, conditioning on $w\sim w_{v,u}$, $\ell(w)=\ell^{\prime}+\ell^{\prime\prime}$ and $w_{\ell^{\prime}}=y$. Second, consider the triplets $\{(u,y,v)\}_{y\in V}$, and let $g_{y}=g_{\ell^{\prime}}^{u,y}*g_{\ell^{\prime\prime}}^{y,v}$ and $p_{y}=p_{\ell^{\prime}}^{u,y}\cdot p_{\ell^{\prime}}{y,v}$. Then by Lemma C.2 we get that $g_{\ell}^{u,v}=\sum_{y\in V}p_{y}\cdot g_{y}$ is the desired $(\epsilon,j_{u,v})$-approximation.

Finally, we prove that $j_{u,v}=O(\log n)$. We will inductively show that the first coordinate $j_{u,v}$ of the output vector from ComputeDistrib$(G,u,v,\ell,\epsilon)$ is at most $k+1$, for $\ell\leq 2^{k}$. For the base case $k=0$, which implies that $\ell=1$ and and the claim trivially holds. Now assume $thatk\geq 1$. Since $\ell\leq 2^{k}$, by construction we get that $\ell^{\prime}\leq 2^{k-1}$ and $\ell^{\prime\prime}\leq 2^{k-1}$. By induction hypothesis, the first coordinates returned by all of the recursion calls are no more than $(k-1)+1=k$. Thus, the returned $j_{u,v}$ at most $k+1=O(\log n)$.

For the running time, note that in all recursion calls of the procedure ComputeDistrib there are at most $n^{2}$ possible pairs $(u,v)$ and $O(\log n)$ possible values of $\ell$. In each of these calls, we invoke the procedure Convolute exactly $n$ times, where each invocation costs $\tilde{O}(\epsilon^{-2})$ by Lemma C.1. Thus the total running time is bounded by $\tilde{O}(n^{3}\epsilon^{-2})$.

∎

We now have all the tools to prove Lemma 4.13.

Proof of Lemma 4.13.

Our algorithm for sampling $s(w)$ is implemented as follows. First, it invokes the procedure $\textsc{computeDistrib}(G,u,v,\ell,\epsilon)$ and obtains the resulting vector $(j_{u,v},g_{\ell}^{u,v},p^{u,v}_{\ell})$. Then it samples from the distribution by choosing the interval $I^{j_{u,v}}_{k}=[\textrm{le},\textrm{ri}]$ with probability $g_{\ell}^{u,v}(k)$, where $\textrm{le}:=(1+O(\frac{\epsilon}{\log n}))^{k}$ and $\textrm{ri}:=(1+O(\frac{\epsilon}{\log n}))^{k+j_{u,v}}$. Finally the algorithm outputs ri. This procedure is summarized in Algorithm 8.

We next argue about the correctness. Note that by property (b) in the definition of approximation $g_{\ell}^{u,v}$ of $f_{s(w),\ell}^{u,v}$, this sampling process can be viewed as sampling the pair $(x,i)$ from the distribution $\boldsymbol{\mathit{H}}_{x,i}$, without knowing $x$. Furthermore, by property (a), each $x$ is sampled with the correct probability $f_{s(w),\ell}^{u,v}(x)$. Since $\epsilon\leq 0.5$ it follows by Lemma C.3 that $\textrm{ri}/\textrm{le}=(1+O(\frac{\epsilon}{\log n}))^{j_{u,v}}\leq(1+\epsilon)$. Thus by property (c) we get that ri is within $[x,(1+\epsilon)x]$ for the (unknown) sampled $x$.

The running time of our sampling procedure is asymptotically dominated by the running time of ComputeDistrib, which is in turn bounded by $O(n^{3}\epsilon^{-2})$, as desired. ∎