Kirchhoff Index As a Measure of Edge Centrality in Weighted Networks: Nearly Linear Time Algorithms
KKirchhoff Index As a Measure of Edge Centrality in WeightedNetworks: Nearly Linear Time Algorithms
Huan Li
School of Computer ScienceFudan University [email protected]
Zhongzhi Zhang
School of Computer ScienceFudan University [email protected]
Abstract
Estimating the relative importance of vertices and edges is a fundamental issue in theanalysis of complex networks, and has found vast applications in various aspects, such as socialnetworks, power grids, and biological networks. Most previous work focuses on metrics of verteximportance and methods for identifying powerful vertices, while related work for edges is muchlesser, especially for weighted networks, due to the computational challenge. In this paper, wepropose to use the well-known Kirchhoff index as the measure of edge centrality in weightednetworks, called θ -Kirchhoff edge centrality. The Kirchhoff index of a network is defined as thesum of effective resistances over all vertex pairs. The centrality of an edge e is reflected in theincrease of Kirchhoff index of the network when the edge e is partially deactivated, characterizedby a parameter θ . We define two equivalent measures for θ -Kirchhoff edge centrality. Both areglobal metrics and have a better discriminating power than commonly used measures, based onlocal or partial structural information of networks, e.g. edge betweenness and spanning edgecentrality.Despite the strong advantages of Kirchhoff index as a centrality measure and its wideapplications, computing the exact value of Kirchhoff edge centrality for each edge in a graph iscomputationally demanding. To solve this problem, for each of the θ -Kirchhoff edge centralitymetrics, we present an efficient algorithm to compute its (cid:15) -approximation for all the m edges innearly linear time in m . The proposed θ -Kirchhoff edge centrality is the first global metric ofedge importance that can be provably approximated in nearly-linear time. Moreover, accordingto the θ -Kirchhoff edge centrality, we present a θ -Kirchhoff vertex centrality measure, as wellas a fast algorithm that can compute (cid:15) -approximate Kirchhoff vertex centrality for all the n vertices in nearly linear time in m . a r X i v : . [ c s . D S ] J a n Introduction
Most real networks (e.g. social networks) are massive and inhomogeneous [New10], where theroles of vertices/edges are often largely different, with peripheral vertices/edges having a limitedeffect on their function, while central vertices/edges having a strong impact on dynamical processes.Thus, it is of paramount importance to design both desirable metrics measuring the centrality orimportance of vertices/edges and fast algorithms identifying vital vertices/edges [LM12]. In pastdecades, a lot of centrality measures have been presented to capture diverse aspects of the informalconcept of importance, and various algorithms for these different metrics have been developedby researchers from interdisciplinary areas, such as computer science [WS03, BV14, BDFMR16],control science [YTQ17], and physics [LCR + + +
13, MGLKT15, HAY16], and current-flow centrality [BF05]. The betweenness of an edge is the fraction of shortest paths between vertexpairs that pass through the edge. The spanning edge centrality of an edge is defined as the probabilitythat it is present in a randomly chosen spanning tree. While the current-flow centrality of an edgedescribes the amount of current flowing through it. Although these edge centrality metrics have beenextensively studied, they themselves are subject to weakness. For example, edge betweenness onlyconsiders shortest paths and ignores those longer paths; spanning edge centrality cannot separatean edge linked to a leaf vertex and another cut edge connecting two large subgraphs, while theimportance of these two edges are obviously different. Moreover, these measures are proposed forunweighted networks, and are either unapplicable to weighted networks, or have high computationalcomplexity when applied to weighted networks.In fact, it is very difficult to rigorously compare different measures of edge centrality, sincethe criteria of edge importance depend on real applications and the problems we are concernedwith [YTQ17]. Hence, it is neither practical nor feasible to propose a universal measure that bestquantifies the importance of edges for all situations. One should define the metric of edge centralityaccording to particular problems. In many real scenarios, the Kirchhoff index [KR93] of a network,defined as the sum of effective resistances over all vertex pairs, can be used as a unifying indicatorto measure the interesting quantities associated with different problems in networks. For example,Kirchhoff index can be used to measure the mean cost of search in a complex network [FQY14],robustness of first-order consensus algorithm in noisy networks [PB14], the global utility of socialrecommender systems [WLC16], among others. Notwithstanding the relevance of the Kirchhoffindex in various applications, there is a disconnect between this notion and efficiency of algorithmsfor estimating it.The main purpose of this paper is to develop an edge centrality notion that not only has gooddiscriminating power, but also can be evaluated using algorithms with good provable performances.For a connected graph, the popular Kirchhoff index follows Rayleigh’s monotonicity law [ESVM + e , we partially deactivate theedge e by changing its weight w ( e ) to θw ( e ), where 0 < θ ≤ / e is reflected in the Kirchhoffindex of the new graph: the larger the Kirchhoff index is, the more important the edge e is. Wedefine two equivalent metrics for edge centrality. One is the Kirchhoff index of the new graph afteredge deactivation, the other is the difference of the Kirchhoff indices between the new graph andthe original graph. For either edge centrality measure, we give a fast algorithm to compute the (cid:15) -approximation for all the m edges in nearly linear time. Furthermore, based on the Kirchhoff edgecentrality index, we propose a vertex importance measure, with the centrality of a vertex beingdefined as the Kirchhoff index of a new graph, where all edges incident to it are deactivated, andprovide an efficient algorithm for estimating this new vertex centrality. Some edge centrality measures and related algorithms have been proposed. Here we give a briefintroduction to these metrics and their computational complexity. Moreover, we simply describesome work or techniques that partially motivate this paper or relate to our algorithm.Edge betweenness is probably the most popular and most studied measure of edge importance.It measures the probablility that a shortest path between two vertices passes through a given edge.A fast algorithm for exact computation of edge betweenness was developed by Brandes [Bra01].For a graph with n vertices and m edge, the complexity for this efficient technique is O ( nm ) and O ( nm + n log n ) for unweighted graphs and weighted graphs, respectively. In order to speed upthe computation, some approximate algorithms have been proposed [BKMM07, BP07, GSS08]. Allthese approximate approaches aim at reducing the computation of shortest paths in different ways,without providing approximation guarantees.Another edge importance measure is spanning edge centrality first introduced in [TMC + O ( mn / ). In order to computespanning edge centrality for massive networks, two fast approximation algorithms [MGLKT15,HAY16] have been designed, both having theoretical guarantees on their accuracy.A third measure for edge importance is current-flow centrality introduced by Brandes andFleischer [BF05]. An edge has relatively significant importance, if it participates in many shortpaths connecting pairs of vertices. Brandes and Fleischer [BF05] provided an algorithm with timecomplexity O ( mn / log n ), which can actually be dropped to O ( mn log n ) as shown in [MGLKT15].Both spanning edge centrality and current-flow centrality are closely related to effective resis-tance [MGLKT15]. In fact, the Kirchhoff edge centrality we propose belongs to the same class of electrica l centrality measures. Moreover, this is the first definition of a global notion of centralitythat can be provably approximated in nearly-linear time, which means that resistance based edgecentrality for graphs may actually be easier to compute than other centrality measures based ondiscrete structures, e.g. triangles or shortest paths.As many previous theoretical studies [BHNT15, KP17], our work is also motivated by graphmining applications. In [BHNT15], a streaming algorithm was developed for analyzing large-scalerapidly-changing graphs, which maintains densest subgraphs in one pass and achieves time andspace efficiency, whereas in our case, effective resistances are persevered under updates. All thenotions (dense subgraph [BHNT15], triangle [KP17], and Kirchhoff centrality) studied before or3n the present paper have vast applications in network analysis, and their related computationalchallenges fall within the scope of computer theory.Our algorithms, in particular the resistance maintenance routines, closely buid upon the sketchingbased inverse maintenance routine from [LSW15] and the computation of multiple partial states ofGaussian eliminations from [DKP + + For a graph G , we write G \ θ e to denote the graph obtained from G by deactivating edge e , i.e.,decreasing the weight of e from w ( e ) to θw ( e ) for some small 0 < θ ≤ /
2. Let L be the Laplacianmatrix of G , and let L \ θ e denote the Laplacian matrix of G \ θ e . Then we can define two metrics for θ -Kirchhoff centrality of an edge e , denoted by C θ ( e ) and C ∆ θ ( e ), respectively. C θ ( e ) is the Kirchhoffindex of the graph G \ θ e , i.e., the sum of effective resistances over all vertex pairs in G \ θ e , while C ∆ θ ( e ) is the difference between the Kirchhoff indices of graph G \ θ e and graph G . Let K ( G ) denotethe Kirchhoff index of graph G , then we have C θ ( e ) = K ( G \ θ e ) and C ∆ θ ( e ) = K ( G \ θ e ) − K ( G ).As has been shown in [ESVM + n times Tr (cid:16) L † (cid:17) , where L † is the pseudoinverse of the graph’s Laplacian matrix. Thus, we have C θ ( e ) = n Tr (cid:16) ( L \ θ e ) † (cid:17) and C ∆ θ ( e ) = n Tr (cid:16) ( L \ θ e ) † (cid:17) − n Tr (cid:16) L † (cid:17) . To compute the exact value of θ -Kirchhoff centrality for eachedge, a naive algorithm would invert the matrix L \ θ e for all e ∈ E . Since a single inversion takes O ( n ω ) time, where ω ≈ .
373 is the matrix multiplication constant [Wil12], the naive algorithmruns in O ( n ω m ) time for all the m edges, which makes it untractable for large networks.In this paper, we consider the scenario in which only approximate values of θ -Kirchhoff centralityare needed. Such approximations are acceptable in many cases because we only need to estimaterelative importance of edges. We give a randomized algorithm EdgeCentComp1 that computes (cid:15) -approximate Kirchhoff edge centrality C θ ( e ) for all the m edges in ˜ O ( m(cid:15) − ) time, and a randomizedalgorithm EdgeCentComp2 that computes (cid:15) -approximate Kirchhoff edge centrality C ∆ θ ( e ) forall the m edges in ˜ O ( mθ − (cid:15) − ) time. The key ingredients of algorithm EdgeCentComp1 areSchur complements and Cholesky factorizations, which have been used in various applications,such as solving linear systems in Laplacians [KLP +
16, KS16] and counting and sampling spanningtrees [DKP +
17, DPPR17]. And the key technique for algorithm
EdgeCentComp2 is the combina-tion of sketching with the Sherman-Morrison formula [SM50] from efficient maintenances of matrixinverses for optimization [LSW15].The performance of the algorithm
EdgeCentComp1 is characterized in the following theorem.
Theorem 1.1.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positive edgeweights w : E → R + , and scalars < θ ≤ / , < (cid:15) ≤ / , the algorithm EdgeCentComp1 ( G =( V, E ) , w, θ, (cid:15) ) returns a set of pairs ˆ C = { ( e, ˆ c e ) | e ∈ E } . With high probability, the following tatement holds: For ∀ e ∈ E , C θ ( e ) ≈ (cid:15) ˆ c e , (1) where C θ ( e ) = X u,v ∈ V R G \ θ e eff ( u, v ) is the sum of effective resistances R G \ θ e eff ( u, v ) over all vertex pairs u and v in graph G \ θ e . The totalrunning time of this algorithm is bounded by O ( m(cid:15) − log m log n polyloglog( n )) . The proof of this theorem appears in Section 4.In Theorem 1.1, θ is arbitrary and can even depend on n , e.g. 1 /n . When θ is constant, wecan give a simpler algorithm EdgeCentComp2 that approximates C ∆ θ -Kirchhoff edge centralityfor all m edges in ˜ O ( mθ − (cid:15) − ) time. The idea is to use the Sherman-Morrison formula, whichgives a fractional expression of the difference between L † and ( L \ θ e ) † , where we can approximatethe numerator by the Johnson-Lindenstrauss lemma, and the denominator by estimating effectiveresistances. The technique is similar to the approach in [LSW15]. Theorem 1.2.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positive edgeweights w : E → R + , and scalars < θ ≤ / , < (cid:15) ≤ / , the algorithm EdgeCentComp2 ( G =( V, E ) , w, θ, (cid:15) ) returns a set of pairs ˆ C = { ( e, ˆ c ∆ e ) | e ∈ E } . With high probability, the followingstatement holds: For ∀ e ∈ E , C ∆ θ ( e ) ≈ (cid:15) ˆ c ∆ e , (2) where C ∆ θ ( e ) = K ( G \ θ e ) − K ( G ) . The total running time of this algorithm is bounded by O ( mθ − (cid:15) − log . n log(1 /(cid:15) ) polyloglog( n )) . The proof of this theorem appears in Section 5. Its advantage is that for moderate values of θ ,it can obtain a more accurate estimate of C ∆ θ ( e ) even if K ( G ) is large. However, when θ is small,the removal of high effective resistance edges can cause a large error in this routine, and we are notguaranteed to even get a good estimate of C θ ( e ) by adding this result to an estimate of K ( G ). As aresult we believe both of our algorithms for estimating Kirchhoff edge centrality are of interest, andcomplement each other.Based on the same idea of the definition for C ∆ θ ( e ), we can define a centrality measure for anyvertex v , which is the difference of Kirchhoff indices between the new graph G \ θ E v and originalgraph G , where G \ θ E v is obtained from G by multiplying the weights of all edges incident with v by θ . We write C ∆ θ ( v ) to denote the θ -Kirchhoff vertex centrality of v . In this situation, the matrixperturbation caused by removing the neighborhood of v is no longer rank 1, and we need to leveragethe approximate Schur complement routines from Section 4.2.2 to compute these intermediatematrices. Specifically, for constant θ , we give an algorithm VertexCentComp that approximates θ -Kirchhoff vertex centrality for all n vertices in ˜ O ( mθ − . (cid:15) − ) time. Theorem 1.3.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positive edgeweights w : E → R + , and scalars < θ ≤ / , < (cid:15) ≤ / , the algorithm VertexCentComp ( G =( V, E ) , w, θ, (cid:15) ) returns a set of pairs ˆ C = { ( v, ˆ c ∆ v ) | v ∈ V } . With high probability, the followingstatement holds: For ∀ v ∈ V , C ∆ θ ( v ) ≈ (cid:15) ˆ c ∆ v , (3)5 ve e Figure 1: Betweenness cannot distinguish between edges e and e . e e K K Figure 2: Spanning edge centrality cannot distinguish between edges e and e . where C ∆ θ ( v ) = K ( G \ θ E v ) − K ( G ) and E v = { ( u, v ) | u ∼ v } is the set of edges incident with v . The total running time of this algorithmis bounded by O ( m ( θ − (cid:15) − log n + θ − . (cid:15) − log n log(1 /(cid:15) )) polyloglog( n )) . The proof of this theorem appears in Section 6.
In addition to the low computational complexity, the θ -Kirchhoff centrality is more discriminatingthan other edge centrality measures, such as edge betweenness centrality and spanning edge centrality.For example, in the graph in Figure 1, the importance of edge e and edge e are different, whichcan be seen by intuition. In fact, we can also understand this difference from the influences whenthe two edges are deleted. If e is removed, the length of shortest path between vertices u and v increases by 6, while the removal of e will increase the length of shortest path between any pair ofvertices by at most 1. However, the betweenness centrality for e and e are the same, being equalto 18, implying that betweenness centrality cannot differentiate e between e . However, these twoedges can be discriminated by the θ -Kirchhoff edge centrality. Exact computation shows that the0 . e and e is C . ( e ) = 132 .
65 and C . ( e ) = 112 .
34, respectively.Thus, e is relatively more important than e , which agrees with our human intuition.We continue to show that the θ -Kirchhoff centrality is also more discriminating than the spanningedge centrality. By intuition, the importance for the two edges e and e of the graph illustrated inFigure 2 are distinct. Unfortunately, spanning edge centrality cannot distinguish these two edges,since their spanning edge centrality is identical, both equalling 1. In contrast, the 0 . e and e is C . ( e ) = 467 .
33 and C . ( e ) = 197 .
33, respectively. Thisimplies that e plays a relatively more significant role than e , which is consistent with our intuition.To further show the capability of our θ -Kirchhoff edge centrality to discriminate between differentedges, we experimentally compare our measure C ∆ θ with other metrics, including edge centrality,spanning edge centrality, and current-flow edge centrality. For each measure, we numerically evaluatethe importance of each edge for some classic real-world networks in Table 1. The data sets are from +
03] 62 159Celegansneural [WS98] 297 2148
Karate Lesmis Adjnoun Dolphins Celegans-neural R e l a t i v e s t a n d a r d d e v i a t i o n BetweennessCurrent-flow centralitySpanning edge centrality0.1-Kirchhoff edge centrality
Figure 3: Relative standard deviation for different edge centrality measures.published data mining related papers [Zac77, Knu93, New06, LSB +
03, WS98]. Based on whichwe then compute the relative standard deviation for each centrality measure (as the authors didin [BWLM16]), where the relative standard deviation is defined as the standard deviation dividedby the average. Figure 3 shows the relative standard deviation for all the centrality measures. It isalways significantly higher for θ -Kirchhoff edge centrality than it is for other measures, meaningthat our measure has a better capability to distinguish between different edges. The remaining part of the paper is organized as follows. In Section 2, we present the backgroundand formulate the problem of computing θ -Kirchhoff centrality. In Section 3, we introduce Schurcomplements and partial Cholesky factorizations, and a lemma with regard to the performanceof the approximate partial Cholesky factorization algorithm in [DKP + EdgeCentComp1 that approximates C θ ( e ). In Section 5, we introduce ouralgorithm EdgeCentComp2 that approximates C ∆ θ ( e ). In Section 6, we introduce our algorithm7 ertexCentComp that approximates C ∆ θ ( v ). In Section 7, we give our conclusion and discusssome directions for future works. We use the notion of (cid:15) -approximation in [PS14].Let a, b ≥ a is an (cid:15) -approximation of b ifexp( − (cid:15) ) a ≤ b ≤ exp( (cid:15) ) a. (4)We write a ≈ (cid:15) b to denote Eq. (4).For two matrices A and B , we write A (cid:22) B to indicate that B − A is positive semidefinite. Wesay A is an (cid:15) -spectral approximation of B ifexp( − (cid:15) ) A (cid:22) B (cid:22) exp( (cid:15) ) A . (5)We write A ≈ (cid:15) B to denote Eq. (5).Note that these two relations are symmetric. Namely, a ≈ (cid:15) b implies b ≈ (cid:15) a and A ≈ (cid:15) B implies B ≈ (cid:15) A .The following facts are basic properties of (cid:15) -approximation: Fact 2.1.
For nonnegative scalars a, b, c, d ≥ , positive semidefinite matrices A , B , C , D ,1. if a ≈ (cid:15) b , then a + c ≈ (cid:15) b + c ;2. if a ≈ (cid:15) b and c ≈ (cid:15) d , then a + c ≈ (cid:15) b + d ;3. if a ≈ (cid:15) b and b ≈ (cid:15) c , then a ≈ (cid:15) + (cid:15) c ;4. if a and b are positive such that a ≈ (cid:15) b , then /a ≈ (cid:15) /b ;5. if a ≈ (cid:15) b , then ac ≈ (cid:15) bc ;6. if A ≈ (cid:15) B , then A + C ≈ (cid:15) B + C ;7. if A ≈ (cid:15) B and C ≈ (cid:15) D , then A + C ≈ (cid:15) B + D ;8. if A ≈ (cid:15) B and B ≈ (cid:15) C , then A ≈ (cid:15) + (cid:15) C ;9. if A and B are positive definite matrices such that A ≈ (cid:15) B , then A − ≈ (cid:15) B − ;10. if A ≈ (cid:15) B and V is a matrix, then V > AV ≈ (cid:15) V > BV . We consider a connected undirected graph G = ( V, E ) with n vertices, m edges, and positive edgeweights w : E → R + . For a pair of vertices u, v ∈ E , we write u ∼ v to denote ( u, v ) ∈ E . TheLaplacian matrix of G is an n × n matrix L with the entry on its u th row and v th column being L ( u, v ) = − w ( u, v ) if u ∼ v, deg( u ) if u = v, , u ) = P u ∼ v w ( u, v ). If A and B are two sets of vertices in G , we write L AB to denote thesubmatrix of L with rows corresponding to A and columns corresponding to B .Let e i denote the i th standard basis vector, and b u,v = e u − e v . We fix an arbitrary orientationof the edges in G . For each edge e ∈ E , we define b e = b u,v , where u and v are head and tail of e ,respectively. It is easy to show that L = P e ∈ E w ( e ) b e b > e . We refer to w ( e ) b e b > e as the Laplacin of e . It is immediate that L is positive semidefinite since x > Lx = x > X e ∈ E w ( e ) b e b > e ! x = X e ∈ E w ( e ) x > b e b > e x = X e ∈ E w ( e ) (cid:16) x > b e (cid:17) ≥ x ∈ R n . Since L is positive semidefinite, we can diagonalize it and write L = n − X i =1 λ i v i v > i , where λ , . . . , λ n − are the nonzero eigenvalues of L and v , . . . , v n − are the corresponding or-thonormal eigenvectors. The pseudoinverse of L is defined as L † = n − X i =1 λ i v i v > i . It is not hard to show that if L and H are Laplacian matrices of connected graphs and L ≈ (cid:15) H ,then L † ≈ (cid:15) H † .We then give the definitions of effective resistance and Kirchhoff index: Definition 2.2 (Effective Resistance) . For a connected undirected graph G = ( V, E ), the effectiveresistance between u and v is defined as R G eff ( u, v ) = b > u,v L † b u,v . Definition 2.3 (Kirchhoff Index) . The Kirchhoff index K ( G ) of a graph G = ( V, E ) is defined asthe sum of effective resistances over all vertex pairs. Namely, K ( G ) = X u,v ∈ V R G eff ( u, v ) . For a graph, Kirchhoff index is a measure of its overall connectedness. A graph with smallerKirchhoff index is better connected on an average. It is known that the Kirchhoff index of a graphequals n times the sum of reciprocals of nonzero eigenvalues of L [ESVM + n times the trace of L † . We give this relation in the following Fact: Fact 2.4.
Let λ , . . . , λ n − be the nonzero eigenvalues of L . The Kirchhoff index of graph G satisfies K ( G ) = n n − X i =1 λ i = n Tr (cid:16) L † (cid:17) . By Rayleigh’s Monotonicity Law [ESVM + Fact 2.5.
The Kirchhoff index of a graph does not decrease when edges are deleted or edge weightsare decreased. .4 Kirchhoff Edge Centrality With Fact 2.5, it is reasonable to measure the importance of an edge e in graph G by the Kirchhoffindex of the new graph in which e is deactivated. We formalize the notion of edge deactivation bydefining θ -deletion of an edge. Definition 2.6 ( θ -Deletion) . Let 0 < θ ≤ / e ∈ E , the θ -deletion of e isto decrease its weight from w ( e ) to θw ( e ). If we use G \ θ e to denote the graph obtained from G by θ -deleting edge e , and L \ θ e to denote the Laplacin matrix corresponding to G \ θ e , we have L \ θ e = L − (1 − θ ) w ( e ) b e b > e . We can also define θ -Deletion of an edge set B ⊂ E , which is to decrease the weight of each edge e ∈ B from w ( e ) to θw ( e ). Similar notations are G \ θ B and L \ θ B = L − (1 − θ ) X e ∈ B w ( e ) b e b > e . We then give the definitions of θ -Kirchhoff edge centrality C θ ( e ) and C ∆ θ ( e ). Definition 2.7 ( θ -Kirchhoff Edge Centrality C θ ) . Let 0 < θ ≤ / e ∈ E ,its θ -Kirchhoff edge centrality C θ ( e ) is defined as the Kirchhoff index of the graph obtained from G by θ -deleting e . Namely, C θ ( e ) = K ( G \ θ e ) . Definition 2.8 ( θ -Kirchhoff Edge Centrality C ∆ θ ) . Let 0 < θ ≤ / e ∈ E ,its θ -Kirchhoff edge centrality C ∆ θ ( e ) is defined as the increase of the Kirchhoff index of the graphupon edge e ’s θ -deletion. Namely, C ∆ θ ( e ) = K ( G \ θ e ) − K ( G ) . Clearly, these two definitions of Kirchhoff edge centrality lead to the same ranking of edges.Following the above two definitions of Kirchhoff edge centrality, we can also define θ -centralityof a vertex. Definition 2.9 ( θ -Kirchhoff Vertex Centrality) . Let 0 < θ ≤ / v ∈ V ,its θ -Kirchhoff vertex centrality C ∆ θ ( v ) is defined as the increase of the Kirchhoff index of the graphupon θ -deletion of its incident edges. Namely, C ∆ θ ( v ) = K ( G \ θ E v ) − K ( G ) , where E v = { ( u, v ) | u ∼ v } is the set of edges incident with v .We now formulate the core problems of approximating θ -Kirchhoff centrality: Problem 1.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, and positiveedge weights w : E → R + , and scalars 0 < θ ≤ /
2, 0 < (cid:15) ≤ /
2, for each e ∈ E , find an (cid:15) -approximation of its θ -Kirchhoff edge centrality C θ ( e ). Problem 2.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, and positiveedge weights w : E → R + , and scalars 0 < θ ≤ /
2, 0 < (cid:15) ≤ /
2, for each e ∈ E , find an (cid:15) -approximation of its θ -Kirchhoff edge centrality C ∆ θ ( e ). Problem 3.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, and positiveedge weights w : E → R + , and scalars 0 < θ ≤ /
2, 0 < (cid:15) ≤ /
2, for each v ∈ V , find an (cid:15) -approximation of its θ -Kirchhoff vertex centrality C ∆ θ ( v ).10 Schur Complements and Partial Cholesky Factorizations
In this section, we introduce Schur complements and partial Cholesky factorizations, which are keytechniques in our algorithm.
We first give the definitions of Schur complements and partial Cholesky factorizations accordingto [KS16, DKP + Definition 3.1 (Schur Complement) . Suppose L is the Laplacian of an undirected positive-weightedconnected graph, and L (: , i ) is the i th column of L . For a vertex v , S (1) def = L − L ( v , v ) L (: , v ) L (: , v ) > is called the Schur complement of L with respect to vertex v . The operation of subtracting L ( v ,v ) L (: , v ) L (: , v ) > from L is called the elimination of vertex v . Suppose we perform asequence of eliminations, where in the i th step, we select a vertex v i ∈ V \ { v , . . . , v i − } andeliminate the vertex v i . We define α i = S ( i − ( v i , v i ) , c i = 1 α i S ( i − (: , v i ) , S ( i ) = S ( i − − α i c i c > i . Then S ( i ) is called the Schur complement with respect to vertices { v , . . . , v i } . Let C = { v , . . . , v i } ,then we also write Sc ( L , C ) = S ( i ) CC to denote the Schur complement of L onto C . Definition 3.2 (Partial Cholesky factorization) . Suppose we eliminate a sequence of vertices v , . . . , v i . Let L be the n × i matrix with c j as its j th column, and D be the i × i diagonal matrix D ( j, j ) = α j , then L = S ( i ) + i X j =1 α j c j c > j = S ( i ) + LDL > . Let us write F = { v , . . . , v i } , and C = V \ F . Let S be the submatrix of S ( i ) with rows andcolumns corresponding to vertices in C , i.e., S = S ( i ) CC . Since S ( i ) CC contains all nonzero entries of S ( i ) , we can write L = L F F L CF ! and L = L F F L CF ! D L F F L CF ! > + F F F C CF S ! = L F F L CF I CC ! D
00 S ! L F F L CF I CC ! > . (6)Here L F F L CF I CC ! is a lower triangular matrix up to row exchanges, and D is diagonal. Eq. (6) isknown as partial Cholesky factorization .It is known that Schur complements of a Laplacian are also Laplacians: Fact 3.3 (Fact 5.1 of [DKP + . The Schur complement of a Laplacian w.r.t. vertices v , . . . , v i is a Laplacian. .2 Commutativity With Edge Deletions According to [KS16], we can write the Schur complement w.r.t. a vertex v as S (1) = X e ∈ E : e v w ( e ) b e b > e + X u ∼ v X v ∼ v w ( u, v ) w ( v, v )deg( v ) b u,v b > u,v , (7)where the first term on rhs is the Laplacian corresponding to the edges not incident with v , and thesecond term on rhs is a Laplacian whose edges are supported on V \ { v } . Thus, S (1) can be seenas a multigraph obtained by adding edges to G [ V \ { v } ], the induced graph of G on V \ { v } . Byinduction, for all i , S ( i ) can be seen as a multigraph obtained by adding edges to G [ V \ { v , . . . , v i } ],the induced graph of G on V \ { v , . . . , v i } . Also, by Eq. (7), edges added to G [ V \ { v } ] to obtain S (1) are fully determined by edges incident with v in the original graph G . By induction, for all i , edges added to G [ V \ { v , . . . , v i } ] to obtain S ( i ) are fully determined by edges incident with { v , . . . , v i } in the original graph G . Thus, deletions (or θ -deletions) performed to edges with bothendpoints in V \ { v , . . . , v i } commute with taking partial Cholesky factorization. Therefore, wehave the following lemma: Lemma 3.4.
Given a connected undirected graph G = ( V, E ) , with positive edge weights w : E → R + ,and associated Laplacian L , a set of vertices F = { v , . . . , v i } ∈ V . Let C = V \ F , and the partialCholesky factorization of L be L = L F F L CF ! D L F F L CF ! > + F F F C CF S ! = L F F L CF I CC ! D
00 S ! L F F L CF I CC ! > . For any edge e whose endpoints are both in C , any ≤ θ < , L \ θ e = L F F L CF ! D L F F L CF ! > + F F F C CF S \ θ e ! = L F F L CF I CC ! D
00 S \ θ e ! L F F L CF I CC ! > . In [DKP + Lemma 3.5 (Lemma 5.7 of [DKP + . There is an algorithm
ApxPartialCholesky ( L , C, (cid:15) ) that when given a connected undirected graph G = ( V, E ) , withpositive edge weights w : E → R + , and associated Laplacian L , a set of vertices C ⊂ V , and a scalar < (cid:15) ≤ / , returns a decomposition ( e L , f D , e S ) . With high probability, the following statement hold: L ≈ (cid:15) ˜L , (8) where F = V \ C and ˜L = e L F F e L CF ! f D e L F F e L CF ! > + F F F C CF e S ! = e L F F e L CF I CC ! f D e S ! e L F F e L CF I CC ! > . (9) Here e S is a Laplacian matrix whose edges are supported on C such that e S ≈ (cid:15) Sc ( L , C ) . Let k = | C | = n − | F | . The total number of non-zero entries in e S is O ( k(cid:15) − log n ) . e L F F e L CF I CC ! is a ower triangular matrix up to row exchanges. The total number of non-zero entries in e L F F e L CF I CC ! is O ( m + n(cid:15) − log n ) . f D is a diagonal matrix.For any vector b ∈ R n , one can evaluate e L F F e L CF I CC ! − b in O ( m + n(cid:15) − log n ) time. Forany vector c ∈ R | F | , one can evaluate (cid:16) f D (cid:17) − c in O ( | F | ) time.The total running time is bounded by O (( m log n + n(cid:15) − log n ) polyloglog( n )) . Comparing to [DKP + δ and just claims high probability. The running time ofthe algorithm ApxPartialCholesky ( L , C, (cid:15) ) in [DKP +
17] is O (( m log n log ( n/δ ) + n(cid:15) − log n log ( n/δ )) polyloglog( n )). To gain high probability, we justset the failure probability δ to 1 /n c for an arbitrary constant c >
0. Then we have the runningtime bounded by O (( m log n + n(cid:15) − log n ) polyloglog( n )).2. emphasizes that inverses of matrices e L F F e L CF I CC ! and f D can both be applied quickly, asthey can be treated as lower triangular matrix and diagonal matrix, respectively.The following lemma shows that edge additions performed within C commute with takingapproximate partial Choleksy factorization: Lemma 3.6.
Given a connected undirected multi-graph G = ( V, E ) , with positive edge weights w : E → R + , and associated Laplacian L , a set of vertices C ⊂ V , and an approximate partialfactorization of L : L ≈ (cid:15) e L F F e L CF I CC ! f D e S ! e L F F e L CF I CC ! > , (10) where F = V \ C . For any edge e (not necessarily in E ) with both endpoints in C and a positivescalar w e > , L + w e b e b > e ≈ (cid:15) e L F F e L CF I CC ! f D e S + (cid:16) w e b e b > e (cid:17) CC ! e L F F e L CF I CC ! > . (11) Proof.
As multiplicative approximations are preserved under additions, by adding w e b e b > e to bothsides of Eq. (10) we have L + w e b e b > e ≈ (cid:15) e L F F e L CF I CC ! f D e S ! e L F F e L CF I CC ! > + w e b e b > e = e L F F e L CF ! f D e L F F e L CF ! > + F F F C CF e S + (cid:16) w e b e b > e (cid:17) CC ! = e L F F e L CF I CC ! f D e S + (cid:16) w e b e b > e (cid:17) CC ! e L F F e L CF I CC ! > . Algorithm for Approximating θ -Kirchhoff Edge Centrality C θ ( e ) L † By Fact 2.4, the Kirchhoff Index of a graph equals n times the trace the Laplacian’s pseudoinverse.Although the explicit pseudoinverse of L is hard to compute, by taking approximate Choleskyfactorizations [KS16, DKP + z > L † z for a z ∈ R n quickly. Thus, we canuse Monte-Carlo methods to estimate trace of L † .The standard Monte-Carlo method for estimating the trace of an implicit matrix A is due toHutchinson [Hut89]. The idea is to estimate the trace of A by M P Mi =1 z > i Az i , where the z i ’s arerandom ± E h z > i Az i i = Tr ( A ), by thelaw of large numbers, M P Mi =1 z > i Az i should be close to Tr ( A ) when M is large. [AT11] gives arigorous bound on the number of Monte-Carlo samples required to achieve a maximum error (cid:15) withprobability at least 1 − δ . Lemma 4.1 (Theorem 7.1 of [AT11], paraphrased) . Let A be a positive semidefinite matrix withrank rank( A ) . Let z , . . . , z M be independent random ± vectors. Let (cid:15), δ be scalars such that < (cid:15) ≤ / and < δ < . For any M ≥ (cid:15) − ln(2rank( A ) /δ ) , the following statement holds withprobability at least − δ : M M X i =1 z > i Az i ≈ (cid:15) Tr ( A ) . Remark 4.2.
We remark that the Hutchinson’s method can be seen as Johnson-LindenstraussLemma [JL84] in some sense. The reason is that since A is positive semidefinite, one can write itstrace as Tr ( A ) = Tr (cid:16) A / A / (cid:17) = (cid:13)(cid:13)(cid:13) A / (cid:13)(cid:13)(cid:13) F , where (cid:13)(cid:13)(cid:13) A / (cid:13)(cid:13)(cid:13) F can be seen as a sum of the squared lengths of the rows of A / . By the discreteversion of Johnson-Lindenstrauss Lemma from [Ach01], we can use a n × k random ± Q ,where k = O ( (cid:15) − log n ), to reduce the dimensions: (cid:13)(cid:13)(cid:13) A / (cid:13)(cid:13)(cid:13) F ≈ (cid:15) k (cid:13)(cid:13)(cid:13) A / Q (cid:13)(cid:13)(cid:13) F . This in turn implies Tr ( A ) ≈ (cid:15) k k X j =1 q > j Aq j , (12)where q j is the j th column of Q . The rhs of (12) can be seen as Hutchinson’s method. Indeed, [AT11]used the discrete Johnson-Lindenstrauss Lemma from [Ach01] to prove their bound.Since L † is positive semidefinite and rank( L † ) = n −
1, by letting δ = 1 /n , we have the followingbound on the number of Monte-Carlo samples required to achieve an (cid:15) -approximation of Tr (cid:16) L † (cid:17) with high probability: 14 emma 4.3. Let L be a Laplacian matrix. Let z , . . . , z M be independent random ± vectors. Let (cid:15) be a scalar such that < (cid:15) ≤ / . For any M ≥ (cid:15) − ln(2 n ) , the following statement holds withprobability at least − /n : M M X i =1 z > i L † z i ≈ (cid:15) Tr (cid:16) L † (cid:17) . A direct conclusion of Lemma 4.3 is that for an edge e ∈ E , its θ -Kirchhoff edge centralitysatisfies C θ ( e ) = K ( G \ θ e ) = n Tr (cid:16) ( L \ θ e ) † (cid:17) ≈ (cid:15) nM M X i =1 z > i ( L \ θ e ) † z i . Thereby, the task of approximating the θ -Kirchhoff edge centrality for all e ∈ E can be dividedinto O ( (cid:15) − log n ) independent tasks, each of which is to compute quadratic forms z > ( L \ θ e ) † z for afixed z ∈ R n for all e ∈ E . We formulate these tasks in the following problem: Problem 4.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positiveedge weights w : E → R + , and associated Laplacian L , a set of edges E Q ⊂ E such that everyvertex in V is incident to some edge e ∈ E Q , a scalar 0 < θ ≤ /
2, and a vector z ∈ R n , find(approximately) z > L † z for all e ∈ E Q . L † Upon Edge Deactivation
The idea of solving Problem 4 is to use recursions based on partial Cholesky factorizations. Wesummarize the key steps in the following enumeration:1. If L only have O (1) vertices, invert L \ θ e to compute z > ( L \ θ e ) † z for all e ∈ E Q and return.2. Divide edges in E Q into E (1) , E (2) with equal sizes.3. Let C denote endpoints of edges in E (1) and F = V \ C .4. By taking (approximate) partial Cholesky factorization of L , find a vector y = y F y C ! , adiagonal matrix D | F |×| F | , and a Laplacian matrix S whose edges are supported on C , such thatfor each edge e ∈ E (1) , z > ( L \ θ e ) † z can be evaluated by computing y > F D − y F + y > C ( S \ θ e ) † y C .5. Compute y > F D − y F by inverting D and y > C ( S \ θ e ) † y C for all e ∈ E (1) by recursion, then use y > F D − y F + y > C ( S \ θ e ) † y C to evaluate z > ( L \ θ e ) † z for all e ∈ E (1) .6. Repeat steps 3 - 5 to E (2) to compute z > ( L \ θ e ) † z for all e ∈ E (2) .The reason that in Step 5 we can compute y > C ( S \ θ e ) † y C for each e ∈ E (1) by recursion is that S is a Laplacian matrix whose edges are supported on C , and hence to compute y > C ( S \ θ e ) † y C forall e ∈ E (1) is just a smaller-sized version of Problem 4 in which L = S and E Q = E (1) .In the rest of this subsection we give first an algorithm that solves Problem 4 exactly and thenan algorithm that solves Problem 4 approximately.15 .2.1 Computing Exact Quadratic Forms of L † Upon Edge Deactivation
We first give an algorithm
ExactQuad ( L , E Q , w, z , θ ) that computes the exact value of z > ( L \ θ e ) † z for a fixed z ∈ R n for all e ∈ E Q (Here w is the edge weight function).In this algorithm, we find y , D , and S in step 4 by eliminating vertices in F and obtain anexact partial Cholesky factorization of L . The following Lemma shows how to find them when anexact partial Cholesky factorization of L is given: Lemma 4.4.
For a graph G = ( V, E ) with associate Laplacian L and a set of vertices C ⊂ V . Let F = V \ C , and the partial Choleksy factorization of L be L = L F F L CF I CC ! D
00 S ! L F F L CF I CC ! > . Let y = y F y C ! = L F F L CF I CC ! − z , then for each edge e ∈ E with both endpoints in C the followingstatement holds: z > ( L \ θ e ) † z = y > F (cid:16) D − (cid:17) y F + y > C ( S \ θ e ) † y C . Proof.
By Lemma 3.4, θ -deletions performed to edges with both endpoints in C commute withtaking partial Cholesky factorization. Thus, for each e ∈ E with both endpoints in C , we have L \ θ e = L F F L CF I CC ! D
00 S \ θ e ! L F F L CF I CC ! > . (13)Inverting both sides of Eq. (13) leads to( L \ θ e ) † = L F F L CF I CC ! −> D − ( S \ θ e ) † ! L F F L CF I CC ! − . Substituting y = y F y C ! = L F F L CF I CC ! − z , we obtain z > ( L \ θ e ) † z = (cid:16) y F y C (cid:17) D − ( S \ θ e ) † ! y F y C ! = y > F (cid:16) D − (cid:17) y F + y > C ( S \ θ e ) † y C . This completes the proof.We give the pseudocode for
ExactQuad in Algorithm 1. Its performance is characterized inLemma 4.5.
Lemma 4.5.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positiveedge weights w : E → R + , and associated Laplacian L , a set of edges E Q ⊂ E such that every vertexin V is incident with some edge e ∈ E Q , a vector z ∈ R n , and a scalar < θ ≤ / , the algorithm ExactQuad ( L , E Q , w, z , θ ) returns a set of pairs N = { ( e, n e ) | e ∈ E Q } , where n e = z > ( L \ θ e ) † z . The total running time of this algorithm is bounded by O ( n ω − m ) . lgorithm 1: ExactQuad ( L , E Q , w, z , θ ) Input : L : A graph Laplacian. E Q : A set of edges supported on vertices in L . w : An edge weight function. z : A vector whose dimension matches the number of vertices in L . θ : The weight of edge e is temporarily changed from w ( e ) to θw ( e ) when it isdeactivated. Output : N = { ( e, n e ) | e ∈ E Q } : n e = z > ( L \ θ e ) † z . Let V denote the vertex set of L . if | V | = 2 then For every edge e ∈ E Q , compute exact n e = z > (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † z , then combinethe results and return N = { ( e, n e ) | e ∈ E Q } . Partition E Q into E (1) , E (2) with (cid:12)(cid:12)(cid:12) E (1) (cid:12)(cid:12)(cid:12) = (cid:22) | E Q | (cid:23) and (cid:12)(cid:12)(cid:12) E (2) (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) − (cid:22) | E Q | (cid:23) . for i = 1 to do Let C denote endpoints of edges in E ( i ) and F = V \ C . Eliminate all vertices in F to get L = L F F L CF I CC ! D
00 S ! L F F L CF I CC ! > . Compute y = y F y C ! = L F F L CF I CC ! − z . Compute f = y > F ( D ) − y F . Call
ExactQuad ( S , E ( i ) , y C , θ ) to compute n ( i ) e = y > C ( S \ θ e ) † y C for all e ∈ E ( i ) andstore ( e, f + n ( i ) e ) in N ( i ) . return N = N (1) ∪ N (2) Proof of Lemma 4.5.
As correctness is clear by Lemma 4.4, we only need to prove the bound ofrunning time. Let T ( m ) denote the running time of ExactQuad ( L , E Q , w, z , θ ), where m = (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) .Let n denote the number of vertices in the original graph, i.e., the graph corresponding to L in the earliest call to ExactQuad . Let n cur denote the number of vertices in L in the currentcall. If n cur = 2, the algorithm goes to Line 3, and hence we have T ( m ) = O (1). Otherwise,the algorithm goes to Lines 4 - 11, among which the most time-consuming work is eliminating F ,inverting L F F L CF I CC ! and D , and recursively calling ExactQuad . The first two both run in O ( n ω cur ) time, and the third runs in 2 T ( m/
2) time. When m > n , we can bound the number ofvertices in the current call by n cur = O ( n ); otherwise when m ≤ n , we can bound the number ofvertices in the current call by n cur = O ( m ). Thus, We have T ( m ) = ( T ( m/
2) + O ( n ω ) , m > n, T ( m/
2) + O ( m ω ) , m ≤ n. (14)Eq. (14) leads to T ( m ) = O ( n ω − m ). 17 .2.2 Approximating Quadratic Forms of L † Upon Edge Deactivation
Clearly, if we only want to approximately compute the quadratic forms, we can use the approximatepartial Choleksy algorithm in Lemma 3.5 to speed up. Thereby, we give an approximation algorithm
QuadEst ( L , E, w, z , θ, (cid:15) ), that computes an (cid:15) -approximation of z > ( L \ θ e ) † z for a fixed z ∈ R n forall e ∈ E Q . We also make a few modifications to maintain the error and further speed up. We listthe modifications in QuadEst below:1. In step 4, instead of computing the exact partial Cholesky factorization of L , we use thealgorithm ApxPartialCholesky in Lemma 3.5 to obtain an approximate partial Choleskyfactorization of L . However, if we pass the whole L to ApxPartialCholesky , it may changethe edges in E (1) , to which we need to perform θ -deletions when deactivating them. Thus,instead, we first delete all edges in E (
1) and pass the resulting L to ApxPartialCholesky ,and then add those edges back to the approximate Schur complement e S returned by it. Thismodification is feasible since adding edges with both endpoints in C commutes with takingapproximate partial Choleksy factorization (Lemma 3.6). This modification is addressed onLines 7 - 9 of Algorithm 2.2. By Lemma 3.5, matrices f D and e L returned by ApxPartialCholesky satisfy that f D isdiagonal, e L is sparse, and e L F F e L CF I CC ! is a lower triangular matrix up to row exchanges.Therefore, by applying inverses of diagonal matrix and lower triangular matrix quickly, wecan compute y = y F y C ! = e L F F e L CF I CC ! − z and y > F (cid:16) f D (cid:17) − y F in linear time of the numberof nonzero entries. This is addressed on Lines 10 - 11 of Algorithm 2.3. Since errors may accumulate among different levels of the recursion, we bound the error by (cid:15)/ log (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) when taking approximate partial Cholesky factorization (Line 8 of Algorithm 2), andbound the error by (cid:15) − (cid:15)/ log (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) when recursively calling QuadEst (Line 12 of Algorithm 2).Thereby, the errors add up to (cid:15) as required, and only an extra log m factor is added to therunning time (see Lemma 4.7 and its proof for details).According to the first modification, in this algorithm, we find y , D , and S in step 4 by takingapproximate partial Cholesky factorization. The following Lemma shows how to find them when anapproximate partial Cholesky factorization is given. Lemma 4.6.
For a graph G = ( V, E ) with associate Laplacian L and a set of vertices C ⊂ V . Let F = V \ C . Let E (1) ⊂ E be a set of edges with both endpoints in C , and H be the Laplacian matrixcorresponding to edges in E (1) , i.e., H = P e ∈ E (1) w ( e ) b e b > e . Clearly L − H is also a Laplacian. Letan approximate partial factorization of L − H be L − H ≈ (cid:15) e L F F e L CF I CC ! f D e S ! e L F F e L CF I CC ! > . (15) Let y = y F y C ! = e L F F e L CF I CC ! − z and e S = e S + H CC , then for each edge e ∈ E (1) the followingstatement holds: z > ( L \ θ e ) † z ≈ (cid:15) y > F (cid:16) f D − (cid:17) y F + y > C (cid:16) e S \ θ e (cid:17) † y C . roof. By Lemma 3.6, adding edges with both endpoints in C commutes with taking approximatepartial Choleksy factorization. Thus, for each edge e ∈ E (1) , by adding first edges in E (1) \ { e } andthen the deactivated edge e (i.e., edge e with weight θw ( e )), we have L − H + ( H \ θ e ) ≈ (cid:15) e L F F e L CF I CC ! f D e S + ( H CC \ θ e ) ! e L F F e L CF I CC ! > . Substituting (cid:16) e S \ θ e (cid:17) = e S + ( H CC \ θ e ) and L \ θ e = L − H + ( H \ θ e ) leads to L \ θ e ≈ (cid:15) e L F F e L CF I CC ! f D (cid:16) e S \ θ e (cid:17)! e L F F e L CF I CC ! > . (16)Note that e S is a Laplacian since it is a sum of two Laplacians. Inverting both sides of Eq. (16)leads to ( L \ θ e ) † ≈ (cid:15) e L F F e L CF I CC ! −> (cid:16) f D (cid:17) − (cid:16) e S \ θ e (cid:17) † e L F F e L CF I CC ! − . (17)Multiplying both sides of Eq. (17) by z > on the left and z on the right and substituting y = y F y C ! = e L F F e L CF I CC ! − z , gives z > ( L \ θ e ) † z ≈ (cid:15) (cid:16) y F y C (cid:17) f D − (cid:16) e S \ θ e (cid:17) † y F y C ! = y > F (cid:16) f D − (cid:17) y F + y > C (cid:16) e S \ θ e (cid:17) † y C . (18)This completes the proof.The pseudocode for QuadEst is given in Algorithm 2. Its performance is characterized inLemma 4.7.
Lemma 4.7.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positiveedge weights w : E → R + , and associated Laplacian L , a set of edges E Q ⊂ E such that every vertexin V is incident with some edge e ∈ E Q , a vector z ∈ R n , and scalars < θ ≤ / , < (cid:15) ≤ / ,the algorithm QuadEst ( L , E Q , w, z , θ, (cid:15) ) returns a set of pairs ˆ N = { ( e, ˆ n e ) | e ∈ E Q } . With highprobability, the following statement holds: For ∀ e ∈ E Q , n e ≈ (cid:15) ˆ n e , (19) where n e = z > ( L \ θ e ) † z . The total running time of this algorithm is bounded by O ( m(cid:15) − log m log n polyloglog( n )) . lgorithm 2: QuadEst ( L , E Q , w, z , θ, (cid:15) ) Input : L : A graph Laplacian. E Q : A set of edges supported on vertices in L . w : An edge weight function. z : A vector whose dimension matches the number of vertices in L . θ : An edge e ’s weight should be temporarily reduced to θw ( e ) whendeactivating it. (cid:15) : Error of the estimates. Output : ˆ N = { ( e, ˆ n e ) | e ∈ E Q } : ˆ n e is an estimate of n e = z > ( L \ θ e ) † z . Let V denote the vertex set of L . if | V | = 2 then For every edge e ∈ E Q , compute exact ˆ n e = n e = z > (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † z , thencombine the results and return ˆ N = { ( e, ˆ n e ) | e ∈ E Q } . Partition E Q into E (1) , E (2) with (cid:12)(cid:12)(cid:12) E (1) (cid:12)(cid:12)(cid:12) = (cid:22) | E Q | (cid:23) and (cid:12)(cid:12)(cid:12) E (2) (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) − (cid:22) | E Q | (cid:23) . for i = 1 to do Let C denote endpoints of edges in E ( i ) and F = V \ C . Let H denote the Laplacian matrix corresponding to edges in E ( i ) , i.e., H ← P e ∈ E ( i ) w ( e ) b e b > e . ( e L , f D , e S ) ← ApxPartialCholesky ( L − H , C, (cid:15)/ log (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) ) Add edges in E ( i ) back to e S and store the resulting Laplacian in e S , i.e., e S ← e S + H CC . Compute y = y F y C ! = e L F F e L CF I CC ! − z in linear time. Compute f = y > F (cid:16) f D (cid:17) − y F in linear time. Call
QuadEst ( e S , E ( i ) , y C , θ, (cid:15) − (cid:15)/ log (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) ) to get an estimate ˆ n ( i ) e of n ( i ) e = y > C (cid:16) e S \ θ e (cid:17) † y C for all e ∈ E ( i ) and store ( e, f + ˆ n ( i ) e ) in ˆ N ( i ) . return ˆ N = ˆ N (1) ∪ ˆ N (2) Proof of Lemma 4.7.
We first prove the error bound (i.e., Eq. (19)) by induction on the size of E Q .For (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) = 1, we have | V | = 2. Hence, the algorithm QuadEst will go into Line 3 and returnsan ˆ n e = n e , which implies that n e ≈ (cid:15) ˆ n e holds for any (cid:15) > ≤ (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) ≤ k , k ≥
1. We now prove that it holds for (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) = k + 1,too. Clearly, by symmetry, it suffices to show that n e ≈ (cid:15) ˆ n e holds for each e ∈ E (1) . By Lemma 3.5,matrices e L , f D , e S on Line 8 satisfy L − H ≈ (cid:15)/ log | E Q | e L F F e L CF I CC ! f D e S ! e L F F e L CF I CC ! > . (20)By Lemma 4.6, we have z > ( L \ θ e ) † z ≈ (cid:15)/ log | E Q | y > F (cid:16) f D − (cid:17) y F + y > C (cid:16) e S \ θ e (cid:17) † y C , (21)20here e S = e S + H CC (Line 9). Since (cid:12)(cid:12)(cid:12) E (1) (cid:12)(cid:12)(cid:12) = (cid:22) | E Q | (cid:23) ≤ k , by inductive assumption, each ˆ n (1) e returned by the recursive call QuadEst on Line 12 satisfies y > C (cid:16) ˆ S \ θ e (cid:17) † y C ≈ (cid:15) − (cid:15)/ log | E Q | ˆ n (1) e , whichwhen adding f = y > F (cid:16) D − (cid:17) y F (Line 11) to its both sides turns into y > F (cid:16) D − (cid:17) y F + y > C (cid:16) ˆ S \ θ e (cid:17) † y C ≈ (cid:15) − (cid:15)/ log | E Q | f + ˆ n (1) e . (22)Combining Eq. (21) and Eq. (22) and substituting n e = z > ( L \ θ e ) † z , we have n e ≈ (cid:15) f + ˆ n (1) e . Thus,Eq. (19) holds for (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) = k + 1, too. By induction, it holds for all (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) .We then prove the running time of the algorithm.Let T ( m, (cid:15) ) denote the running time of QuadEst ( L , E Q , w, z , θ, (cid:15) ), where m = (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) and (cid:15) isthe error of estimates. Let n denote the number of vertices in the original graph, i.e., the graphcorresponding to L in the earliest call to QuadEst . Let n cur and m cur denote the number of verticesand the number of edges in L in the current call, respectively. In each call other than the earliestcall, the Laplacian L equals e S on Line 9 of the parent call (i.e., the call that invoked current call),where we have the total number of edges in S being the number of edges in e S plus the number ofedges in H . Since H is the Laplician corresponding to edges in E ( i ) , which is precisely E Q in thecurrent call, we have the number of edges in H equaling m = (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) . By Lemma 3.5, the number ofedges in e S is O ( n cur (cid:15) − log m log n ), where there is an extra log m factor because the error is setto (cid:15)/ log (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) when calling ApxPartialCholesky on Line 8. Hence, the number of edges in L inthe current call is bounded by m cur = O ( m + n cur (cid:15) − log m log n ).If n cur = 2, the algorithm goes to Line 3, and hence we have T ( m, (cid:15) ) = O (1). Otherwise, thealgorithm goes to Lines 4 - 13, among which the most time-consuming work can be divided intothree parts:1. The first part is computing e L F F e L CF I CC ! − z and y > F (cid:16) f D (cid:17) − y F . Since e L F F e L CF I CC ! isa lower triangular matrix up to row exchanges and f D is diagonal, their inverse can beapplied in linear time of the number of nonzero entries. By Lemma 3.5, this part runs in O ( m cur + n cur ( (cid:15)/ log m ) − log n ) = O ( m + n cur (cid:15) − log m log n ) time.2. The second part is taking approximate partial Cholesky factorization, which by Lemma 3.5runs in O (( m cur log n + n cur ( (cid:15)/ log m ) − log n ) polyloglog( n )) = O ( m log n + n cur (cid:15) − log m log n polyloglog( n )) time.3. The third part is recursively calling QuadEst ( e S , E ( i ) , y C , θ, (cid:15) − (cid:15)/ log (cid:12)(cid:12)(cid:12) E Q (cid:12)(cid:12)(cid:12) ), which runs in2 T ( m/ , (cid:15) − (cid:15)/ log m ) time.The first two parts add up to a running time of O ( m log n + n cur (cid:15) − log m log n polyloglog( n )).When m > n , we can bound the number of vertices in the current call by n cur = O ( n ); otherwisewhen m ≤ n , we can bound the number of vertices in the current call by n cur = O ( m ). Thus, Wehave T ( m, (cid:15) ) = ( T ( m/ , (cid:15) − (cid:15)/ log m ) + O ( m log n + n(cid:15) − log m log n polyloglog( n )) , m > n T ( m/ , (cid:15) − (cid:15)/ log m ) + O ( m log n + m(cid:15) − log m log n polyloglog( n )) , m ≤ n . (23)Eq. (23) leads to T ( m, (cid:15) ) = O ( m(cid:15) − log m log n polyloglog( n )).21 .3 Approximating C θ ( e ) We are now ready to give the algorithm
EdgeCentComp1 ( G = ( V, E ) , w, θ, (cid:15) ), which computesan (cid:15) -approximation of the θ -Kirchhoff edge centrality C θ ( e ) for all e ∈ E . The pseudocode for EdgeCentComp1 is given in Algorithm 3. Its performance is characterized in Theorem 1.1.
Algorithm 3:
EdgeCentComp1 ( G = ( V, E ) , w, θ, (cid:15) ) Input : G = ( V, E ), w : A connected undirected graph with positive edgesweights w : E → R + . θ : An edge e ’s weight should be temporarily reduced to θw e whendeactivating it. (cid:15) : Error of the centrality estimate per edge. Output : ˆ C = { ( e, ˆ c e ) | e ∈ E } : ˆ c e is an estimate of C θ ( e ), the θ -Kirchhoffedge centrality of e . Let z , . . . , z M be independent random ± M = (cid:6) (cid:15) − ln(2 n ) (cid:7) . for i = 1 to M do Call
QuadEst ( L G , E, w, z i , θ, (cid:15) ) to get an estimate ˆ n e of n e = z > i ( L \ θ e ) † z i for each e ∈ E , and store each ˆ n e in ˆ n ( i ) e . For each e ∈ E compute ˆ c e = nM M P i =1 ˆ n ( i ) e and return ˆ C = { ( e, ˆ c e ) | e ∈ E } . Proof of Theorem 1.1.
The running time is the total cost of O ( (cid:15) − log n ) calls to QuadEst , each of which runs in O ( m(cid:15) − log m log n polyloglog( n )) time according to Lemma 4.7.Since M = (cid:6) (cid:15) − ln(2 n ) (cid:7) ≥ (cid:0) (cid:15) (cid:1) − ln(2 n ), by Lemma 4.3, for each e ∈ E , there is1 M M X i =1 z > i ( L \ θ e ) † z i ≈ (cid:15) Tr ( L \ θ e ) . Multiplying both sides by n and substituting C θ ( e ) = n Tr ( L \ θ e ), we have nM k X i =1 z > i ( L \ θ e ) † z i ≈ (cid:15) C θ ( e ) . (24)By Lemma 4.7, each ˆ n ( i ) e on Line 3 satisfies n ( i ) e ≈ (cid:15) ˆ n ( i ) e , (25)where n ( i ) e = z > i ( L \ θ e ) † z i . Summing Eq. (25) over i = 1 , . . . , M and multiplying both sides by nM lead to nM M X i =1 z > i ( L \ θ e ) † z i ≈ (cid:15) nM M X i =1 ˆ n ( i ) e . Combining with Eq. (24) and substituting ˆ c e = nM M P i =1 ˆ n ( i ) e (Line 4) lead to C θ ( e ) ≈ (cid:15) ˆ c e .22 Algorithm for Approximating θ -Kirchhoff Edge Centrality C ∆ θ ( e ) By Sherman-Morrison formula, for an edge e ∈ E and a scalar 0 < θ <
1, we have( L \ θ e ) † = (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † = L † + (1 − θ ) w ( e ) L † b e b > e L † − (1 − θ ) w ( e ) b > e L † b e . (26)Since the off-the-shelf Sherman-Morrison formula is for full rank matrices, we give the detailed proofof Equation (26) in Appendix A.Since by our definition C ∆ θ ( e ) = K ( G \ θ e ) − K ( G ) = n (cid:16) Tr (cid:16) ( L \ θ e ) † (cid:17) − Tr (cid:16) L † (cid:17)(cid:17) , it follows that C ∆ θ ( e ) = n (1 − θ ) w ( e )Tr (cid:16) L † b e b > e L † (cid:17) − (1 − θ ) w ( e ) b > e L † b e . (27)The numerator of (27) is the trace of an implicit matrix, and hence can be approximated byHutchinson’s [AT11, Hut89] Monte-Carlo method. To apply L † , we can utilize nearly-linear timesolvers for Laplacian systems [ST14, CKM + + Lemma 5.1 (Theorem 1.1 of [CKM + . There is an algorithm y = LaplSolve ( L G , z , δ ) which takes a Laplacian matrix L G of a graph G with n vertices and m edges, a vector z ∈ R n , and a scalar δ > , and returns a vector y ∈ R n such that with highprobability the following statement holds: (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) L ≤ δ (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L , where k x k L = √ x > Lx . The algorithm runs in expected time O ( m log . n log(1 /δ ) polyloglog( n )) . To track the error for the solver, we will need the following two lemmas, whose proofs aredeferred to Appendix C.1.
Lemma 5.2.
Let L be the Laplacian of a graph with all weights in the range [1 , U ] , and z beany vector such that k z k ≤ n . Suppose y is a vector such that (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) L ≤ δ (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L for some < δ < . For any edge e of the graph, we have (cid:12)(cid:12)(cid:12) y > b e b > e y − z > L † b e b > e L † z (cid:12)(cid:12)(cid:12) ≤ δn U . (28) Lemma 5.3.
Let L be the Laplacian of a graph with all weights in the range [1 , U ] . For any edge e of the graph, we have Tr (cid:16) L † b e b > e L † (cid:17) ≥ n U . The denominator of (27) is just 1 − (1 − θ ) w ( e ) R eff ( e ). Since w ( e ) R eff ( e ) is between 0 and 1and θ is positive, (1 − θ ) w ( e ) R eff ( e ) is strictly bounded away from 1. Thus, we can multiplicativelyapproximate the denominator by approximating R eff ( e ), for which we can use the random projectionin [SS11]. By Using the solvers from [CKM +
14] in the effective resistance estimation procedureof [SS11], we immediately have the following lemma:
Lemma 5.4.
There is an algorithm
EREst ( G, (cid:15) ) that when given a graph G = ( V, E ) , returns anestimate ˆ r e of R eff ( e ) for all e ∈ E in O ( m(cid:15) − log . n polyloglog( n )) time. With high probability, ˆ r e ≈ (cid:15) R eff ( e ) holds for all e ∈ E .
23s these estimates are approximate, we will need to bound their approximations when subtractedfrom 1. Here we use the fact that 0 < θ <
1, and that the weight times effective resistance of anedge, w ( e ) b Te L † b e is between 0 and 1. Since we will also need the matrix version of this type ofapproximation propagation when computing Kirchhoff vertex centralities in Section 6, we will statethe more general version here. Lemma 5.5. If A and B are matrices such that (cid:22) A (cid:22) I , and A ≈ (cid:15) B for some < (cid:15) ≤ / ,then for any < θ ≤ / such that (cid:15)/θ ≤ / , we have I − (1 − θ ) A ≈ (cid:15)/θ I − (1 − θ ) B . The proof is deferred to Appendix B.We then give an algorithm
EdgeCentComp2 to approximate the θ -edge Kirchhoff centrality C ∆ θ for all m ∈ E . The pseudocode for EdgeCentComp2 is given in Algorithm 4. The performanceof
EdgeCentComp2 is characterized in Theorem 1.2.
Algorithm 4:
EdgeCentComp2 ( G, θ, (cid:15) ) Input : G : A graph. θ : a scalar between 0 and 1 / (cid:15) : the error parameter. Output : ˆ C ∆ = { ( e, ˆ c ∆ e ) | e ∈ E } . Let z , . . . , z M be independent random ± M = (cid:6) (cid:15) − ln(2 n ) (cid:7) . for i = 1 to M do y i ← LaplSolve ( L G , z i , (cid:15)n − U − ) for each e ∈ E do Compute ˆ c ∆( i ) e = y > i b e b > e y i . ˆ r e ← EREst ( G, θ(cid:15)/ Compute ˆ c ∆ e = (1 − θ ) nM w ( e ) M P i ˆ c ∆( i ) e / (1 − (1 − θ ) w ( e )ˆ r e ) for each e . Proof of Theorem 1.2.
The running time is the total cost of O ( (cid:15) − log n ) calls to LaplSolve eachof which runs in O ( m log . n log(1 /(cid:15) ) polyloglog( n )) time, and a call to EREst which runs in O ( mθ − (cid:15) − log . n polyloglog( n )) time.Since M = (cid:6) (cid:15) − ln(2 n ) (cid:7) ≥
48 ( (cid:15)/ − ln(2 n ), by Lemma 4.3, we have1 M M X i =1 z > i L † b e b > e L † z i ≈ (cid:15)/ Tr (cid:16) L † b e b > e L † (cid:17) . (29)By Lemma 5.3, we have Tr (cid:16) L † b e b > e L † (cid:17) ≥ n U , and hence 1 M M X i =1 z > i L † b e b > e L † z i ≥ exp( − (cid:15)/
3) 2 n U ≥ n U , (30)where the second inequality follows by 0 < (cid:15) ≤ / δ = (cid:15)n − U − when invoking LaplSolve , by Lemma 5.1, (cid:13)(cid:13)(cid:13) y i − L † z i (cid:13)(cid:13)(cid:13) L ≤ (cid:15)n − U − (cid:13)(cid:13)(cid:13) L † z i (cid:13)(cid:13)(cid:13) L holds for each i . Then, by Lemma 5.2, we have that (cid:12)(cid:12)(cid:12) y > i b e b > e y i − z > i L † b e b > e L † z i (cid:12)(cid:12)(cid:12) ≤ (cid:15)n − U − holds for each i . We then have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M M X i =1 y > i b e b > e y i − M M X i =1 z > i L † b e b > e L † z i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ M M X i =1 (cid:12)(cid:12)(cid:12) y > i b e b > e y i − z > i L † b e b > e L † z i (cid:12)(cid:12)(cid:12) ≤ (cid:15)n − U − ≤ (cid:15) M M X i =1 z > i L † b e b > e L † z i ! , where the last inequality follows by (30). Thus,(1 − (cid:15)/
6) 1 M M X i =1 z > i L † b e b > e L † z i ≤ M M X i =1 y > i b e b > e y i ≤ (1 + (cid:15)/
6) 1 M M X i =1 z > i L † b e b > e L † z i , which implies 1 M M X i =1 z > i L † b e b > e L † z i ≈ (cid:15)/ M M X i =1 y > i b e b > e y i . (31)By Lemma 5.4, we have ˆ r e ≈ θ(cid:15)/ R eff ( e ), which by Lemma 5.5 implies1 − (1 − θ ) w ( e )ˆ r e ≈ (cid:15)/ − (1 − θ ) w ( e ) R eff ( e ) . (32)Combining Equations (29), (31) and (32), we have M M P i =1 y > i b e b > e y i − (1 − θ ) w ( e )ˆ r e ≈ (cid:15) Tr (cid:16) L † b e b > e L † (cid:17) − (1 − θ ) w ( e ) R eff ( e ) , which together with Equation (27) proves this theorem. θ -Kirchhoff Vertex Centrality We now combine the projection based approximation algorithm from Section 5 with the recursiveSchur complement approximation algorithm to produce a routine for estimating Kirchhoff vertexcentrality as defined in Definition 2.9 in nearly-linear time.25 .1 Turning to Low-rank Updates
We will treat the θ -deletion of a vertex as θ -deleting a batch of edges from the graph, which in turncorresponds to a high rank update to the graph Laplacian. Specifically, we can define the matrix B E v as the deg u ( v ) × n edge-vertex incidence matrix containing the edges incident to v , and W E v asthe corresponding deg u ( v ) × deg u ( v ) diagonal edge weight matrix. Here we use E v def = { ( u, v ) | u ∼ v } to denote the set of edges incident with v , and deg u ( v ) def = | E v | to denote the number of edgesincident with v . Then the graph Laplacian with edges in E v θ -deleted is L \ θ E v = L − (1 − θ ) B > E v W E v B E v . By Lemma 4.3, our goal becomes solving Problem 4 on the difference between the pseudoinverses ofthese matrices. Specifically, computing the value of z > (cid:18)(cid:16) L − (1 − θ ) B > E v W E v B E v (cid:17) † − L † (cid:19) z for a vector z . For this we once again turn to low-rank updates, specifically the Woodbury formula. Lemma 6.1 (Derived from Woodbury formula) . Given an edge set T ⊂ E supported on vertex set C , and a scalar < θ < . Let B T be the | T | × n edge-vertex incidence matrix corresponding toedges in T , and W T be the | T | × | T | diagonal edge weight matrix corresponding to edges in T . Thefollowing statement holds: ( L \ θ T ) † = L † + (1 − θ ) L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T L † . (33)Since the off-the-shelf Woodbury formula is for full rank matrices, we give the detailed proof ofEquation (33) in Appendix A.Note that this formula applies to any subset of T edges. The only property of evaluatingKirchhoff vertex centrality we need is that the total size of such T s over all vertices is 2 m .This means just as in Section 5, the problem reduces to estimating z > L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T L † z . Furthermore, since we can compute L † z to high accuracy via a single solver to a linear system in agraph Laplacian, and W / T B T is | T | × n matrix with 2 | T | nonzero entries, we can compute foreach set T the vector W / T B T L † z in O ( | T | ) time (after ˜ O ( m ) preprocessing time to compute anapproximation to L † z ). To track the error for the solver, we will need to following two lemmas,which we prove in Appendix C.2. Lemma 6.2.
Let L be the Laplacian of a graph with all weights in the range [1 , U ] , and z beany vector such that k z k ≤ n . Suppose y is a vector such that (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) L ≤ δ (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L for some < δ < . For any edge set T ⊂ E , we have: (cid:12)(cid:12)(cid:12)(cid:12) z > L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T L † z − y > B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T y (cid:12)(cid:12)(cid:12)(cid:12) ≤ θ − δn U . (34) Lemma 6.3.
Let L be the Laplacian of a graph with all weights in the range [1 , U ] . For any edgeset T ⊂ E of the graph, we have Tr (cid:18) L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T L † (cid:19) ≥ | T | n U . .2 Approximating Quadratic Forms Since we can utilize nearly-linear time solvers for Laplacian linear systems to compute high accuracyapproximations to the vector L † z , the problem is further reduced to estimating quadratic forms of (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − . Since the edges in T form a subgraph of L , we have B > T W T B (cid:22) L , and in turn0 (cid:22) W / T B T L † B > T W / T (cid:22) I This coupled with the assumption that 0 < θ means that the eigenvalues of matrix(1 − θ ) W / T B T L † B > T W / T are bounded away from 1. Therefore, we can use iterative methodsto solve the resulting system. As we work entirely with matrix approximations, we will use thefollowing matrix-based version of Chebyshev iteration. More details on these iterative methods canbe found in Section 11.2 of [GVL12]. Lemma 6.4 (Chebyshev iteration) . There is an algorithm
ChebSolve ( P , κ, (cid:15), b ) such that for anypositive definite matrix P along with κ such that κ I (cid:22) P (cid:22) I , ChebSolve ( P , κ, (cid:15), b ) correspondsto a linear operator on b such that the matrix Z ChebSolve realizing this operator satisfies Z ChebSolve ≈ (cid:15) P − , and the cost of the algorithm is O ( √ κ log(1 /(cid:15) )) matrix-vector multiplications involving P . Therefore the main difficulty becomes finding the matrix W / T B T L † B > T W / T . Note that while B T has up to n columns, most of these column are 0s. So it means that we canonly consider the entries corresponding to V ( T ), the set of vertices incident to at least one edge in T , using the following fact about Schur complements. Fact 6.5.
Let L be a Laplacian matrix, and C be a subset of vertices. Then, we have (cid:16) L † (cid:17) CC = Sc ( L , C ) † . However, we only have approximate Schur complements. To bound this also, we once againinvoke the bound about preservations of approximations when subtracting matrices from I fromLemma 5.5. Lemma 6.6.
There is an algorithm e x = QuadSolve ( B T , W T , b , θ, (cid:15), e S ) which takes an edge-vertex incidence matrix B T corresponding to edges in T ⊂ E with edge weight matrix W T supportedon vertex set V ( T ) , a vector b ∈ R n , scalars < θ ≤ / and < (cid:15) < / , and a Laplacian matrix e S whose edges are supported on V ( T ) such that e S ≈ (cid:15)θ/ Sc ( L , V ( T )) , and returns a value e x satisfying e x ≈ (cid:15) b > (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − b . The algorithm runs in time O (nnz( e S ) θ − . log n log(1 /(cid:15) )+ | T | θ − . (cid:15) − log n log(1 /(cid:15) ) polyloglog( n )) ,where nnz( e S ) is the number of nonzero entries in e S . roof. We will invoke preconditioned Chebyshev iteration as stated in Lemma 6.4 to estimate thequantity b > (cid:18) I − (1 − θ ) W / T B T,V ( T ) e S † B > T,V ( T ) W / T (cid:19) − b . Since B T is only non-zero on the entries corresponding to V ( T ), Fact 6.5 gives W / T B T L † B > T W / T = W / T B T,V ( T ) Sc ( L , V ( T )) † B > T,V ( T ) W / T , and hence Fact 2.1 Part 10 gives W / T B T L † B > T W / T ≈ (cid:15)θ/ W / T B T,V ( T ) e S † B > T,V ( T ) W / T . (35)Also, since T is a subset of edges, (cid:16) W / T B T (cid:17) > (cid:16) W / T B T (cid:17) (cid:22) L , which in turn implies W / T B T L † B > T W / T (cid:22) I , and θ I (cid:22) I − (1 − θ ) W / T B T,V ( T ) Sc ( L , V ( T )) † B > T,V ( T ) W / T (cid:22) I . Combining this with the approximation factor above from Equation (35) and Lemma 5.5 thengives I − (1 − θ ) W / T B T L † B > T W / T ≈ (cid:15)/ I − (1 − θ ) W / T B T,V ( T ) e S † B > T,V ( T ) W / T . (36)To apply e S † , we can invoke the algorithm ApxPartialCholesky ( e S , { v } , (cid:15)θ/
9) in Lemma 3.5 foran arbitrary vertex v to get an (cid:15)θ/ complete Cholesky factorization of e S andthen apply its inverse quickly. Suppose the Cholesky factorization returned is e L f D e L > ≈ (cid:15)θ/ e S , thenagain by Lemma 5.5 we have I − (1 − θ ) W / T B T,V ( T ) e S † B > T,V ( T ) W / T ≈ (cid:15)/ I − (1 − θ ) W / T B T,V ( T ) (cid:18) e L f D e L > (cid:19) † B > T,V ( T ) W / T . (37)Combining Equation (36) and (37) leads to I − (1 − θ ) W / T B T L † B > T W / T ≈ (cid:15)/ I − (1 − θ ) W / T B T,V ( T ) (cid:18) e L f D e L > (cid:19) † B > T,V ( T ) W / T , (38)which also means that all the eigenvalues of I − (1 − θ ) W / T B T,V ( T ) (cid:18) e L f D e L > (cid:19) † B > T,V ( T ) W / T arebetween exp( − (cid:15)/ θ and 1. Therefore by Lemma 6.4, we can access a linear operator Z Solve suchthat Z ChebSolve ≈ (cid:15)/ I − (1 − θ ) W / T B T,V ( T ) (cid:18) e L f D e L > (cid:19) † B > T,V ( T ) W / T ! − , (39)28hose cost is O ( θ − . log(1 /(cid:15) )) matrix-vector multiplications involving I − (1 − θ ) W / T B T,V ( T ) (cid:18) e L f D e L > (cid:19) † B > T,V ( T ) W / T . Here I , W / T and B T,V ( T ) can all be appliedin O ( | T | ) time. By Lemma 3.5, (cid:18) e L f D e L > (cid:19) † can be applied in O (nnz( e S ) + | T | θ − (cid:15) − log n ) time,and ApxPartialCholesky ( e S , { v } , (cid:15)θ/
9) runs in O (nnz( e S ) log n + | T | θ − (cid:15) − log n polyloglog( n ))time.Inverting both sides of Equation (38) and then combining it with Equation (39) gives Z ChebSolve ≈ (cid:15) (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − , so we can set e x = b > Z ChebSolve b and return it as our overall estimate.Thus, the problem becomes efficiently approximating Schur complements onto subsets of edges.We give an algorithm QuadApprox that first computes approximate Schur complements ontoneighbors of each vertex and then uses the algorithm
QuadSolve in Lemma 6.6 to compute z > B > E v W / E v (cid:16) I − (1 − θ ) W / E v B E v L † B > E v W / E v (cid:17) − W / E v B E v z for some vector z . The pseudocode for QuadApprox is given in Algorithm 5. Note that in thepseudocode we use G [ C ] to denote G ’s induced graph on a subset of vertices C , z C to denote a | C | -dimensional vector obtained from z by taking entries corresponding to vertices in C , and deg u ( v )to denote the number of edges incident with v . The performance of QuadApprox is characterizedin Lemma 6.7.
Lemma 6.7.
Given a connected undirected graph G = ( V, E ) with n vertices, m edges, positiveedge weights w : E → R + , and associated Laplacian L , a set of vertices V Q ⊂ V such that V = n N ( v ) | v ∈ V Q o ∪ V Q , a vector z ∈ R n , and scalars < θ ≤ / , < (cid:15) ≤ / , the algorithm QuadApprox ( G, L , V Q , z , θ, (cid:15)θ/ , (cid:15) ) returns a set of pairs ˆ N ∆ = { ( v, ˆ n ∆ v ) | v ∈ V Q } . With highprobability, the following statement holds: For ∀ v ∈ V Q , n ∆ v ≈ (cid:15) ˆ n ∆ v , (40) where n ∆ v = z > B > E v W / E v (cid:16) I − (1 − θ ) W / E v B E v L † B > E v W / E v (cid:17) − W / E v B E v z , and E v = { ( u, v ) | u ∼ v } is the set of edges incident with v . The total running time of this algorithmis bounded by O ( m ( θ − (cid:15) − log n + θ − . (cid:15) − log n log(1 /(cid:15) )) polyloglog( n )) .Proof of Lemma 6.7. Let volume( V Q ) denote the quantity vol on Line 5. We first observe thatevery time we recursively call QuadApprox , one of the following two events occurs:1. volume( V Q ) becomes no more than its 3 / (cid:12)(cid:12)(cid:12) V Q (cid:12)(cid:12)(cid:12) becomes 1 (Lines 11).When (cid:12)(cid:12)(cid:12) V Q (cid:12)(cid:12)(cid:12) = 1, the algorithm will go to Lines 2 - 4, and hence the recursion depth is only 1. Then,as we set V Q = V in the earliest call to QuadApprox , we have that the total recursion depth isno more than log volume( V ) = log m . 29 lgorithm 5: QuadApprox ( G, S , V Q , z , θ, (cid:15) , (cid:15) ) Input : G = ( V, E ): A graph. S : A graph Laplacian whose edges are supported on V . V Q ⊂ V : a set of vertices z ∈ R | V | : a vector. θ : a scalar between 0 and 1/2. (cid:15) : the error parameter for Schur complement. (cid:15) : the error parameter for QuadSolve . Output : ˆ N ∆ = { ( v, ˆ n ∆ v ) | v ∈ V Q } . if (cid:12)(cid:12)(cid:12) V Q (cid:12)(cid:12)(cid:12) = 1 then Let n and m be the number of vertices and edges in G , respectively. Let B be the m × n edge-vertex incidence matrix of G , and W be the m × m diagonaledge weight matrix of G . Let ˆ n ∆ v = QuadSolve ( B , W , W / Bz , θ, (cid:15) , S ) and return n ( v, ˆ n ∆ v ) o for the only vertex v ∈ V Q . Let vol = P v ∈ V Q deg u ( v ), and set (cid:15) schur = (cid:15) / log vol. Let V be vertices in V Q with deg u ( v ) ≥ vol / if V = ∅ then for each v ∈ V do Let C denote v and its neighbors. ( e L , f D , e S ) ← ApxPartialCholesky ( S , C, (cid:15) schur ) ˆ N ∆( v ) ← QuadApprox ( G [ C ] , e S , { v } , z C , θ, (cid:15) − (cid:15) schur , (cid:15) ) Let C denote vertices in V Q \ V and their neighbors. ( e L , f D , e S ) ← ApxPartialCholesky ( S , C, (cid:15) schur ) return QuadApprox ( G [ C ] , e S , V Q \ V , z C , θ, (cid:15) − (cid:15) schur , (cid:15) ) ∪ (cid:16)S v ∈ V ˆ N ∆( v ) (cid:17) Divide V Q into two parts V (1) and V (2) such that both P v ∈ V (1) deg u ( v ) and P v ∈ V (2) deg u ( v )are in the range h vol , vol i . for i = 1 to do Let C denote vertices in V ( i ) and their neighbors. ( e L , f D , e S ) ← ApxPartialCholesky ( S , C, (cid:15) schur ) ˆ N ∆( i ) ← QuadApprox ( G [ C ] , e S , V ( i ) , z C , θ, (cid:15) − (cid:15) schur , (cid:15) ) return ˆ N ∆(1) ∪ ˆ N ∆(2) . 30e then give guarantees for our approximations. Note that we set (cid:15) schur = (cid:15) / log vol (Line 5),and when recursively calling QuadApprox we set the (cid:15) of the recursive call to (cid:15) − (cid:15) schur (Line 11, 14and 19). Then, since we set (cid:15) = (cid:15)θ/ QuadApprox , we have that (cid:15) schur ≤ (cid:15)θ / log m always holds. Coupled with the fact that the total recursion depth is no more thanlog m , on Line 4 we have that S ≈ (cid:15)θ/ Sc ( L G , V ( E v ))holds for the only vertex v ∈ V Q , where L G is the Laplacian matrix of the graph in the earliest callto QuadApprox (i.e., the original graph). Then by Lemma 6.6, ˆ n ∆ v on Line 4 satisfiesˆ n ∆ v ≈ (cid:15) z > B > E v W / E v (cid:16) I − (1 − θ ) W / E v B E v L † B > E v W / E v (cid:17) − W / E v B E v z . We now analyze the running time.Let T ( m ) denote the running time of QuadApprox ( G, S , V Q , z , θ, (cid:15) , (cid:15) ), where m = volume( V Q ).Let n cur and m cur denote the number of vertices and the number of edges in L in the current call,respectively. We first assume QuadSolve to be an O (1) operation, and hence T ( m ) = O (1) for (cid:12)(cid:12)(cid:12) V Q (cid:12)(cid:12)(cid:12) = 1. For (cid:12)(cid:12)(cid:12) V Q (cid:12)(cid:12)(cid:12) >
1, We consider the set V on Line 6:1. If V is not empty, the algorithm goes to Lines 7 - 14. Since there are at most 4 vertices in V ,and by our assumption the recursive calls to QuadApprox on Line 11 all run in O (1) time, wehave by Lemma 3.5 Lines 7 - 13 runs in total O (( m cur log n + n cur (cid:15) − log n ) polyloglog( n ))time. Since V is not empty, we have volume( V \ V ) ≤ volume( V ). Hence, the running timeof the recursive call to QuadApprox on Line 14 is at most T (3 m/ V is empty, the algorithm goes to Lines 15 - 20. By Lemma 3.5, the calls to ApxPartialCholesky on Line 18 run in total O (( m cur log n + n cur (cid:15) − log n ) polyloglog( n ))time. The running time of the recursive calls to QuadApprox on Line 19 is T ( m )+ T ( m − m ),where m = volume( V (1) ) is in the range [ m/ , m/ (cid:15) schur = O ( θ(cid:15)/ log n ), n cur = O ( m ),‘ and by Lemma 3.5 m cur = O ( n cur (cid:15) − log n ) = O ( mθ − (cid:15) − log n ), we have in the worst case T ( m ) = T (3 m/
4) + T ( m/
4) + O ( mθ − (cid:15) − log n polyloglog( n )) , which gives T ( m ) = O ( mθ − (cid:15) − log n polyloglog( n )).Note that we get this running time under the assumption that QuadSolve is an O (1) operation.Thus, we also need to analyze the total running time of the calls to QuadSolve on Line 4.By Lemma 6.6, the
QuadSolve ( B , W , W / Bz , θ, (cid:15) , S ) on Line 4 runs in O (nnz( S ) θ − . log n log(1 /(cid:15) )+deg u ( v ) θ − . (cid:15) − log n log(1 /(cid:15) ) polyloglog( n )) time, where v indicatesthe only vertex in V Q and deg u ( v ) is the number of edges incident to v . By Lemma 3.5, we havennz( S ) = O (deg u ( v ) (cid:15) − log n ) = O (deg u ( v ) θ − (cid:15) − log n ). Then, summing this running timeover all vertices gives O ( mθ − . (cid:15) − log n log(1 /(cid:15) ) polyloglog( n )), which plus T ( m ) gives the overallrunning time of this algorithm. C ∆ θ ( v ) We give the pseudocode of the algorithm
VertexCentComp which approximates θ -Kirchhoffvertex centrality C ∆ θ ( v ) for all v ∈ V in Algorithm 6. Note that in this algorithm we once againinvoke the Laplacian solver of [CKM + VertexCentComp is characterizedin Theorem 1.3. Analyzing this algorithm gives the main result for estimating vertex centralities.31 lgorithm 6:
VertexCentComp ( G = ( V, E ) , w, θ, (cid:15) ) Input : G = ( V, E ), w : A connected undirected graph with positive edgesweights w : E → R + . θ : A scalar between 0 and 1/2. (cid:15) : Error of the centrality estimate per vertex. Output : ˆ C ∆ = n ( v, ˆ c ∆ v ) | v ∈ V o . Let z , . . . , z M be independent random ± M = (cid:6) (cid:15) − ln(2 n ) (cid:7) . for i = 1 to M do y i ← LaplSolve ( L G , z i , θ(cid:15)n − U − ) (cid:16) ˆ N ∆( i ) = { ( v, ˆ n ∆( i ) v ) | v ∈ V } (cid:17) ← QuadApprox ( L G , V, y i , θ, θ(cid:15)/ , (cid:15)/ For each v ∈ V compute ˆ c ∆ v = (1 − θ ) nM M P i =1 ˆ n ∆( i ) v and return ˆ C ∆ = n ( v, ˆ c ∆ v ) | v ∈ V o . Proof of Theorem 1.3.
The running time is the total cost of O ( (cid:15) − log n ) calls to LaplSolve eachof which runs in O ( m log . n log( (cid:15)θ ) polyloglog( n )) time, and O ( (cid:15) − log n ) calls to QuadApprox each of which runs in O ( m ( θ − (cid:15) − log n + θ − . (cid:15) − log n log(1 /(cid:15) )) polyloglog( n )) time.In the rest of this proof, we will use the matrix C v , defined as C v def = B > E v W / E v (cid:16) I − (1 − θ ) W / E v B E v L † B > E v W / E v (cid:17) − W / E v B E v , to simplify notation.Since M = (cid:6) (cid:15) − ln(2 n ) (cid:7) ≥
48 ( (cid:15)/ − ln(2 n ), by Lemma 4.3, we have1 M M X i =1 z > i L † C v L † z i ≈ (cid:15)/ Tr (cid:16) L † C v L † (cid:17) . (41)By Lemma 6.3, we have Tr (cid:16) L † C v L † (cid:17) ≥ n U , and hence 1 M M X i =1 z > i L † C v L † z i ≥ exp( − (cid:15)/
3) 2 n U ≥ n U , (42)where the second inequality follows by 0 < (cid:15) ≤ / δ = θ(cid:15)n − U − when invoking LaplSolve , by Lemma 5.1, (cid:13)(cid:13)(cid:13) y i − L † z i (cid:13)(cid:13)(cid:13) L ≤ θ(cid:15)n − U − (cid:13)(cid:13)(cid:13) L † z i (cid:13)(cid:13)(cid:13) L , holds for each i . Then, by Lemma 6.2, we have that (cid:12)(cid:12)(cid:12) y > i C v y i − z > i L † C v L † z i (cid:12)(cid:12)(cid:12) ≤ (cid:15)n − U − i . We then have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M M X i =1 y > i C v y i − M M X i =1 z > i L † C v L † z i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ M M X i =1 (cid:12)(cid:12)(cid:12) y > i C v y i − z > i L † C v L † z i (cid:12)(cid:12)(cid:12) ≤ (cid:15)n − U − ≤ (cid:15) M M X i =1 z > i L † C v L † z i ! , where the last inequality follows by (42). Thus,(1 − (cid:15)/
6) 1 M M X i =1 z > i L † C v L † z i ≤ M M X i =1 y > i C v y i ≤ (1 + (cid:15)/
6) 1 M M X i =1 z > i L † C v L † z i , which implies 1 M M X i =1 z > i L † C v L † z i ≈ (cid:15)/ M M X i =1 y > i C v y i . (43)By Lemma 6.7, we have ˆ n ∆( i ) v ≈ (cid:15)/ y > i C v y i . (44)Combining Equation (41), (43) and (44), we have1 M M X i =1 ˆ n ∆( i ) v ≈ (cid:15) Tr (cid:16) L † C v L † (cid:17) , which coupled with the fact that C ∆ θ ( v ) = K ( G \ θ E v ) − K ( G ) by definition= n (cid:16) Tr (cid:16) ( L \ θ E v ) † (cid:17) − Tr (cid:16) L † (cid:17)(cid:17) by Fact 2.4= n (1 − θ )Tr (cid:16) L † C v L † (cid:17) by Equation (33)gives the guarantee of our approximation. The Kirchhoff index arises in many applications such as noisy consensus problems [PB14] and socialrecommender systems [WLC16]. It is a global index, and any changes of network structure, e.g.weight of edges, can be reflected in this popular index. In this paper, we proposed to use Kirchhoffindex as a global metric of the importance of edges in weighted undirected networks. For anynetwork, when the weight of any edge e is changed from w ( e ) from θw ( e ), the Kirchhoff index ofthe resulting graph will strictly increase, with the increase deciphering the importance of edge e .We used the Kirchhoff index of the new graph, or its increment with respect to the original graph,33s the θ -Kirchhoff edge centrality. We demonstrated experimentally that this new global measure ofedge centrality has a more discriminating power than edge betweenness, spanning edge centrality,and current-flow centrality.However, the time cost of exactly computing the θ -Kirchhoff edge centrality is prohibitive. Toovercome this weakness, we introduced two approaches that estimate the θ -Kirchhoff edge centralityfor all edges in nearly linear time. Our proposed centrality metrics are the first global measure ofcentrality that can be estimated in nearly linear time. Our algorithms combine techniques fromseveral recent works on graph algorithms [LSW15, DKP + θ -Kirchhoff vertex centrality, as well as estimating the Kirchhoffedge centrality to a set of edges. This raises the possibility of designing highly efficient algorithmsthat can detect the set of k most influential edges, that is, the k edges whose θ -deletion leads to thelargest increase of the Kirchhoff index.Despite the advantages of our algorithms, their theoretical performance still has much room forimprovement , both in the overhead of logarithmic factors and the dependency on θ . The latteris particularly interesting because our two algorithms for estimating edge centrality can performbetter under different regimes of θ . On the other hand, the importance of centrality measuresin graph mining means it is just as, if not more, interesting to study the practical behaviors ofour algorithms. Specifically, to see if they are reasonably fast and accurate on massive networkswith millions of vertices and edges. Recent packages for solving large scale linear systems andrelated tasks [LB12, KMT11, SSM14, Spi17] should greatly facilitate such a study. Moreover, thesignificantly higher deviations from our experiments suggest the question of whether it is possible totheoretically model the advantages/disadvantages of the many centrality measures.Finally, it should be mentioned that as an application of the introduced edge centrality, westudied the vertex centrality based on the idea of the definition for C ∆ θ ( e ). Actually, we can alsodefine the centrality of a vertex v as the Kirchhoff index of the graph G \ θ E v , the algorithm for the (cid:15) -approximation of which is similar to EdgeCentComp1 . We thus omit the algorithmic details ofthis version of vertex centrality for the lack of space. Another reason for ignoring this algorithm isthat our main focus is the edge centrality. 34 eferences [Ach01] Dimitris Achlioptas. Database-friendly random projections. In
Proceedings of the20th ACM SIGACT-SIGMOD-SIGART Symposium on Principles of Database Systems(PODS) , 2001.[ADK +
16] Ittai Abraham, David Durfee, Ioannis Koutis, Sebastian Krinninger, and Richard Peng.On fully dynamic graph sparsifiers. In
Proceedings of IEEE 57th Annual Symposiumon Foundations of Computer Science (FOCS) , pages 335–344, 2016.[AT11] Haim Avron and Sivan Toledo. Randomized algorithms for estimating the trace ofan implicit symmetric positive semi-definite matrix.
Journal of ACM , 58(2):8:1–8:34,2011.[BCH14] Daniel Bienstock, Michael Chertkov, and Sean Harnett. Chance-constrained optimalpower flow: Risk-aware network control under uncertainty.
SIAM Review , 56(3):461–495, 2014.[BDFMR16] Francesco Bonchi, Gianmarco De Francisci Morales, and Matteo Riondato. Centralitymeasures on big graphs: Exact, approximated, and distributed algorithms. In
Proceed-ings of the 25th International Conference Companion on World Wide Web (WWW) ,pages 1017–1020, 2016.[BF05] Ulrik Brandes and Daniel Fleischer. Centrality measures based on current flow. In
Proceedings of the 22nd Annual Symposium on Theoretical Aspects of Computer Science(STACS) , volume 3404, pages 533–544, 2005.[BHNT15] Sayan Bhattacharya, Monika Henzinger, Danupon Nanongkai, and CharalamposTsourakakis. Space-and time-efficient algorithm for maintaining dense subgraphs onone-pass dynamic streams. In
Proceedings of the 47th annual ACM Symposium onTheory of Computing (STOC) , pages 173–182, 2015.[BKMM07] David A Bader, Shiva Kintali, Kamesh Madduri, and Milena Mihail. Approximat-ing betweenness centrality. In
Proceedings of the 5th International Conference onAlgorithms and Models for the Web-Graph (WAW) , volume 4863, pages 124–137, 2007.[BLHL01] Tim Berners-Lee, James Hendler, and Ora Lassila. The semantic web.
ScientificAmerican , 284(5):28–37, 2001.[BP07] Ulrik Brandes and Christian Pich. Centrality estimation in large networks.
Interna-tional Journal of Bifurcation and Chaos , 17(07):2303–2318, 2007.[Bra01] Ulrik Brandes. A faster algorithm for betweenness centrality.
Journal of MathematicalSociology , 25(2):163–177, 2001.[BV14] Paolo Boldi and Sebastiano Vigna. Axioms for centrality.
Internet Mathematics ,10(3-4):222–262, 2014.[BWLM16] Elisabetta Bergamini, Michael Wegner, Dimitar Lukarski, and Henning Meyerhenke.Estimating current-flow closeness centrality with a multigrid laplacian solver. In
Proceedings of the 7th SIAM Workshop on Combinatorial Scientific Computing (CSC) ,pages 1–12, 2016. 35CKM +
14] Michael B. Cohen, Rasmus Kyng, Gary L. Miller, Jakub W. Pachocki, Richard Peng,Anup Rao, and Shen Chen Xu. Solving SDD linear systems in nearly m log / n time.In Proceedings of the 46th annual ACM Symposium on Theory of Computing (STOC) ,pages 343–352, 2014.[DKP +
17] David Durfee, Rasmus Kyng, John Peebles, Anup B. Rao, and Sushant Sachdeva.Sampling random spanning trees faster than matrix multiplication. In
Proceedings ofthe 49th annual ACM Symposium on Theory of Computing (STOC) , pages 730–742,2017.[DPPR17] David Durfee, John Peebles, Richard Peng, and Anup B. Rao. Determinant-preservingsparsification of SDDM matrices with applications to counting and sampling spanningtrees.
CoRR , abs/1705.00985, 2017.[DSD +
11] Li Ding, Dana Steil, Brandon Dixon, Allen Parrish, and David Brown. A relationcontext oriented approach to identify strong ties in social networks.
Knowledge-BasedSystems , 24(8):1187–1195, 2011.[ESVM +
11] W Ellens, FM Spieksma, P Van Mieghem, A Jamakovic, and RE Kooij. Effectivegraph resistance.
Linear Algebra and its Applications , 435(10):2491–2506, 2011.[FQY14] Minyu Feng, Hong Qu, and Zhang Yi. Highest degree likelihood search algorithmusing a state transition matrix for complex networks.
IEEE Transactions on Circuitsand Systems I: Regular Papers , 61(10):2941–2950, 2014.[GN02] Michelle Girvan and Mark EJ Newman. Community structure in social and biologicalnetworks.
Proceedings of the National Academy of Sciences , 99(12):7821–7826, 2002.[GSS08] Robert Geisberger, Peter Sanders, and Dominik Schultes. Better approximation ofbetweenness centrality. In
Proceedings of the Meeting on Algorithm Engineering &Expermiments , pages 90–100, 2008.[GVL12] Gene H Golub and Charles F Van Loan.
Matrix computations , volume 3. JHU Press,2012.[HAY16] Takanori Hayashi, Takuya Akiba, and Yuichi Yoshida. Efficient algorithms for spanningtree centrality. In
Proceedings of the 25th International Joint Conference on ArtificialIntelligence (IJCAI) , pages 3733–3739, 2016.[Hut89] MF Hutchinson. A stochastic estimator of the trace of the influence matrix for Lapla-cian smoothing splines.
Communications in Statistics-Simulation and Computation ,18(3):1059–1076, 1989.[HX16] Nicholas JA Harvey and Keyulu Xu. Generating random spanning trees via fastmatrix multiplication. In
Proceedings of Latin American Symposium on TheoreticalInformatics , pages 522–535, 2016.[JL84] William B Johnson and Joram Lindenstrauss. Extensions of Lipschitz mappings intoa Hilbert space.
Contemporary Mathematics , 26(189-206):1, 1984.[KLP +
16] Rasmus Kyng, Yin Tat Lee, Richard Peng, Sushant Sachdeva, and Daniel A. Spielman.Sparsified Cholesky and multigrid solvers for connection Laplacians. In
Proceedings of he 48th annual ACM Symposium on Theory of Computing (STOC) , pages 842–850,2016.[KMT11] Ioannis Koutis, Gary L. Miller, and David Tolliver. Combinatorial preconditioners andmultilevel solvers for problems in computer vision and image processing. ComputerVision and Image Understanding , 115(12):1638–1646, 2011.[Knu93] Donald Ervin Knuth.
The Stanford GraphBase: a Platform for Combinatorial Com-puting , volume 37. Addison-Wesley Reading, 1993.[KP17] John Kallaugher and Erie Price. A hybrid sampling scheme for triangle counting.In
Proceedings of the 28th Annual ACM-SIAM Symposium on Discrete Algorithms(SODA) , pages 1778–1797, 2017.[KPPS17] Rasmus Kyng, Jakub Pachocki, Richard Peng, and Sushant Sachdeva. A framework foranalyzing resparsification algorithms. In
Proceedings of the 28th Annual ACM-SIAMSymposium on Discrete Algorithms (SODA) , pages 2032–2043, 2017.[KR93] Douglas J Klein and Milan Randić. Resistance distance.
Journal of MathematicalChemistry , 12(1):81–95, 1993.[KS16] Rasmus Kyng and Sushant Sachdeva. Approximate gaussian elimination for lapla-cians - fast, sparse, and simple. In
Proceedings of IEEE 57th Annual Symposium onFoundations of Computer Science (FOCS) , pages 573–582, 2016.[LB12] Oren E. Livne and Achi Brandt. Lean algebraic multigrid (LAMG): fast graph laplacianlinear solver.
SIAM Journal on Scientific Computing , 34(4), 2012.[LCR +
16] Linyuan Lü, Duanbing Chen, Xiao-Long Ren, Qian-Ming Zhang, Yi-Cheng Zhang, andTao Zhou. Vital nodes identification in complex networks.
Physics Reports , 650:1–63,2016.[LM12] Amy N Langville and Carl D Meyer.
Who’s .Princeton University Press, 2012.[LS15] Yin Tat Lee and Aaron Sidford. Efficient inverse maintenance and faster algorithms forlinear programming. In
Proceedings of IEEE 56th Annual Symposium on Foundationsof Computer Science (FOCS) , pages 230–249, 2015.[LS17] Yin Tat Lee and He Sun. An SDP-based algorithm for linear-sized spectral sparsifica-tion. In
Proceedings of the 49th Annual ACM Symposium on Theory of Computing(STOC) , pages 678–687, 2017.[LSB +
03] David Lusseau, Karsten Schneider, Oliver J Boisseau, Patti Haase, Elisabeth Slooten,and Steve M Dawson. The bottlenose dolphin community of doubtful sound featuresa large proportion of long-lasting associations.
Behavioral Ecology and Sociobiology ,54(4):396–405, 2003.[LSW15] Yin Tat Lee, Aaron Sidford, and Sam Chiu-wai Wong. A faster cutting plane methodand its implications for combinatorial and convex optimization. In
Proceedings ofIEEE 56th Annual Symposium on Foundations of Computer Science (FOCS) , pages1049–1065, 2015. 37MGLKT15] Charalampos Mavroforakis, Richard Garcia-Lebron, Ioannis Koutis, and EvimariaTerzi. Spanning edge centrality: Large-scale computation and applications. In
Proceedings of the 24th International Conference on World Wide Web (WWW) , pages732–742, 2015.[New06] Mark EJ Newman. Finding community structure in networks using the eigenvectorsof matrices.
Physical review E , 74(3):036104, 2006.[New10] Mark Newman.
Networks: An introduction . Oxford university press, 2010.[PB14] Stacy Patterson and Bassam Bamieh. Consensus and coherence in fractal networks.
IEEE Transactions on Control of Network Systems , 1(4):338–348, 2014.[PS14] Richard Peng and Daniel A Spielman. An efficient parallel solver for SDD linearsystems. In
Proceedings of the 46th annual ACM Symposium on Theory of Computing(STOC) , pages 333–342, 2014.[San04] Piotr Sankowski. Dynamic transitive closure via dynamic matrix inverse. In
Proceedingsof 45th Annual IEEE Symposium on Foundations of Computer Science (FOCS) , pages509–517, 2004.[SM50] Jack Sherman and Winifred J Morrison. Adjustment of an inverse matrix correspondingto a change in one element of a given matrix.
The Annals of Mathematical Statistics ,21(1):124–127, 1950.[Spi17] Dan Spielman. Laplacians.jl. https://github.com/danspielman/Laplacians.jl ,2017.[SS11] Daniel A. Spielman and Nikhil Srivastava. Graph sparsification by effective resistances.
SIAM Journal of Computing , 40(6):1913–1926, 2011.[SSM14] Christian Staudt, Aleksejs Sazonovs, and Henning Meyerhenke. Networkit: Aninteractive tool suite for high-performance network analysis.
CoRR , abs/1403.3005,2014.[ST14] D. Spielman and S. Teng. Nearly linear time algorithms for preconditioning and solvingsymmetric, diagonally dominant linear systems.
SIAM Journal on Matrix Analysisand Applications , 35(3):835–885, 2014.[TMC +
13] Andreia Sofia Teixeira, Pedro T Monteiro, João A Carriço, Mário Ramirez, andAlexandre P Francisco. Spanning edge betweenness. In
Proceedings of 11st Workshopon Mining and Learning with Graphs , volume 24, pages 27–31, 2013.[Wil12] Virginia Vassilevska Williams. Multiplying matrices faster than Coppersmith-Winograd. In
Proceedings of the 44th annual ACM Symposium on Theory of Computing(STOC) , pages 887–898, 2012.[WLC16] Felix Ming Fai Wong, Zhenming Liu, and Mung Chiang. On the efficiency of socialrecommender networks.
IEEE/ACM Transactions on Networking , 24(4):2512–2524,2016.[WS98] Duncan J Watts and Steven H Strogatz. Collective dynamics of ‘small-world’ networks.
Nature , 393(6684):440–442, 1998. 38WS03] Scott White and Padhraic Smyth. Algorithms for estimating relative importance innetworks. In
Proceedings of the 9th ACM International Conference on KnowledgeDiscovery and Data Mining (KDD) , pages 266–275, 2003.[YTQ17] Keyou You, Roberto Tempo, and Li Qiu. Distributed algorithms for computation ofcentrality measures in complex networks.
IEEE Transactions on Automatic Control ,62(5):2080–2094, 2017.[Zac77] Wayne W Zachary. An information flow model for conflict and fission in small groups.
Journal of Anthropological Research , 33(4):452–473, 1977.39
Proofs of Our Version of Sherman-morrision and WoodburyFormulas
In this section, we give detailed proofs for the Sherman-Morrision and Woodbury formulas we used,i.e., Equations (26) and (33).In the proofs, we will use the matrix Π defined as Π def = LL † = I − n , where is the vector with all entries being 1. Proof of Equation (26).
First, we have b e (cid:16) − (1 − θ ) w ( e ) b > e L † b e (cid:17) = b e − (1 − θ ) w ( e ) b e b > e L † b e = (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) L † b e , (45)where the second equality follows by b e = Π b e = LL † b e .Since θ <
1, we have that L − (1 − θ ) w ( e ) b e b > e is a Laplacian matrix and 1 − (1 − θ ) w ( e ) b > e L † b e is strictly positive. Thus, Equation (45) implies (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † b e = L † b e − (1 − θ ) w ( e ) b > e L † b e . Then, we have L † = (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) L † = (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † (cid:16) Π − (1 − θ ) w ( e ) b e b > e L † (cid:17) = (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † − (1 − θ ) w ( e ) (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † b e b > e L † = (cid:16) L − (1 − θ ) w ( e ) b e b > e (cid:17) † − (1 − θ ) w ( e ) L † b e b > e L † − (1 − θ ) w ( e ) b > e L † b e , which implies Equation (26). Proof of Equation (33).
First, we have B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) = B > T W / T − (1 − θ ) B > T W T B T L † B > T W / T = (cid:16) L − (1 − θ ) B > T W T B T (cid:17) L † B > T W / T , (46)where the second equality follows by B > T = Π B > T = LL † B > T .Since θ <
1, we have that I − (1 − θ ) W / T B T L † B > T W / T is positive definite and L − (1 − θ ) B > T W T B T is a Laplacian matrix. Thus, Equation (46) implies (cid:16) L − (1 − θ ) B > T W T B T (cid:17) † B > T W / T = L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − . L † as (cid:16) L − (1 − θ ) B > T W T B T (cid:17) † (cid:16) L − (1 − θ ) B > T W T B T (cid:17) L † = (cid:16) L − (1 − θ ) B > T W T B T (cid:17) † (cid:16) Π − (1 − θ ) B > T W T B T L † (cid:17) = (cid:16) L − (1 − θ ) B > T W T B T (cid:17) † − (1 − θ ) (cid:16) L − (1 − θ ) B > T W T B T (cid:17) † B > T W T B T L † = (cid:16) L − (1 − θ ) B > T W T B T (cid:17) † − (1 − θ ) L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T L † , which implies Equation (33). B Approximations When Subtracted From Identity Matrix
In this section, we bound the transfer of approximations between A and B to approximationsbetween I − (1 − θ ) A and I − (1 − θ ) B . Proof of Lemma 5.5.
The given condition with the approximation can be written as:(1 − (cid:15) ) A (cid:22) B (cid:22) (1 + 2 (cid:15) ) A , which implies I − (1 − θ ) (1 + 2 (cid:15) ) A (cid:22) I − (1 − θ ) B (cid:22) I − (1 − θ ) (1 − (cid:15) ) A . Since 0 (cid:22) A (cid:22) I , we have the following lower bound: I − (1 − θ ) (1 + 2 (cid:15) ) A = (cid:18) − (cid:15)θ (cid:19) ( I − (1 − θ ) A ) + 2 (cid:15)θ I − (1 − θ ) (cid:18) (cid:15) + 2 (cid:15)θ (cid:19) A (cid:23) (cid:18) − (cid:15)θ (cid:19) ( I − (1 − θ ) A ) + (cid:20) (cid:15)θ − (1 − θ ) (cid:18) (cid:15) + 2 (cid:15)θ (cid:19)(cid:21) A . The coefficient on the trailing A in turn simplifies to 2 (cid:15) − (1 − θ )2 (cid:15) ≥ I − (1 − θ ) (1 − (cid:15) ) A = (cid:18) (cid:15)θ (cid:19) ( I − (1 − θ ) A ) − (cid:15)θ I + (1 − θ ) (cid:18) (cid:15) + 2 (cid:15)θ (cid:19) A (cid:22) (cid:18) (cid:15)θ (cid:19) ( I − (1 − θ ) A ) . Then the final bound involving exp(3 (cid:15)/θ ) follows from the condition of (cid:15)/θ being small.
C Error Tracking for Laplacian Solvers
In this section, we provide more details on error tracking for Laplacian solvers in Section 5 and 6 ina way similar to Section 4 of [SS11].We first give bounds on eigenvalues of L . Let L be the Laplacian matrix of a graph G = ( V, E )with n vertices, m edges and edge weights all in the range [1 , U ]. Let 0 = λ ≤ λ ≤ . . . ≤ λ n bethe eigenvalues of L , and 0 = ν ≤ ν ≤ . . . ≤ ν n be the eigenvalues of the normalized Laplacianmatrix, N def = D − / LD − / , of G . It is easy to verify that ν i ≤ λ i ≤ nU ν i holds for all i . Let φ G = min S ⊂ V | ∂ ( S ) | min ( d ( S ) , d ( V \ S ))41e the conductance of G , where | ∂ ( S ) | denotes the total weights of edges with one endpoint in S and the other endpoint in V \ S , and d ( S ) denotes the total degree of vertices in S . Then, we canbound λ by λ ≥ ν ≥ φ G / ≥ (cid:18) n U (cid:19) / , U ]= 12 n U . (47)We then bound λ n using the fact that L G (cid:22) U L K n , where K n is the complete graph of n vertices.Thus, λ Gn ≤ λ K n n U = nU. (48)From (47) and (48) it is immediate that12 n U Π (cid:22) L (cid:22) nU I and 1 nU Π (cid:22) L † (cid:22) n U I hold, where Π def = LL † = I − n > .We will also need to use the inequality (cid:12)(cid:12)(cid:12) x − y (cid:12)(cid:12)(cid:12) ≤ (2 | y | + | x − y | ) | x − y | for scalars x, y , which follows by (cid:12)(cid:12)(cid:12) x − y (cid:12)(cid:12)(cid:12) ≤ ( | x | + | y | ) | x − y | ≤ ( | y | + | y + ( x − y ) | ) | x − y | ≤ (2 | y | + | x − y | ) | x − y | . C.1 Error Tracking for the Laplacian Solver in Section 5
Proof of Lemma 5.2.
The lhs of inequality (28) can be written as (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12) . We first bound the value (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12) by (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) b e b > e by the triangle inequality of norms ≤ (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) L since b e b > e (cid:22) L ≤ δ (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L = δ p z > L † z ≤ δn . s z > L † zz > z since k z k ≤ n ≤ √ δn . U since L † ≤ n U I .42e then use the inequality (cid:12)(cid:12) x − y (cid:12)(cid:12) ≤ (2 | y | + | x − y | ) | x − y | to bound (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12) : (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:18) (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e + (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) (cid:12)(cid:12)(cid:12)(cid:12) k y k b e b > e − (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) b e b > e (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:16) (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L + √ δn . U (cid:17) √ δn . U by b e b > e (cid:22) L and the above bound ≤ (cid:16) √ n . U + √ δn . U (cid:17) √ δn . U since k z k ≤ n and L † ≤ n U I ≤ δn U by δ < Proof of Lemma 5.3. Tr (cid:16) L † b e b > e L † (cid:17) = b > e L † L † b e by cyclicness of trace= 2 b > e (cid:16) L † (cid:17) b e b > e b e since k b e k = 2 ≥ n U since (cid:16) L † (cid:17) (cid:23) n U Π . C.2 Error Tracking for the Laplacian Solver in Section 6
Proof of Lemma 6.2.
The lhs of inequality (34) can be seen as the difference between the followingtwo values: k y k B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T . By the triangle inequality of norms, the difference between the square roots of these two values is atmost (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T ≤ θ − . (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) B > T W T B T since (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − (cid:22) θ I ≤ θ − . (cid:13)(cid:13)(cid:13) y − L † z (cid:13)(cid:13)(cid:13) L since B > T W T B T (cid:22) L ≤ θ − . δ (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L = θ − . δ p z > L † z ≤ θ − . δn . s z > L † zz > z since k z k ≤ n ≤√ θ − . δn . U since L † ≤ n U I .43hen by the inequality (cid:12)(cid:12) x − y (cid:12)(cid:12) ≤ (2 | y | + | x − y | ) | x − y | , the lhs of (34) is at most (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T + √ θ − . δn . U √ θ − . δn . U ≤ (cid:16) θ − . (cid:13)(cid:13)(cid:13) L † z (cid:13)(cid:13)(cid:13) L + √ θ − . δn . U (cid:17) √ θ − . δn . U since B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T (cid:22) θ L ≤ θ (cid:16) √ n . U + √ δn . U (cid:17) √ δn . U since k z k ≤ n and L † (cid:22) n U I ≤ θ − δn U by δ < Proof of Lemma 6.3. Tr (cid:18) L † B > T W / T (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − W / T B T L † (cid:19) ≥ Tr (cid:16) L † B > T W T B T L † (cid:17) since (cid:16) I − (1 − θ ) W / T B T L † B > T W / T (cid:17) − (cid:23) I =Tr L † X e ∈ T w ( e ) b e b > e ! L † ! = X e ∈ T w ( e )Tr (cid:16) L † b e b > e L † (cid:17) ≥ | T | n U by Lemma 5.3 and w ( e ) ≥≥