Fully Dynamic Betweenness Centrality Maintenanceon Massive Networks
Takanori Hayashi

The University of Tokyo
Tokyo, 113-0033, Japan

thayashi@is.s.u-tokyo.ac.jp

Takuya Akiba

Yuichi Yoshidatakiba@nii.ac.jp

yyoshida@nii.ac.jp

National Institute of Informatics, National Institute of Informatics
CIT STAIR Lab & JST PRESTOTokyo, 101-8430, Japan

ABSTRACT

work analysis software such as SNAP [26], WebGraph [9],Gephi, NodeXL, and NetworkX.
Betweenness centrality measures the importance of eachvertex in terms of the number of shortest paths passing
through the vertex. For static settings, considerable effort
has been devoted to developing fast, exact, and approximate algorithms [3, 4, 10, 11, 18, 30, 32, 36]. Moreover, recent
approximation algorithms can accurately estimate the betweenness centrality of all vertices in near-linear time, which
is sufficiently fast even for today’s very large networks withbillions of edges.
However, certain real-world networks today are highly dynamic and are evolving rapidly. Examples include, but notlimited to, online social networks. Furthermore, an individual’s activities are often bursty with regard to temporal
locality [5]. Therefore, in order to realize and improve realtime network-aware applications (e.g., burst detection, trend
analysis, and socially-sensitive search ranking), instantly reflecting graph changes into centrality values is crucial.
In this regard, even if we use the recent sophisticated algorithms, the naive approach of re-running the static algorithms after each update would be too expensive for thelarge networks of today. Hence, considerable research is
being conducted to develop efficient dynamic methods [7,20, 22, 23, 25]. However, current dynamic methods are still
limited in terms of scalability. Exact dynamic methods require Ω(n2 ) preprocessing time and Ω(n2 ) data structure
size [20, 22, 23, 25]. The only dynamic approximate method
by Bergamini et al. [7] is semi-dynamic, i.e., it does not support any removal. Moreover, although its scalability is fairly
better than those of exact dynamic methods, it can barely
processes graphs with tens of millions of edges, and it cannot handle larger graphs with billions of edges. The method
by Bergamini et al. is based on shortest-path sampling byRiondato et al. [32].

Measuring the relative importance of each vertex in a network is one of the most fundamental building blocks in
network analysis. Among several importance measures, betweenness centrality, in particular, plays key roles in many
real applications. Considerable effort has been made fordeveloping algorithms for static settings. However, real networks today are highly dynamic and are evolving rapidly,
and scalable dynamic methods that can instantly reflect
graph changes into centrality values are required.
In this paper, we present the first fully dynamic method
for managing betweenness centrality of all vertices in a large
dynamic network. Its main data structure is the weighted
hyperedge representation of shortest paths called hypergraphsketch. We carefully design dynamic update procedure withtheoretical accuracy guarantee. To accelerate updates, we
further propose two auxiliary data structures called two-ballindex and special-purpose reachability index. Experimental
results using real networks demonstrate its high scalability
and efficiency. In particular, it can reflect a graph change
in less than a millisecond on average for a large-scale webgraph with 106M vertices and 3.7B edges, which is several
orders of magnitude larger than the limits of previous dynamic methods.

1.

INTRODUCTION

Measuring the relative importance of each vertex in a network is one of the most fundamental tasks in network analysis. Among several importance measures, betweenness centrality [2, 16], in particular, plays key roles in many applications, e.g., finding important actors in terrorist networks [13,24], estimating lethality in biological networks [14, 21], detecting hierarchical community structures [19, 29]and identifying crucial intersections in road networks [28]. Owing to
its popularity, almost all descriptive literature on network
analysis introduces betweenness centrality [28, 37]. Moreover, betweenness centrality is often adopted in various net-

1.1

Contributions

In this paper, we present the first fully dynamic approximate betweenness centrality algorithm, which maintains the
betweenness centrality of all vertices along with the internal data structures. It is fully dynamic in the sense that it
can instantly reflect insertions and deletions of vertices andedges. Our method has several orders of magnitude better
scalability, and unlike previous methods, it can handle networks with billions of edges, while maintaining a low error
rate. Moreover, even on very large networks, the proposed
method can update the internal data structure and centrality values of all the affected vertices in milliseconds, which isan order of magnitude faster than Bergamini’s method [7].

This work is licensed under the Creative Commons AttributionNonCommercial-NoDerivatives 4.0 International License. To view a copyof this license, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. For
any use beyond those covered by this license, obtain permission by emailing
info@vldb.org.
Proceedings of the VLDB Endowment, Vol. 9, No. 2
Copyright 2015 VLDB Endowment 2150-8097/15/10.

48

Our method is also based on sampling. However, unlike
Bergamini’s method [7], which selects a single shortest path
between each sampled pair of vertices, our focus is to maintain all the shortest paths between each sampled pair of vertices. Thereby, the necessary number of pairs of vertices are
saved. We represent the shortest paths in a data structure
named hypergraph sketch, which is inspired from a previous
static method [36]. The betweenness centrality of any vertex
can be easily obtained from the hypergraph sketch.
Accordingly, our main technical challenge is correctly andefficiently updating the hypergraph sketch. This is nontrivial because the part of the hypergraph sketch that requires modifications cannot be detected locally. Moreover,even if we overlook efficient updating, it is far from trivial
to correctly update the hypergraph sketch without reconstruction, because it is important to consider the probability
distribution of sampled pairs for accurate centrality estimation.
We summarize our main technical contributions for addressing these challenges as follows:

2.

RELATED WORK

2.1

Algorithms for Static Networks

As explained in the previous section, many algorithms
have been proposed for static settings. The asymptoticallyfastest exact betweenness centrality algorithm was proposedby Brandes [10]. This algorithm computes, for each vertex
v, the contribution of v to the centrality of all other vertices by conducting a single BFS and a backward aggregation. It runs in Θ(nm) time on unweighted networks andΘ(nm + n2 log n) time on weighted networks, and requires
Θ(n+m) space. There have been a few attempts to speedup
Brandes’ algorithm by exploiting the characteristics of complex networks [4,15,30], but these methods also need Θ(nm)time to compute betweenness centrality in the worst case.
Therefore, computing exact betweenness centrality on large
real-world networks is prohibitively expensive.
To handle large real-world networks, many approximate
algorithms have been proposed [3,11,18,32,36]. These methods are based on random sampling and allow us to estimate
betweenness centralites of vertices with acceptable accuracy
in almost linear time. Brandes and Pich [11] and Geisberger et al. [18] estimate the betweenness centralities bysampling a small set of vertices and accumulating their contributions to the betweenness centralities of other vertices.
Bader et al. [3] estimate the centrality of one specified vertex by adaptive sampling. Riondato and Kornaropoulos [32]approximate the centrality values by sampling peer-to-peer
shortest paths instead of vertices. Yoshida [36] estimatesthe centrality values by sampling vertex pairs. In algorithms
by [11, 36], in order to guarantee that the error of estimatedcentrality of every vertex is less than  with probability 1−δ,
we need to sample Ω( 12 log nδ ) vertices or vertex pairs. In
contrast, the method due to Riondato et al. only have tosample Ω( 12 log V D(G)) shortest paths, where V D(G) is theδdiameter of a graph G when we ignore the weight of edges.
We note that V D(G) is often much smaller than n on complex networks.

• Hypergraph Sketch: Our main data structure is theweighted hyperedge representation of shortest paths called
hypergraph sketch. We carefully design an update procedure that probabilistically replaces vertex pairs using
non-uniform distribution. We theoretically guaranteethat it ensures estimation accuracy against any sequenceof graph changes.
• Two-Ball Index: To efficiently detect the parts of hypergraph sketches that require modification, we propose
an auxiliary data structure named two-ball index. It basically maintains subtrees of shortest path trees and is designed to exploit small-world property of networks [27,34]
to achieve a small data structure size.
• Special-Purpose Reachability Index: Although solely
using the two-ball index is sufficiently efficient for (almost) strongly connected graphs, its performance degrades terribly when handling unreachable pairs. To
cope with this issue, as another auxiliary data structure,
we propose a special-purpose dynamic reachability index,
which maintains the reachability of fixed pairs of vertices
rather than all pairs.

2.2

Algorithms for Dynamic Networks

In recent years, several studies have proposed on exact
dynamic algorithms [20, 22, 23, 25]. These algorithms store
auxiliary data structures to avoid re-computing the betweenness centralities from scratch. Lee et al. [25] decompose agraph into subgraphs and reduce the search space for each
update. They achieve several times faster updates in comparison with Brandes’ exact algorithm. Kas et al. [22] extend the dynamic all-pair shortest path algorithm by Ramalingam et al. [31] to enable dynamic updates of betweenness
centrality. Green et al. [20] accelerates dynamic updates
by maintaining a shortest path tree (SPT) rooted at eachvertex in a graph and updating distance and the numberof shortest paths from a root. Kourtellis et al. [23] further improve the scalability by storing SPTs on disk in acompact format and distributing the computation. Among
exact dynamic algorithms, only this algorithm can handle
million-node networks, although it uses 100 machines and
16TB of space. However, all of these exact algorithms require Ω(n2 ) space to store their data structures and it isdifficult to handle large networks in memory.
The first semi-dynamic approximate algorithm was proposed by Bergamini et al. [7], which can handle each edge

We conducted extensive experiments using real large networks for demonstrating the performance of our method andfor empirically analyzing our new techniques. The results
validate the high scalability of our method, in that the proposed method can handle large-scale complex networks withover 1B edges in memory. In addition, each update is processed in several milliseconds on all the datasets that were
examined in this study. Moreover, the proposed method’s
processing time is faster than that of previous state-of-theart algorithms for dynamic networks.
Organization. The remainder of this paper is organized asfollows. We discuss related work in Section 2. In Section 3,
we explain definitions and notations. We present our hypergraph sketches, two-ball index and special-purpose reachability index in Sections 4, 5 and 6, respectively. In Section 7,
we present experimental results. We analyze the time andspace complexity of our method in Section 8. Finally in
Section 9, we present our conclusion.

49

insertion on million-node networks in less than ten milliseconds and yield further speedups by processing batches ofedge insertions. The betweenness centrality estimation methodof this algorithm is based on [32]. For each shortest path,
this algorithm stores an SPT rooted at one endpoint of it.
By maintaining the shortest path distance and the numberof shortest paths from the root of each SPT, the contribution of each shortest path is kept up to date. However, Ω(n)space is required to store the information in each SPT and
thousands of SPTs might be used to estimate betweenness
centrality. As a result, this algorithm faces difficulty in handling billion-scale real-world networks in memory.

3.

PRELIMINARIES

Notations.
In this study, a network is represented byan unweighted directed graph G = (V, E), where V is aset of n vertices and E is a set of m edges. For a vertex v ∈ V , the set of out-neighbors and in-neighbors are−→←−denoted by N (v) and N (v), respectively. Let d(u, v) bethe (shortest path) distance from a vertex u to a vertex
v. If there is no path from u to v, d(u, v) is considered as∞. For two vertices s and t, let P (s, t) be the set of vertices on at least one of the shortest paths from s to t, i.e.,
P (s, t) = {v ∈ V | d(s, v) + d(v, t) = d(s, t)}.
Betweenness Centrality. The number of shortest paths
from s to t is denoted by σ(s, t), and the number of shortest
paths from s to t passing through v is denoted by σ(s, t | v).
If s = v or t = v, let σ(s, t | v) be zero. The betweenness
centrality of a vertex v on a graph G is
X σ(s, t | v).
C(v) =σ(s, t)s,t∈V

3.1

eH (v)] = C(v) for any
Theorem 1. [36] We have EH [Cv ∈ V . Moreover for , δ > 0, by choosing M = Ω( 12 log 1δ ),eH (v) − C(v) ≥ n2 ] ≤ δ for any v ∈ V .
we have Pr[ CIn our method, we always keep the value of wH (v) withrespect to the latest graph for every v ∈ V . Hence, we cananswer the estimated betweenness centrality of any vertexin O(1) time.
We also note that complex networks have vertices thoughwhich most vertex pairs have a shortest path pass [6]. In
other words, betweenness centralities of those vertices areΩ(n2 ). Since we are interested in only those vertices in many
applications, the guarantee given by Theorem 1 is sufficientfor practical purpose.

4.2

Dynamic Graphs and Problem Definitions

4.2.1

Algorithm Description

Suppose that a graph Gτ is obtained by updating a graphGτ −1 . We want to efficiently obtain the hypergraph sketch
for Gτ using the hypergraph sketch H for Gτ −1 .
Adding or deleting an edge is trivial. We simply need
to update the distance information of hyperedges. What is
non-trivial is efficiently updating it, and we will look at thisissue in Sections 5 and 6.
When adding or deleting vertices, we sometimes need toreplace vertex pairs in the sampled set S since we need a setof pairs that are uniformly sampled from the current vertex
set Vτ . The pseudo-code of our updating procedures is given
in Algorithm 1.

HYPERGRAPH SKETCH AND ITS EXTENSION TO THE DYNAMIC SETTING

In this section, we first introduce the data structure, called
the hypergraph sketch, which is used to estimate the betweenness centrality of vertices. Then, we explain how to
update the hypergraph sketch when the graph is dynamically updated. For simplicity, we only explain a high-leveldescription of our method in this section, which is sufficient
to verify its correctness. Although we use several auxiliary
data structures to achieve scalability, the detailed implementation is deferred to Sections 5 and 6.

4.1

Dynamic Updates of Hypergraph Sketch

In this section, we describe how to dynamically update
the hypergraph sketch and verify its correctness.

When we consider dynamic networks, G denotes the latest
graph, and Gτ denotes the graph at time τ . For simplicity,
we assume that time is described by positive integers (i.e.,graph snapshots are G1 , G2 , . . .). We add subscripts to indicate the time we consider. For example, dτ (u, v) denotes
the distance between u and v in the graph Gτ .
In this paper, we study methods that, given a dynamic
graph, construct and manage a data structure to manage
the approximate value of betweenness centrality Cτ (v) for
each vertex v ∈ V , where τ denotes the latest time.

4.

chosen vertex pair (s, t) ∈ S, a hyperedge with an auxiliary information on each vertex is added to a hypergraph
H. Specifically, est := {(v, σ(s, v), σ(v, t)) | v ∈ P (s, t)} is
added to H. We call H a hypergraph sketch. In the following sections, we simply call this set a hyperedge. For a
hyperedge est , V (est ) denotes the set of vertices containedin e, that is, P (s, t).
By conducting a (bidirectional) BFS from s to t, the set
P (s, t) and the number of shortest paths σ(s, v) and σ(v, t)
for every v can be computed in O(m) time. Hence, we
can construct a hypergraph sketch with M hyperedges inO(M m) time.
We define the weight of v on a hypergraph sketch H byPP
= (s,t)∈S σ(s,v)·σ(v,t). Given awH (v) := (s,t)∈S σ(s,t|v)
σ(s,t)
σ(s,t)eH (v) = n2 wH (v) as an estimation ofvertex v, we output CM
eH (v) is an unbiased estimator of C(v)C(v). The value C
that is well concentrated:

4.2.2

Correctness

Let H1 , H2 , . . . be the distribution of hypergraph sketches
for G1 , G2 , . . ., respectively, where Gτ is obtained from Gτ −1
by applying the update method described in the previous
section. For any τ and a hypergraph sketch H for Gτ coneH (v) is an unbiased
structed by our algorithm, the value C
estimator of Cτ (v) that is well concentrated:

Hypergraph Sketch for Static Graphs

Theorem 2. For any τ , the following hold: We haveeH (v)] = Cτ (v) for any v ∈ V . Moreover for , δ >EH∼Hτ [CeH (v) −0, by choosing M = Ω( 12 log 1δ ), we have PrH∼Hτ [ C2Cτ (v) ≥ n ] ≤ δ for any v ∈ V .

We now describe the hypergraph sketch for static graphs.
Given a graph G = (V, E) and an integer M , we sample
a set S of M vertex pairs uniformly at random. For each50

Algorithm 1 Vertex operations
1: procedure AddVertex(H, v)
2:
Let Gτ be obtained from Gτ −1 by adding v.
3:for each est ∈ E(H) do
4:
continue with probability |Vτ −1 |2 /|Vτ |2 .
5:
Sample (s0 , t0 ) ∈ (Vτ × Vτ ) \ (Vτ −1 × Vτ −1 ).
6:
Replace est by the hyperedge es0 t0 made from (s0 , t0 ).

hypergraph sketch. In this section, we describe our two-ballindex (TB-index ), which addresses this issue.

5.1

7: procedure RemoveVertex(H, v)
8:
Let Gτ be obtained from Gτ −1 by deleting v.
9:for each est ∈ E(H) do
10:
if s 6= v and t 6= v then continue.
11:
Sample (s0 , t0 ) ∈ Vτ × Vτ uniformly at random.
12:
Replace est by the hyperedge es0 t0 made from (s0 , t0 ).
In words, for any τ and a vertex v ∈ V , the probability that
CH (v) is far apart from Cτ (v) can be made arbitrarily small
by choosing M sufficiently large. In particular, the accuracyof our method does not deteriorate over time.
Let Sτ be the distribution of the set S of vertex pairs attime τ , that is, H ∼ Hτ is obtained by sampling S ∼ Sτ .
Then compute a hyperedge for each pair in S. To prove
Theorem 2, by Theorem 1, it suffices to show that S is a
uniform distribution over the sets of M vertex pairs in Vτ .
Proof of Theorem 2. We use the induction on τ . When
τ = 1, the distribution Sτ is clearly uniform from the construction.
Let Gτ be the current graph, and assume that the distribution Sτ −1 is a uniform distribution over sets of M vertex
pairs in Vτ −1 . Suppose that we are to modify G. When
adding or deleting an edge, we do not modify the pairs, andhence Sτ remains a uniform distribution. When adding ordeleting a vertex, Sτ remains a uniform distribution from
Lemmas 3 and 4, as given below.

These sets have important roles in our method. We describe−→←−
B (v, dv ) and B (v, dv ) balls whose centers are v and radius
are dv . The set of vertices whose distance from v is less than−→−→or equal to dv is denoted by V ( B (v, dv )). Let B (v, dv )[w] =−→−→d(v, w) if w ∈ V ( B (v, dv )). Otherwise, B (v, dv )[w] = ∞.
←−←−We define V ( B (v, dv )) and B (v, dv )[w] in the same way as−→−→
V ( B (v, dv )) and B (v, dv )[w].
Then, our data structures are as follows.

Lemma 3. Suppose that Gτ = (Vτ , Eτ ) is obtained fromGτ −1 = (Vτ −1 , Eτ −1 ) by adding a vertex v, and Sτ −1 is a
uniform distribution over sets of M vertex pairs in Vτ −1 .
Then, the distribution Sτ is also a uniform distribution oversets of M vertex pairs in Vτ .
Proof. Since Sτ −1 is a uniform distribution, the processof sampling S ∼ Sτ can be regarded as follows. Each time
we sample a pair (s, t) ∈ Vτ −1 × Vτ −1 uniformly at random,
we keep it as is with probability |Vτ −1 |2 /|Vτ |2 , and replace itwith a pair sampled from (Vτ ×Vτ )\(Vτ −1 ×Vτ −1 ) uniformly
at random with the remaining probability. Hence, Sτ is alsoa uniform distribution.

• A set ∆s,t = {δs (v) | v ∈ P (s, t)}. Here δs (v) is supposed1 to be d(s, v). This information is not needed toestimate the betweenness centrality, but it is useful todetect the change of shortest paths from s to t when anedge or a vertex is deleted.
−→←−• Two sets β s and β t of vertices with distance informa−→tion, which are supposed to be the balls B (s, ds ) and←−B (t, dt ), respectively. Here, two non-negative integers
ds and dt satisfy ds + dt + 1 = d(s, t) − x. The valuesof ds and dt are determined when the two balls are constructed. The newly formed shortest paths from s to t
by edge insertions are detected by using these balls.

Lemma 4. Suppose that Gτ = (Vτ , Eτ ) is obtained fromGτ −1 = (Vτ −1 , Eτ −1 ) by deleting a vertex v, and Sτ −1 isa uniform distribution over sets of M vertex pairs in Vτ −1 .
Then, the distribution Sτ is also a uniform distribution oversets of M vertex pairs in Vτ .
Proof. Since Sτ −1 is a uniform distribution, the processof sampling S ∼ Sτ can be regarded as follows. Each time
we sample a pair (s, t) ∈ Vτ −1 × Vτ −1 uniformly at random,
we keep it as is if (s, t) ∈ Vτ × Vτ and replace it with apair sampled from Vτ × Vτ uniformly at random otherwise.
Hence, Sτ is also a uniform distribution.

5.

Data Structure

For each sampled vertex pair (s, t) ∈ S, the hyperedge est
must be updated when the set of shortest paths from s tot is changed by dynamic updates. Hence, quickly detecting
the change of shortest paths from s to t is important toefficiently update the approximate betweenness centralities.
A straightforward way of detecting the change of shortest
paths from s to t is conducting a (bidirectional) BFS froms to t. Obviously, traversing the whole graph is inefficientfor computing shortest paths between thousands of vertex
pairs. Another approach is to maintain the shortest pathtree (SPT) rooted at the vertex s. Several dynamic algorithms [7, 20, 23] are based on variants of this approach because the incremental updates of an SPT can be processed
quickly. However, keeping thousands of SPTs on billionscale networks requires huge amount of space, and it is difficult to handle such networks in memory.
To achieve both high performance and scalability, we construct and store data structures using balls. We define a set−→←−of pairs B (v, dv ) and B (v, dv ) as follows.
−→
B (v, dv ) = {(w, d(v, w)) | d(v, w) ≤ dv }←−B (v, dv ) = {(w, d(w, v)) | d(w, v) ≤ dv }The data structures are parameterized by a non-negative
integer x, which allows the trade-off between the index size
and the update time: As the value of x increases, the index size and the update time of edge deletions decrease
whereas the update time of edge insertions increases. Two−→←−balls B (s, ds ) and B (t, dt ) never overlap. When x = 0, twoballs are directly connected by an edge on a shortest path
between s and t, i.e., the distance between the two balls is
1. By increasing x, two balls become smaller, and the twoballs become distant.

TWO-BALL INDEXIn order to compute approximate betweenness centrality
in a dynamic setting, we need to be able to update the hyperedge est efficiently for each vertex pair (s, t) ∈ S in the1

We use the expression “supposed to be” here because,
while δs (v) usually corresponds to d(s, v), during updateprocedure, δs (v) may temporarily differ from d(s, v).

51

−
→
Algorithm 2 Update B (s, ds ) after edge (u, v) is inserted

The set of triplets
n
o
−
→ ←
−
I = (∆s,t , β s , β t ) | (s, t) ∈ S−→

is the TB-index for the set S of sampled vertex pairs.
−→ ←−The worst-case space complexity of each triplet (∆s,t , β s , β t )−→←−is O(n). However, the size of B (s, ds ) and B (t, dt ) tends
to be much smaller than n on real-world complex networks
in practice. This is because these networks have the smallworld property, i.e., the average distance between two vertices is often O(log n) [34]. Hence the radius of two balls
ds and dt tend to be quite small. Moreover, the size of thehyperedge est is further smaller than the size of two balls.
The compactness of the TB-index on real-world networks is
empirically shown in Section 7.

5.2

1: procedure InsertEdgeIntoBall(u, v, β s )2:
Q ← An empty FIFO queue.
−
→
−
→
3:
if β s [v] > β s [u] + 1 then
−
→
−
→
4:
β s [v] ← β s [u] + 1; Q.push(v).
5:while not Q.empty() do
6:v ← Q.pop().
−→7:
if β s [v] = ds then continue.
8:
for each (v, c) ∈ E do
−
→
−
→
9:
if β s [c] > β s [v] + 1 then
−
→
−
→
10:
β s [c] ← β s [v] + 1; Q.push(c).

5.3Incremental Update

In this subsection, we present how to efficiently update−→ ←−each triplet (∆s,t , β s , β t ) in the TB-index as well as the hyperedge est when the graph is dynamically updated. Throughout this section, τ denotes the latest time.

Index Construction

In this subsection, we explain how to efficiently construct−→ ←−a triplet (∆s,t , β s , β t ) for a fixed vertex pair (s, t). Then,the construction of I is simply performed by computing thetriplet for each sampled vertex pair (s, t) ∈ S. This procedure is used when we start dealing with a new dynamic
graph, as well as when sampled pairs are changed due tovertex addition and removal (see Section 4.2).
−→←−
First, two balls B (s, ds ) and B (t, dt ) are computed basedon bidirectional BFS. We initialize ds and dt to be zero.
−→←−Then, as long as B (s, ds ) and B (t, dt ) share no vertex, we−→
increment ds if the size of B (s, ds ) is smaller than that of−→←−
B (t, dt ), and increment dt otherwise. Then, we set β s =
←
−
−
→
←
−
−
→
B (s, ds ) and β t = B (t, dt ). If we find B (s, ds + 1) =
−
→
←
−
←
−
B (s, ds ) or B (t, dt +1) = B (t, dt ) along the way, then there−→−→
is no path from s to t. In this case, we simply set β s = β t =∅ to save space.
−→←−
When we find that B (s, ds ) and B (t, dt ) share a vertex,−→
the vertex set P (s, t) is computed by using B (s, ds ) and
←
−
−
→
←
−
B (t, dt ) . Let C = V ( B (s, ds )) ∩ V ( B (t, dt )) be the setof shared vertices. Let Ps and Pt be the set of vertices on
the shortest paths from s to C and the shortest paths fromC to t, respectively. Note that P (s, t) equals Ps ∪ Pt ∪ C,
because C = {v ∈ P (s, t) | d(s, v) = ds } and every shortest
path from s to t passes through C. Two vertex sets Ps andPt can be easily obtained by conducting BFSes from C to s
and t. After the vertex set P (s, t) is obtained, we decrease
the radius ds and dt one by one as long as ds + dt + 1 >d(s, t) − x holds in order to save space. More specifically,in each iteration, the radius of the larger ball is decreased.
Finally, the set ∆s,t is computed by conducting a BFS on−→←−
P (s, t) (when β s and β t are empty, we simply set ∆s,t = ∅to save the space). Since the number of vertices in P (s, t) ismuch smaller than n, the BFS on P (s, t) can be conducted
quickly. Similarly, we update the distance information in est
by conducting BFSes on P (s, t).
−→
Since the sizes of P (s, t) and the two balls B (s, ds ) and←−B (t, dt ) are often much smaller than n on real-world networks, we use hash tables to keep σ(s, v) and σ(v, t) for each
−
→
−
→
←
−
v ∈ V (est ), β s [v] for each v ∈ V ( β s ), and β t [v] for each←−−→v ∈ V ( β t ) in order to save the space. As with balls, β s [v]denotes d(s, v), which is stored in the hash table.

5.3.1

Edge Insertion−→
Suppose that an edge (u, v) is inserted. Two balls β s←−and β t are updated first; then, we update ∆s,t and thehyperedge est . The value of ds and dt are not modified
although we may have a new shortest path between s and t.
Algorithm 2 shows the pseudo-code for updating each ballafter edge insertion. Balls are updated in a similar manneras in the case of previous methods for maintaining SPTs [1,−→−→17, 33]. In the case that β s [u] + 1 < β s [v], the distancefrom s to v becomes smaller by passing through the new−→−→edge (u, v), and β s [v] is updated to β s [u] + 1. Then werecursively examine out-neighbors of v. In other words, the−→ball β s is updated by conducting a BFS from v on vertices←−
whose distance from s is at most ds . The ball β t is updated
in the same manner as above.
−→
The set ∆s,t is updated by using the updated balls β s =
←
−
−
→
←
−
B τ (s, ds ) and β t = B τ (t, dt ). The change of ∆s,t occurs
when new shortest paths are formed by adding (u, v). To
quickly detect newly formed shortest paths, two more balls←−−→
B τ (u, du ) and B τ (v, dv ) are used. The radii of two balls←−
du and dv are initialized to zero. The radius of B τ (u, du )is increased one by one by conducting a BFS from u while−→←−B τ (s, ds ) and B τ (u, du ) share no vertex, and the radius of−→−→B τ (v, dv ) is increased one by one as long as B τ (v, dv ) and←−B τ (t, dt ) share no vertex. However, if vertices u and v are
distant from vertices s and t, a large part of the graph might←−−→
be visited before the whole B τ (u, du ) and B τ (v, dv ) are obtained even when there is no newly formed shortest paths
passing through the edge (u, v). To avoid such an unnecessary computation, we carefully introduce upper bounds ru
and rv on du and dv , respectively, and stop BFSes when theradii reach these upper bounds. We consider the following
four cases:−→←−
1. The case where u ∈ V ( B τ (s, ds )) and v ∈ V ( B τ (t, dt )):In this case, we can immediately detect that a new shortest path from s to t is formed.
−→←−
2. The case where u ∈ V ( B τ (s, ds )) and v 6∈ V ( B τ (t, dt )):−→In this case, we already know that d(s, u) = B τ (s, ds )[u] =−→β s [u]. Thus, the following lemma holds.

52

−
→
Algorithm 3 Update B (s, ds ) after edge (u, v) is deleted.

Lemma 5. If there exists a shortest path from s tot passing through the edge (u, v), then there exists dv ≤−→−→dτ −1 (s, t)− B τ (s, ds )[u]−dt −1 such that V ( B τ (v, dv ))∩←−V ( B τ (t, dt )) 6= ∅.

−→1: procedure HasParent(w, β s )2:for each (p, w) ∈ E do
−
→
−
→
3:
if β s [p] + 1 = β s [w] then return true.
4:
return false.

−→
Therefore, we can set rv = dτ −1 (s, t) − β s [u] − dt − 1.
Note that the index has the information of dτ −1 (s, t).
−→←−
3. The case where u 6∈ V ( B τ (s, ds )) and v ∈ V ( B τ (t, dt )):As in the previous case, we can set ru = dτ −1 (s, t) −←−β t [v] − ds − 1.
−→←−
4. The case where u 6∈ V ( B τ (s, ds )) and v 6∈ V ( B τ (t, dt )):In this case, the following lemma holds.

−→

5: procedure DeleteEdgeFromBall(u, v, β s )
−
→
−
→
−
→
6:
if {u, v} 6⊆ V ( β s ) ∨ β s [u] + 1 6= β s [v] then return7:
Step 1: Collect vertices w with d(s, w) increased.
8:D ← ∅.
9:Q ← An empty FIFO queue.
−→10:
if not HasParent(v, B (s, ds )) then11:
D ← {v}.
−→
12:
β s [v] ← ∞; Q.push(v).
13:while not Q.empty() do
14:w ← Q.pop().
15:
for each (w, c) ∈ E do
−
→
−
→
16:
if β s [c] < ∞∧ not HasParent(c, β s ) then17:
D ← D ∪ {c}.
−→18:β s [c] ← ∞; Q.push(c).
−→19:Step 2: Update β s [w] for each w ∈ D.
0
20:
Q ← An empty min-based priority queue.
21:
for each w ∈ D do n
o
−
→
22:
dw
β s [p] | (p, w) ∈ E .
min ← minp∈V

Lemma 6. If there exists a shortest path from s to t
passing through the edge (u, v), then there exist du , dv ≤−→←−x such that V ( B τ (s, ds )) ∩ V ( B τ (u, du )) 6= ∅ and−→←−V ( B τ (v, dv )) ∩ V ( B τ (t, dt )) 6= ∅.
Proof. Assume that there are no such du and dv .
−→←−If V ( B τ (s, ds )) ∩ V ( B τ (u, du )) = ∅, then we havedτ (s, u) + dτ (u, v) + dτ (v, t) ≥ dτ (s, u) + 1 + dt
>ds + x + 1 + dt = dτ −1 (s, t)

23:
24:
25:
26:
27:
28:
29:
30:
31:
32:

Hence, there is no newly formed shortest path from s
to t, which is a contradiction. Similarly, the case where−→←−V ( B τ (v, dv )) ∩ V ( B τ (t, dt )) = ∅ reaches a contradiction.
From the lemma above, we can set ru = rv = x.
←−
If there exist du ≤ ru and dv ≤ rv such that V ( B (u, du ))∩
−
→
−
→
←
−
V ( B (s, ds )) 6= ∅ and V ( B (v, dv )) ∩ V ( B (t, dt )) 6= ∅, thenthere may exist a newly formed shortest path passing throughthe edge (u, v). In this case, the set of vertex Pτ (s, t) is contained in the vertex set P := V (est ) ∪ Pτ (s, u) ∪ Pτ (v, t).
Therefore, ∆s,t and est can be efficiently updated by conducting BFSes on P . We note that P is often small becauseof the small average distance.

5.3.2

−
→
0
β s [w] ← dw
min + 1; Q .push(w).
while not Q0 .empty() dow ← Q0 .pop().
−→d ← β s [w].
if d = ds then continue.
for each (w, c) ∈ E do−→
if β s [c] > d + 1 then−→β s [c] ← d + 1; Q0 .push(c).
for each w ∈ D do−→−→
Remove w from β s if β s [w] > ds .

Step 2: In this step, for each w ∈ U , the distance dτ (s, w)−→
is computed so as to determine β s [w]. For each vertex−→
w
w
w ∈ U ,n β s [w] is initializedo to dmin + 1, where dmin =−→minp∈V β s [p] | (p, v) ∈ E and pushed into a min-basedpriority queue Q0 .
−→Suppose that a vertex w with the minimum β s [w] is0popped from Q . The distance dτ (s, w) is determined foreach w ∈ U by conducting Dijkstra’s algorithm using thepriority queue Q0 . After the algorithm ends, for each w ∈ U ,−→−→the vertex w is deleted from β s if β s [w] is larger than ds .
Now we describe how to update the hyperedge est . We
only use the information of ∆s,t to update est . Note that est
should be updated only if the set of shortest paths from s
to v is changed, i.e., {u, v} ⊆ V (est ) ∧ ∆s,t (u) + 1 = ∆s,t (v)hold. There are two cases to consider:

Edge Deletion

Suppose that an edge (u, v) is deleted from the graph. We−→
first present a procedure for updating the ball β s . The ball←−
β t can be updated similarly.
−→−→
Ball B (s, ds ) changes only if {u, v} ⊆ V ( B τ −1 (s, ds )) and
−
→
−
→
−
→
B τ −1 (s, ds )[u]+1 = B τ −1 (s, ds )[v] hold. In this case, β s is
updated in the same way as the procedure for maintaining
SPTs [17]. Algorithm 3 shows the pseudo-code for updating
a ball after an edge deletion. This algorithm mainly consistsof two steps.
Step 1: The set U of vertices whose distances from s are
increased is obtained by conducting a BFS from v. We first
push v into a queue Q.
Let w be the vertex popped from Q. Each neighbor cof w is examined and checked whether−→ there still−→exists avertex p that satisfies (p, c) ∈ E and β s [p] + 1 = β s [c]. If
there is no such vertex, it implies that the distance d(s, w)
is increased after the edge deletion. In this case, w is added
to U and enqueued to Q, and neighbors of w are examinedrecursively.

1. dτ (s, t) 6= dτ −1 (s, t): This case happens if and only ifστ −1 (s, u) = στ −1 (v, t) = 1 holds. This condition can beeasily checked by looking at variables σ(s, u) and σ(v, t)−→ ←−in est . If this happens, the triplet (∆s,t , β s , β t ) as wellas the hyperedge est are constructed from scratch again.
2. Otherwise: In this case, although dτ (s, t) = dτ −1 (s, t)holds, the number of shortest paths may have changed.
Since the vertex set Pτ (s, t) is contained in V (est ), we

53

can obtain est after the edge deletion by conducting a
BFS from s to t on V (est ).

5.3.3

𝑠

𝑟

Vertex Operations
𝑡

𝑡

𝑡

(a)

(b)

(c)

Figure 1: An example of index size reduction. Vertex color indicates its status: Blue is a vertex thatshould be maintained for vertex pair (s, t), Red is avertex that can be reached from a vertex r, Gray
is a vertex that can reach r, and Yellow is a vertex
contained in the same SCC with r.
reachable vertices from a single vertex suffices. However,
SPTs can be updated more easily than other alternatives.
−→←−
Let r be a vertex contained in the SCC. Let T r and T r
be the SPT rooted at r on the original graph and the graphwith each edge reversed, respectively.
Tr
We define the reduced reachable vertices Rst
from s as
follows.

−
→
←
−

if s ∈ V ( T r ) ∧ t ∈ V ( T r ),
∅
−
→
−
→
Tr
Rst = Rs \ V ( T r ) if t ∈/ V ( T r ),Rsotherwise.

SPECIAL-PURPOSEREACHABILITY INDEXThe algorithm explained in the previous section works wellon social networks. However, on web graphs, the performance of our method becomes much worse than that on social networks. This is because most web graphs are directed
and there are many pairs of vertices that are unreachable
from each other. Since we only keep a triplet of three empty
sets for each unreachable vertex pair (s, t) ∈ S in the TBindex in order to save space, we have to check whether there
exists a path from s to t after each edge insertion. Checking
the reachability of thousands of vertex pairs is not acceptable for handling dynamic updates efficiently.
Therefore, to efficiently process the update after edge insertions, we need to construct an index that quickly answers
and maintains reachability of each vertex pair (s, t) ∈ S. In
this section, we present an index for this purpose. By exploiting the structure of web graphs [12], we construct ahighly scalable index. Our index outperforms in terms of
the performance and scalability compared to state-of-theart methods that are designed to answer the reachability ofany vertex pair [35, 38].

6.1

𝑠
𝑟

First, we update the vertex pair as described in Section 4.
If the vertex pair is replaced, then we construct the hyperedge as well as the triplet from scratch.
When adding a vertex, we have nothing to do. Whenremoving a vertex, the update procedure is similar to the
update procedure when an edge is deleted. We first update−→←−
two balls β s and β t . Then, we detect the set U of vertices such that the distances from the centers of balls are
increased, and then, the distances from the centers are determined using Dijkstra’s algorithm from U . We also update
the hyperedge est if v is contained in V (est ). If the distancefrom s to t is changed after the vertex removal, the triplet−→ ←−(∆s,t , β s , β t ) is constructed from scratch. Otherwise, theupdated hyperedge est is obtained by conducting BFS froms and t on V (est ).

6.

𝑠

−
→ ←
−
Tr
Using T r , T r and Rst, the reachability from s to t can beeasily answered.
Lemma 7. There exists a path from s to t if and only if
−
→
←
−
Tr
s ∈ V ( T r ) ∧ t ∈ V ( T r ) or t ∈ Rst
holds.
Figure 1 shows an example of index size reduction using
our technique. Let us consider a graph and a vertex pair(s, t) as shown in Figure 1a. In this graph, all the vertices in
the graph can be reached from a vertex s. Hence, without−→←−
pruning by using T r and T r , all the vertices in the graphshould be maintained for (s, t). Let us choose a vertex r in←−−→
the graph as the root of two SPTs T r and T r . As shown in−→
Figure 1b, the vertex t is not contained in V ( T r ), and re−→Tr
duced reachable vertices Rst
becomes Rs \V ( T r ). Thus, weTr
only have to maintain four vertices included in Rst
insteadof ten vertices that can be reached from s (Figure 1c).
Since the SCC will change over time, it is difficult for us toquickly check which vertex is contained in the SCC. Hence,
we randomly choose a set of K vertices {r1 , r2 , . . . , rK } asroots of SPTs. If we set the value K to a moderate value,at least one of the chosen vertices will be contained in the
SCC with high probability because of a large number of thevertices belong to the SCC [12]. We define the set T of pairsof SPTs as follows.
o
n−
→ ←
−
−
→ ←
−
−
→
←
−
T = ( T r1 , T r1 ), ( T r2 , T r2 ), . . . , ( T rK , T rK )Data Structure

To quickly answer the reachability of vertex pairs, we wantto manage a vertex set Rs for each vertex pair (s, t), where
Rs is the set of all vertices that can be reached from s. The
reachability from s to t can be easily answered by checking
whether t ∈ Rs or not. Obviously, keeping all the reachable
vertices spoils the scalability of our method owing to thelarge number of reachable vertices. To reduce the numberof vertices stored in the index, we exploit the structure of
web graphs. On web graphs, the size of vertex set Rs tends
to be large when the largest strongly connected componentof the graph, called the SCC [12], can be reached from s.
We introduce another data structure with which we canreduce the size of the index size for a vertex pair (s, t) suchthat there is a path from s to t passing through the SCC.
To maintain the reachability of such vertex pairs, we use
SPTs rooted at vertices in the SCC. As distance information
is not necessary here, any data structure that can maintainWe denote the set of rootsthat ←
can
from
n
−
→
− be reached ←− oavertex v by Tin (v) = ri | ( T ri , T ri ) ∈ T , v ∈ V ( T ri ) .
We denotenthe set of roots that can reachoa vertex v by
−
−
→
−
→ ←
Tout (v) = ri | ( T ri , T ri ) ∈ T , v ∈ V ( T ri ) . The following lemma holds from the definition of Tin and Tout .
Lemma 8. There exists a path from s to t passing throughat least one of the vertices in {r1 , r2 , . . . , rK } if and only if
Tin (s) ∩ Tout (t) 6= ∅ holds.

54

Tr
Now we define, by generalizing the notion of Rst, the setTof reachable vertices Rst reduced by T for each vertex pair(s, t) as follows:(∅
if Tin (s) ∩ Tout (t) 6= ∅,
T
S
−
→
Rst =
→
−
Rs \ 1≤i≤K:t∈V
V(
T)
otherwise.
ri/ (T r )iT
We note that Rst
is represented as an SPT to support increT
mental updates, i.e., not only vertices in Rst
but also theirdistances from s are stored in memory.
In the case where Tin (s) ∩ Tout (t) = ∅ holds, the reducedTreachable vertices Rst
is computed by conducting a pruned
BFS from s using the information of T . Most procedures inpruned BFS are the same as standard BFS. The difference
between them is that, when we examine a neighbor w ofa vertex v, we skip processing w if Tout (w) ∩ T out (t) 6= ∅holds, where T out (t) = {r1 , r2 , . . . , rK } \ Tout (t). We notethat Tout (w) ∩ T out (t) 6= ∅ implies there is no path from w
to t since otherwise there exists a path from ri ∈ Tout (w) ∩T out (t) to t via w, which contradicts ri ∈ T out (t).
T
By means of T and Rst
we can easily check the reachability from s to t.

7.

EXPERIMENTS

We conducted experiments on many real-world complex
networks. The index size, the performance, and the accuracy of our method were compared with other state-of-theart algorithms. In this section, we refer to the algorithm
proposed by Bergamini et al. [7] as BMS. Since exact dynamic methods [20,22,23,25] require Ω(n2 ) space and it wasdifficult to examine them on large networks, only BMS wasused as our competitor for performance comparison. Unlessotherwise mentioned, our indices are constructed with 1,000
chosen vertex pairs and ten roots of SPTs, that is, M = 1000
and K = 10.
Environment. All the algorithms were implemented with
C++11 and Google Sparse Hash and compiled with gcc 4.8.3
with -O3 option, which will be available from our websites.
All experiments were conducted with a single thread on aLinux server with Intel Xeon E5-2670 2.60GHz CPU and
512GB memory.
Datasets. All experiments were conducted on real-worldsocial networks and web graphs that are publicly available.
The information of datasets is shown in Table 1. The average distance is estimated by 10,000 random vertex pairs.
Five networks twitter-2010, in-2004, indochina-2004, it-2004
and uk-2007 were downloaded from Laboratory for Web Algorithmics [8, 9]. The other networks used in the experiments were downloaded from Stanford Network Analysis
Project [26].
Scenario. We first randomly deleted 200 edges from thegraph. Second, we constructed indices to measure the construction time and data structure size. Then, we inserted
the 200 edges into the graph one by one, and measured theaverage edge insertion time. Next, we deleted these edges
from the graph one by one and measured the average edge
deletion time. After that, we added 200 new vertices to theindex constructed from each dataset and measured the average vertex insertion time. Finally, we deleted 200 randomly
chosen vertices from the graph and measured the averagevertex deletion time. Unless otherwise mentioned, we repeated this procedure five times and the average values arereported.

Lemma 9. There exists a path from s to t if and only ifT
Tin (s) ∩ Tout (t) 6= ∅ or t ∈ Rst.
 T
Let R = Rsi ti | 1 ≤ i ≤ M . Then, the index for reachability query is the pair of two sets (T , R). The set T is
obtained by conducting BFSes 2K times, and the set R is
obtained by conducting pruned BFSes M times. Since thevalue K does not have to be large, say less than 30, the condition Tout (w) ∩ T out (t) = ∅ can be checked in O(1) time byrepresenting sets of roots Tout (w) and Tout (t) with a bit vector. Thus, the time complexity of pruned BFS is also O(m),and the overall time complexity of our reachability index isO((K + M )m). Our index consumes O((K + M )n) space,but the actual memory usage is much smaller in practice.

6.2

deletion. U2 is obtained in the same way as Step 1 of Algorithm 3. Let U = U1 ∪ U2 . For each w ∈ U , it is pushed
into a min-based priority queue in the same way as Step 2of Algorithm 3. The distance from s to each vertex w ∈ U isdetermined through Dijkstra’s algorithm by using the pruning method in pruned BFS. For each vertex v ∈ U , v isT
added to Rst
if d(s, v) < ∞.
Vertex Operations.
When adding a vertex, we havenothing to do. When removing a vertex, we first remove itfrom each SPT in T . For each sampled vertex pair (s, t) ∈S, if Tout (t) is changed after vertex removal, we computeT
Rstfrom scratch. Otherwise, we collect the vertex set U
Twhose distances from s are changed, and we update Rst
asis done when edges are removed. When some root vertex
ri is removed from a graph, we choose a new root vertex
uniformly at random.

Incremental Update

When an edge is inserted or deleted, all the SPTs in T are
updated first. Each SPT is updated in the same way as the
update of balls except that the distance from the root is notTlimited. Then, for each vertex pair (s, t), the vertex set Rst
is re-computed. If the set of roots Tout (t) is changed after theTupdate, the vertex set Rst
is also significantly changed. InTthis case, reduced reachable vertices Rst
are computed fromscratch by conducting pruned BFSes again. Otherwise, weT
update Rst
in a way described below.
Edge Insertion. Assume that an edge (u, v) is insertedT
into the graph. Some vertices will be newly added to Rst
because of the appearance of paths passing through the edgeT(u, v). These vertices is detected and added to Rst
by conducting a pruned BFS from v in the same manner as theTupdate of balls. In some cases, there exists a vertex w ∈ Rst
satisfying Tout (w) ∩ T out (t) 6= ∅ owing to the change inT
Tout (w), i.e., the vertices included in Rst
becomes a superT
set of the original definition of Rst . However, our reachability index still answers the reachability between each vertex
pair correctly.
Edge Deletion.
Assume that an edge (u, v) is deleted
from the graph. A vertex w ∈ V may be reached from s
with pruned BFS if Tout (w) ∩ T out (t) = ∅ holds after the
update. Let U1 be the set of such vertices. Let U2 be theset of vertices whose distance from s increased after the edge

7.1

Construction Time

Table 2 shows that our index can be efficiently constructed
on all datasets. Since the index construction time is notstrongly affected by the value of the trade-off parameter x,

55

social(d)social(u)social(d)social(d)social(d)social(u)social(d)social(u)NotreDame
GoogleBerkStan
in-2004
indochina-2004
it-2004
uk-2007web(d)
web(d)
web(d)
web(d)
web(d)
web(d)
web(d)

m

d

D

35K
37K
77K
1.6M
4.8M
3.1M
42M
66M

422K
368K
905K
31M
69M
117M
1.5B
1.8B

12.2
10.0
11.7
18.8
14.2
38.1
35.3
27.5

11.7
4.0
4.1
5.3
5.9
4.2
4.5
5.0

326K
876K
685K
1.4M
7.4M
41M
106M

1.5M
5.1M
7.6M
17M
194M
1.2B
3.7B

4.6
5.8
11.1
12.2
26.2
27.9
35.3

Insertion time per edge (ms)

HepPh
Enron
Slashdot0811
Pokec
LiveJournal
Orkut
twitter-2010
Friendster

n

8
7
6
5
4
3
2
1
0

BMS
EI0
EI1

10

0

10

1

Batch size

10

2

10

3

7
BMS
EI0
EI1

6
5
4
3
2
1
0

100

101

Batch size(a) Slashdot0811

11.4
11.7
13.7
15.3
15.7
15.0
15.4

10
BMS
EI0
EI1

8
6
4
2
0

100

101

Batch size

102(c) BerkStan

102

103

(b) Pokec

Insertion time per edge (ms)Type

Insertion time per edge (ms)

DatasetInsertion time per edge (ms)Table 1: Information of datasets. The average degree and average distance are denoted by d and D,
respectively.

103

6
BMS
EI0
EI1

5
4
3
2
1
0

100

101

Batch size

102

103

(d) in-2004

Figure 2: The processing time of batch updates.
the maximum average construction time among x = 0, 1, 2
is reported.
The index can be constructed in about 40 minutes forFriendster, which is the largest social networks used in our
experiments. On other networks, our index can be constructed in 15 minutes. Thus, the construction time of ourindex is highly efficient in comparison with the computation
time of exact betweenness centrality.

7.2

time of the index. It can also be seen from Table 2 thatthe edge insertion time increases as the parameter x getslarger. This bears the trade-off of the index size and theedge insertion time.
The average edge deletion time is less than one millisecondon all datasets when x = 0. Moreover, when we increase theparameter x, the edge deletion time becomes even faster.
This is because deleted edges affect balls less often when
the balls are small.
Vertex insertion/deletion times are also shown in Table 2.
Since the update times for vertex insertions and deletions are
not strongly affected by the value of the trade-off parameterx, the maximum average update time among x = 0, 1, 2 isreported.
Vertex insertion time ranges from 0.1 milliseconds to 35.0
milliseconds. Vertex deletion time ranges from 0.2 milliseconds to 7.3 milliseconds. Although these update times are
slower than edge insertion/deletion times, this performanceis acceptable because the number of vertex operations tends
to be fewer than the number of edge operations.
Batch Update. As BMS can gain update speed by handling multiple edge insertions as a batch, we also evaluatedthe insertion time when multiple edges are inserted at onceand handled as batches. In this experiment, our method stillprocesses each edge insertion one by one. Figure 2 shows theinsertion time per edge for different sizes of batches of multiple edge insertion. The results of our method with x = 2
are not shown because the processing time is slower than the
single edge update of BMS. The processing time for a single
edge update of our method with x = 0 is at least five times
faster than that of BMS. Even when the size of a batch is
1,024, on social networks, the processing time per edge ofour method is still faster than that of BMS.

Data Structure Size

The total index size of our method and the average size ofdata structures stored to support efficient updates of eachhyperedge are listed in Table 3. We also report the memory
usage for naturally representing each graph, i.e., each edge
appears twice in forward and backward adjacency lists, andrepresented by four-byte integers.
As shown in Table 3, we successfully constructed indiceson both billion-scale social networks and web graphs. Especially on uk-2007, which is the largest network used in our
experiments, the index size becomes smaller than the original graph size. In comparison with the index of BMS, ourindex is 20 times smaller on all datasets even when x = 0.
We can see that the size of data structures for each sampled vertex pair in our method is several orders of magnitude
smaller than the size of each SPT stored in BMS. This indicates the high scalability of our method with regard to
hyperedges. For example, even if we increase the number ofhyperedges M from 1,000 to 10,000, the index size of our
method on uk-2007 would be smaller than 60 GB. As shownin Table 3, we can further reduce the data size for each hyperedge by increasing the parameter x. Especially on social
networks, it becomes 10 times smaller by increasing x fromzero to two, which have a great impact when M is large.

7.3

Update Time

Average update times for insertions and deletions of edges
and vertices are also listed in Table 2. As shown in thetable, each edge insertion is processed in two milliseconds
on almost all datasets when x = 0. This processing timeis three orders of magnitude faster than full re-computation

7.4

Parameter M and Accuracy

We compared the accuracy of our method and BMS under various number of samples M . We used BMS only since56

Table 2: The timing results of our index on real-world datasets. IT denotes the index construction time. EIx
and EDx denote the average edge insertion/deletion time, respectively, where x is the trade-off parameter.
VI and VD denote the average vertex insertion/deletion time. For IT, VI and VD of our method, maximum
values among x = 0, 1, 2 are reported. DNF means it did not finish in one hour or ran out memory.
Ours

DatasetBMS [7]

IT

EI0

EI1

EI2

ED0

ED1

ED2

VI

VD

IT

EI

Slashdot0811
Pokec
LiveJournal
Orkut
twitter-2010
Friendster

626.4 ms
18.8 s27.5 s87.1 s931.1 s2211.4 s

0.9 ms
0.4 ms0.3 ms
1.1 ms219.6 ms
0.8 ms

6.8 ms
3.4 ms
3.3 ms
22.4 ms
1.0 s
17.6 ms

205.5 ms
82.3 ms
164.4 ms
3.3 s
DNF
14.2 s0.2 ms0.2 ms0.2 ms
0.5 ms0.3 ms
0.4 ms0.2 ms0.2 ms
0.1 ms
0.4 ms0.2 ms0.3 ms

0.1 ms0.2 ms
0.1 ms
0.4 ms
1.2 ms0.3 ms

0.1 ms
0.7 ms
1.5 ms
0.9 ms
14.7 ms
24.4 ms

0.3 ms0.2 ms0.2 ms
1.0 ms
7.3 ms
1.3 ms21.8 s810.1 s2491.1 s
DNF
DNF
DNF

6.1 ms
7.3 ms
8.0 ms———

NotreDame
GoogleBerkStan
in-2004
indochina-2004
it-2004
uk-2007564.3 ms
3.1 s
4.3 s
3.1 s
23.5 s135.6 s419.6 s0.3 ms
0.4 ms
0.8 ms
0.5 ms
0.7 ms
0.9 ms
1.7 ms

0.9 ms
1.4 ms
5.2 ms
3.3 ms51.7 ms
16.0 ms
31.8 ms

3.5 ms
3.4 ms
16.5 ms
35.7 ms
12.0 s346.2 ms
1.5 s

0.2 ms0.2 ms
0.4 ms0.3 ms
0.4 ms
0.5 ms
0.9 ms0.2 ms0.2 ms0.3 ms0.3 ms0.3 ms
0.4 ms
0.8 ms

0.1 ms0.2 ms0.3 ms0.2 ms0.3 ms
0.4 ms
0.7 ms

0.1 ms
0.5 ms
0.4 ms
0.5 ms
2.4 ms12.9 ms
35.0 ms

0.5 ms0.2 ms
0.4 ms0.2 ms
0.6 ms
0.5 ms
0.6 ms4.9 s102.2 s71.1 s85.8 s777.2 s
DNF
DNF

3.2 ms
8.3 ms
14.3 ms
3.9 ms
6.3 ms——Table 3: The index size on real-world datasets. ISx denote the total index size, and HSx denote the data
structure size for each vertex pair, where x is the trade-off parameter. PS denotes the size of each SPT in
BMS [7]. The memory usage for representing each graph is denoted by G.
Ours

DatasetBMS [7]

G

IS1

IS2

HS0

HS1

HS2

IS

PS

Slashdot0811
Pokec
LiveJournal
Orkut
twitter-2010
Friendster

23.0 MB394.6 MB
1.1 GB
867.3 MB9.3 GB
14.7 GB

20.1 MB368.7 MB
1.1 GB
704.2 MB9.2 GB
14.4 GB

19.9 MB362.8 MB
1.1 GB
683.3 MB9.2 GB
14.4 GB

6.0 KB
36.0 KB
55.5 KB
192.6 KB
108.4 KB
313.5 KB

3.1 KB
10.1 KB
14.2 KB
29.5 KB
22.9 KB
34.8 KB

2.9 KB
4.3 KB
5.6 KB
8.6 KB
23.5 KB
14.0 KB

542.7 MB11.4 GB
34.0 GB———

541.6 KB
11.4 MB33.9 MB———

7.2 MB
245.0 MB552.0 MB937.5 MB11.7 GB14.4 GB

NotreDame
GoogleBerkStan
in-2004
indochina-2004
it-2004
uk-200790.4 MB248.8 MB290.9 MB514.9 MB
2.2 GB
10.3 GB
26.4 GB89.3 MB230.1 MB259.8 MB483.5 MB
2.1 GB
9.8 GB
25.5 GB

88.0 MB213.1 MB229.1 MB451.2 MB
1.9 GB
9.5 GB
24.7 GB18.9 KB
56.5 KB
140.5 KB
211.2 KB
609.6 KB
1.2 MB
3.3 MB

17.8 KB
37.8 KB
109.3 KB
179.9 KB
453.3 KB
704.4 KB
2.4 MB

16.4 KB
20.8 KB
78.6 KB
147.5 KB
312.5 KB
406.8 KB
1.6 MB

2.3 GB
6.1 GB
4.8 GB
9.7 GB
52.0 GB——

2.3 MB
6.1 MB
4.8 MB
9.7 MB
51.9 MB——12.0 MB40.8 MB60.8 MB135.3 MB
1.6 GB
9.2 GB
29.9 GB

Maximum absolute error [×10−2]

other estimation methods are developed for the static setting.
For each number of samples M = 250, 500, 1, 000, 2, 000,
4, 000, 8, 000, 16, 000, we compared the approximate centralities of vertices with the exact values obtained by Brandes’ exact algorithm [10]. Due to the high cost of exact
betweenness centrality computation, the comparisons wereconducted on small networks HepPh and Enron. The maximum and average absolute errors of estimated centrality
values were evaluated, where we define the error of an estie
right aftermated centrality of a vertex v as n12 C(v) − C(v)the index construction in order to compare the accuracy on
networks of varying size. In Figure 3a and Figure 3b, theaverages of these values for ten trials are shown. As shownin Figure 3a, both methods show better accuracy by increasing the number of samples, as suggested by their theoreticalguarantees. Figure 3a and Figure 3b show that our method
estimates centralities more accurately than BMS does in almost all settings, although BMS has better theoretical guarantee. It may be because hypergraph sketch uses more information for each sampled vertex pair (s, t), i.e., the numberof shortest paths σ(s, v) and σ(v, t) for each v ∈ P (s, t), to

3.0

HepPh (Ours)
HepPh (BMS)
Enron (Ours)
Enron (BMS)

2.5
2.0
1.5
1.0
0.5
0.0

103

Number of samples

104(a) Maximum error

Average absolute error [×10−4]

IS0

1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

HepPh (Ours)
HepPh (BMS)
Enron (Ours)
Enron (BMS)103

Number of samples

104(b) Average error

Figure 3: Betweenness centrality estimation error
1eC(v) − C(v)on small networks.
n2estimate betweenness centrality.

7.5

Parameter K and Reachability Index

We now show the index size and construction time of our
special-purpose reachability index under various number ofroot vertices K. The results show the effectiveness of ourtechnique in reducing the index size and how many SPTs

57

107

Time (s)

T
Rst

∗

105
104
103
102

10

3

10

2

NotreDame
GoogleBerkStan
in-2004
indochina-2004

Vertex

101

Edge

100
10−1

101
100

Table 4: Running time of each update operation.

104

NotreDame
GoogleBerkStan
in-2004
indochina-2004

106

0

5

10

Number of SPTs

15

20

T
(a) Average size of Rst
.

10−2

0

5

10

Number of SPTs

15(b) Construction time.

RemovalO((M + K)T D )Insertion

O(M Ex + (M + K)T EI )DeletionO((M + K)T D )that kBx k and kRK k were at most 0.02m on million-node
networks in our experiments.
Since we need O(M Bx ), O(Kn), O(M RK ), and O(n)space to store balls, SPTs, the set of reduced reachable vertices, and the centrality values of vertices and other variables, respectively, the overall space requirement of our data
structure is O(M ( Bx + RK ) + Kn). Bx and RK aremuch smaller than n in practice; They were at most 0.01n on
million-node networks in our experiments. The space complexity of our data structure becomes significantly smaller
in comparison with the algorithm proposed by Bergamini etal. [7], which stores M SPTs.

are required to sufficiently reduce the index size.
TIndex Size.
The sizes of Rst
over 1,000 vertex pairs
were evaluated under various numbers of SPTs up to 20.
The averages over ten trials on directed graphs are shown in
Figure 4a. Sampling ten root vertices seems to be enoughTfor reducing the number of vertices stored in Rst
except for
T
NotreDame. By sampling ten root vertices, on average, Rst
becomes at least 3,000 times smaller on indochina-2004 andat least 100,000 times smaller on Google.
Construction Time. The construction time of the index
for maintaining the reachability of 1,000 vertex pairs werealso evaluated. The average construction time over ten trials
are shown in Figure 4b. As the number of SPTs increases,
the construction time at first quickly drops due to the deT
creasing size of Rst, but it gradually increases owing to theincreasing construction time of SPTs.

8.2

Updating TimeIn our data structure, we store 2M SPTs in the TB-index,
2K SPTs in the reachability index, and M SPTs for reduced
reachable vertices. Updating these SPTs is the bottleneckof our update procedures. Hence, we first consider the timecomplexity of updating SPTs.
Let C be the set of vertices whose distances from the rootof an SPT are modified by the update operation, kCk be
the number of edges incident to C, and D be the average
number of vertices whose
P distance from v ∈ C is less thanor equal to two, i.e.,v∈C {w | d(v, w) ≤ 2} / C . When
adding a vertex, we only have to allocate a new space forthe newly added vertex. When inserting an edge, we only
need to look at neighbors of vertices whose distances fromthe root decreased. When removing a vertex or an edge,
we need O(D C ) time to compute the vertex set C, andO(kCk log kCk) time to update the distance from the root
to each vertex in C. Hence, the time complexity of vertex addition, edge insertion, and vertex/edge deletion are
O(1), T EI := O(kCk), and T D := O(D C + kCk log kCk),
respectively. We remark that C is much smaller than V inpractice; The average of C on million-node networks was
at most 40 when removing a vertex from an SPT on thewhole graph and at most four in other cases.
Now we show the running time of each update operationin Table 4. Except vertex addition, the bottleneck of the
update procedure is updating SPTs, and the time complexity of each operation is O((M + K)T D ). Precisely speaking,we may replace several vertex pairs when processing vertex
operations, but we ignore its time complexity in our analysis because the probability of occurring such a replacement
is only O( n1 ). When inserting an edge, we have to check
whether new shortest paths are formed. Due to the small
size of balls, the probability that a randomly inserted edge
has an endpoint in a ball is quite small, and we only consider the case that the inserted edge has no endpoint in aball. Thus, when an edge is inserted, O(M Ex ) additionaltime is required to traverse a graph from this edge, where
Ex is the number of vertices whose distances from endpointsof the inserted edge are less than or equals to x.

COMPLEXITY ANALYSIS

In this section, we compute the construction time andspace requirement of our data structure. We also roughly
analyze the running time of each update operation. For
the sake of simplicity, we assume that a hash table can beaccessed in a constant time and the size of a hash table islinear in the number of items stored in it. Although theanalysis below lacks some rigidity for simplicity, it would
explain why our method is efficient on real-world datasets.

8.1

O(M + K)

20

Figure 4: The statistical information of the index for
maintaining the reachability of 1,000 vertex pairs.

8.

Addition

Construction Time and Space Requirement

Before starting analysis, we introduce several notations.
For an integer x, let Bx be the average size of two balls
P
−
→
←
−
for a vertex pair, i.e., (s,t)∈S ( B (s, ds ) + B (t, dt ) )/ S
under the parameter x that is given when the index is constructed, and kBx k be the average number of edges incident−→←−
to a vertex in B (s, ds ) or B (t, dt ). Let RK be the averageTsize of reduced reachable vertices Rst
and kRK k be the avT
erage number of edges incident to Rst
when the size of the
set T of SPT pairs is K.
When constructing a TB-index, the bottleneck is performing M bidirectional BFSes to compute the distance between
each vertex pair and obtaining two balls, and its running
time is O(M kB0 k). When constructing our special-purpose
reachability index, we need to perform 2K BFSes to obtain the set of K pairs of SPTs and M pruned BFSes tocompute the reduced reachable vertices for each vertex pair,
and its running time is O(Km+M kRK k). Thus, the overall
construction time is O(M (kB0 k + kRK k) + Km). We note

58

9.

CONCLUSION[15] D. Erdös, V. Ishakian, A. Bestavros, and E. Terzi. A
divide-and-conquer algorithm for betweenness centrality. In
SDM, 2015. to appear.
[16] L. C. Freeman. A Set of Measures of Centrality Based on
Betweenness. Sociometry, 40(1):35–41, 1977.
[17] D. Frigioni, A. Marchetti-Spaccamela, and U. Nanni. Fully
Dynamic Algorithms for Maintaining Shortest Paths Trees.
J. Algorithms, 34(2):251–281, 2000.
[18] R. Geisberger, P. Sanders, and D. Schultes. Better
Approximation of Betweenness Centrality. In ALENEX,
pages 90–100, 2008.
[19] M. Girvan and M. E. J. Newman. Community structure insocial and biological networks. P. Natl. Acad. Sci.,
99(12):7821–6, 2002.
[20] O. Green, R. McColl, and D. A. Bader. A fast algorithm
for streaming betweenness centrality. In
SocialCom/PASSAT, pages 11–20, 2012.
[21] H. Jeong, S. P. Mason, A.-L. Barabási, and Z. N. Oltvai.
Lethality and centrality in protein networks. Nature,411(6833):41–42, 2001.
[22] M. Kas, M. Wachs, K. M. Carley, and L. R. Carley.
Incremental Algorithm for Updating Betweenness
Centrality in Dynamically Growing Networks. In
ASONAM, pages 33–40, 2013.
[23] N. Kourtellis, G. D. F. Morales, and F. Bonchi. Scalable
Online Betweenness Centrality in Evolving Graphs.
Transactions on Knowledge and Data Engineering, 2014.
to appear.
[24] V. E. Krebs. Mapping networks of terrorist cells.
Connections, 24(3):43–52, 2002.
[25] M.-J. Lee, J. Lee, J. Y. Park, R. H. Choi, and C.-W.
Chung. QUBE: a Quick algorithm for Updating
BEtweenness centrality. In WWW, pages 351–360, 2012.
[26] J. Leskovec and A. Krevl. SNAP Datasets: Stanford large
network dataset collection.
[27] S. Milgram. The Small World Problem. Psychology today,
2(1):60–67, 1967.
[28] M. E. J. Newman. Networks: An Introduction. Oxford
University Press, 2010.
[29] M. E. J. Newman and M. Girvan. Finding and evaluating
community structure in networks. Phys. Rev. E, 69:026113,
2004.
[30] R. Puzis, P. Zilberman, Y. Elovici, S. Dolev, andU. Brandes. Heuristics for Speeding Up Betweenness
Centrality Computation. In SocialCom/PASSAT, pages
302–311, 2012.
[31] G. Ramalingam and T. W. Reps. On the Computational
Complexity of Incremental Algorithms. University of
Wisconsin-Madison, Computer Sciences Department, 1991.
[32] M. Riondato and E. M. Kornaropoulos. Fast
Approximation of Betweenness Centrality ThroughSampling. In WSDM, pages 413–422, 2014.
[33] K. Tretyakov, A. Armas-cervantes, L. Garcı́a-ba nuelos,
and M. Dumas. Fast Fully Dynamic Landmark-based
Estimation of Shortest Path Distances in Very Large
Graphs. In CIKM, pages 1785–1794, 2011.
[34] D. J. Watts and S. H. Strogatz. Collective dynamics of’small-world’ networks. Nature, 393(6684):440–442, 1998.
[35] H. Yildirim, V. Chaoji, and M. J. Zaki. DAGGER: A
Scalable Index for Reachability Queries in Large Dynamic
Graphs. CoRR, abs/1301.0:11, 2013.
[36] Y. Yoshida. Almost Linear-Time Algorithms for Adaptive
Betweenness Centrality using Hypergraph Sketches. In
KDD, pages 1416–1425, 2014.
[37] R. Zafarani, M. A. Abbasi, and H. Liu. Social Media
Mining: An Introduction. Cambridge University Press,
2014.
[38] A. D. Zhu, W. Lin, S. Wang, and X. Xiao. Reachability
Queries on Large Dynamic Graphs: A Total OrderApproach. In SIGMOD, pages 1323–1334, 2014.

In this paper, we present the first approximate betweenness centrality algorithm for fully dynamic networks. Our
TB-index efficiently detects the change of shortest paths between vertex pairs and updates approximate betweenness
centralities of vertices after each update. The size of the TBindex is maintained small enough to handle large scale complex networks. We also proposed a fully dynamic specialpurpose reachability index for maintaining the reachabilityof some fixed vertex pairs in order to further improve theperformance of the TB-index on directed graphs. The set of
SPTs rooted at vertices in the core substantially reduces thenumber of vertices maintained in our reachability index. Experimental results showed the TB-index is compact enoughto handle dynamic updates on very large complex networks
with tens of thousands of vertices and over one billion edges
in memory. Moreover, each update after edge insertion andedge deletion can be processed in several milliseconds.

Acknowledgments. T. H., T. A., and Y. Y. are supported
by JST, ERATO, Kawarabayashi Large Graph Project. T. A.
is supported by JSPS Grant-in-Aid for Research Activity Startup (No. 15H06828) and PRESTO, JST. Y. Y. is supported by
JSPS Grant-in-Aid for Young Scientists (B) (No. 26730009) andMEXT Grant-in-Aid for Scientific Research on Innovative Areas(No. 24106003).

10.

REFERENCES[1] T. Akiba, Y. Iwata, and Y. Yoshida. Dynamic andhistorical shortest-path distance queries on large evolving
networks by pruned landmark labeling. In WWW, pages
237–248, 2014.
[2] J. M. Anthonisse. The Rush In A Directed Graph.
Technical Report BN 9/71, Stichting MathematischCentrum, 1971.
[3] D. A. Bader, S. Kintali, K. Madduri, and M. Mihail.
Approximating Betweenness Centrality. In WAW, pages
124–137, 2007.
[4] M. Baglioni, F. Geraci, M. Pellegrini, and E. Lastres. Fast
Exact Computation of Betweenness Centrality in Social
Networks. In ASONAM, pages 450–456, 2012.
[5] A.-L. Barabasi. The origin of bursts and heavy tails inhuman dynamics. Nature, 435:207–211, 2005.
[6] M. Barthélemy. Betweenness centrality in large complex
networks. The European Physical Journal B - Condensed
Matter and Complex Systems, 38(2):163–168, 2004.
[7] E. Bergamini, H. Meyerhenke, and C. L. Staudt.
Approximating Betweenness Centrality in Large Evolving
Networks. In ALENEX, pages 133–146, 2015.
[8] P. Boldi, M. Rosa, M. Santini, and S. Vigna. Layered Label
Propagation: A MultiResolution Coordinate-Free Ordering
for Compressing Social Networks. In WWW, 2011.
[9] P. Boldi and S. Vigna. The WebGraph Framework I:
Compression Techniques. In WWW, pages 595–601, 2004.
[10] U. Brandes. A Faster Algorithm for Betweenness
Centrality. J. Math. Sociol., 25(2):163–177, 2001.
[11] U. Brandes and C. Pich. Centrality Estimation in Large
Networks. Int. J. Bifurcat. Chaos, 17(07):2303–2318, 2007.
[12] A. Broder, R. Kumar, F. Maghoul, P. Raghavan,
S. Rajagopalan, R. Stata, A. Tomkins, and J. Wiener.
Graph Structure in the Web. Comput. Netw.,
33(1-6):309–320, 2000.
[13] T. Coffman, S. Greenblatt, and S. Marcus. Graph-based
technologies for intelligence analysis. Commun. ACM,
47(3):45–47, 2004.
[14] A. Del Sol, H. Fujihashi, and P. O’Meara. Topology ofsmall-world networks of protein–protein complex
structures. Bioinformatics, 21(8):1311–1315, 2005.

59