| {z }
=:C
ρ . (3.1)
Matrix C will play a crucial role in computing centralities.
3.3 Current-Flow Measures of Centrality
Two of the most widely used centrality measures are based on a model of non-splitting information transmission along shortest paths. Note that in the following, distances may well be defined in terms of an edge length (or resistance) r: V −→R+.
Definition 3.3.1 ([3, 52]) (Shortest-path) betweenness cB is defined by cB: V −→R≥0
v 7−→cB(v) := 1 nB
X
s,t∈V
ns,v,t ns,t ,
where ns,t denotes the number of shortest paths from s to t, ns,v,t denotes the number of shortest paths from s to t with v as an inner vertex, and nB := (n−1)(n−2)is a normalizing constant (nB:=n(n−1) ifv may also be a start or end vertex).
It thus measures the degree to which a vertex is participating in the commu-nication between pairs of other vertices.
Definition 3.3.2 ([5]) (Shortest-path) closeness centrality cC is defined by cC: V −→R+
v 7−→cC(v) := nC P
t6=v
dist(v, t) ,
where dist(v, w) denotes the length of a shortest path between v and w and nC :=n−1 is a normalizing constant.
It thus measures the degree to which a vertex is close to other vertices (on average).
Both measures assume that information (or whatever else is being mod-eled) flows along shortest paths, and does not split. We next describe two alternative measures that build on the same intuition, but let information flow and split like current in an electrical network.
3.3.1 Current-Flow Betweenness Centrality
In electrical networks, the analog of the fraction of shortest st-paths passing through a vertex (or an edge) is the fraction of a unit st-current flowing through that vertex (or edge). Given a supply ρ, we therefore define the throughput of a vertex v ∈V to be
τ(v) := 1
2 −|ρ(v)|+X
e∈v
|I(eσ)|
! ,
where the term −|ρ(v)| accounts for the fact that only inner vertices are considered in the definition of shortest-path betweenness centrality. To in-clude start and end vertex, it should be replaced by +|ρ(v)|. Accordingly, the throughput of an edge e∈E is defined as
τ(e) =|I(eσ)|.
Let τs,t denote the throughput in case of an st-current.
Definition 3.3.3 ([98]) LetG= (V, E, c)be an electrical network. Current-flow betweenness centrality cCB: V −→R≥0 is defined by
cCB(v) := 1 nB
X
s,t∈V
τs,t(v) .
For the following reason, current-flow betweenness is also calledrandom-walk betweenness. A simple random st-walk is a random walk that starts at s, ends in t, and continues at vertex v 6= t by picking an incident edge e ∈ E with probability c(e)/P
e03vc(e0). Then, given an st-current, the amount of current flowing through a particular edge eσ equals the expected difference of the number of times that the simple random st-walk passes edge eσ along and against its orientation (see, e.g., [6]).
3.3.2 Current-Flow Closeness Centrality
Similar to the above variation of betweenness centrality, we utilize the ana-log of shortest-path distance in electrical networks to introduce a variant of closeness centrality.
Definition 3.3.4 Let G = (V, E, c) be an electrical network. Current-flow closeness centrality cCC: V −→R+ is defined by
cCC(s) := nC P
t6=s
Us,t(s)−Us,t(t) for all s∈V .
If we want to assure that cCC is bounded into ]0,1], then we have to redefine nC := 2(n−1)/n. Current-flow closeness centrality is well-defined, because any two absolute potentials differ only by an additive constant. Since we only consider unit st-currents, the term Us,t(s)−Us,t(t) corresponds to the effective resistance Reff(s, t), which can be seen as an alternative measure of distance between s and t.
In [109], Information centrality is introduced by means of the statistical design of experiments. Information is defined as the inverse of the vari-ance of some given stochastic signals. One can show that the varivari-ance then corresponds to effective resistance, and information corresponds to current.
Information centrality cI: V −→R+ is defined by
cI(s)−1 :=nCs,sI + tr(CI)−2/n , (3.2) where CI = (L+J)−1 with Laplacian matrix L and J :=11T. Information centrality is often referred to, but not frequently used; most likely because its underlying intuition is not widely understood.
Theorem 3.3.5 Current-flow closeness equals information centrality.
Proof: We first note that (3.2) can be rewritten as cI(s)−1 =X
t∈V
Cs,sI +Ct,tI −Cs,tI −Ct,sI . (3.3) On the other hand, current-flow closeness can be rewritten into the same form, though with matrix C introduced in (3.1),
cCC(s)−1 =X
t6=s
Us,t(s)−Us,t(t) = X
t6=s
Cs,s−Cs,t−(Ct,s−Ct,t)
=X
t∈V
Cs,s+Ct,t−Cs,t−Ct,s .
We show that these terms are actually equal using a matrix C, such thatb CI =C+C, that contributes zero to the common summation scheme.b
For j = 1, . . . , n let ej := (0, . . . ,0,1,0, . . . ,0)T, with 1 at position j, and let C·,j and C·,jI denote the column j of the corresponding matrix. The columns of C are uniquely determined by
LC·,j =ej −e1 and C1,j = 0 and those of CI satisfy
(L+J)C·,jI =ej (3.4)
by definition. Projecting (3.4) onto the kernel of L, i.e., multiplying both sides with 11T/n from the left, yields
J C·,jI = (1TC·,jI )1 = 1 n1 .
Equation (3.4) is therefore equivalent to
It is easily verified thatCbcontributes zero when subjected to the summation
scheme of (3.3).
3.4 Improved Exact Computation
For comparison note that shortest-path betweenness and closeness can be computed in O(nm+n2logn) time and O(n+m) space using an efficient implementation of Dijkstra’s algorithm, see [10].
For current-flow betweenness centrality, matrix C defined in (3.1) is de-termined by inverting the restricted Laplacian L. Sincee Us,t = Cρs,t and I((v, w)) = (U(w)−U(v))·c({v, w}), we can use the (weighted) incidence
to compute st-currents −Is,t = BTCρs,t. From the entries of current-flow matrix F :=BTC the centrality scores are then determined via
cCB(v) = 1
where vj < vk if j < k. The total time for current-flow betweenness is thus inO(M(n−1) +mn2), see [98], whereM(n)∈ O(n3) is the time required to compute the inverse of an n×n-matrix. Note that M(n) ∈ Ω(n2logn) for arbitrary real matrices.
This can be improved as follows (see Algorithm2).
Algorithm 2: Current-flow betweenness
Input: electrical networkG= (V, E, c) with vertices v1, . . . , vn Output: current-flow betweenness cCB: V −→R≥0
begin cCB ←0 C←
0 0T 0Le−1
for e∈E do
Fe,· ← BTC
e,·
sort rowFe,· in non-increasing order
2.1
pos(e, v)← rank of Fe,v in sorted rowFe,·
for j = 1, . . . , ndo
increase cCB(source(eσ)) by (j−pos(e, vj))·Fe,vj
increase cCB(target(eσ)) by (n+ 1−j−pos(e, vj))·Fe,vj for j = 1, . . . , ndo
cCB(vj)←(cCB(vj)−j + 1)·2/nB 2.2
end
Theorem 3.4.1 The computation time of current-flow betweenness is in O(M(n−1) +mnlogn).
Proof: We refer to Algorithm 2. We can compute cCB(v) by summing up only the inflows, i.e., positive current on an edge directed to v or negative current on an edge leaving v, as follows. Note that for every non-outlet the inflow is equal to the outflow by KCL. We will later take care of the outlets.
The total inflow τin intov equals
Inflows include the vanishing unit current whenever v is the sink. In the summation over all pairs s < t this will be the case j −1 times, namely when v =vj. Note that the inflow of the source is always zero. Subtracting the vanishing currents from the total inflow yields half of the current-flow betweenness. The relation
cCB(vj) = 2(τin(vj)−j+ 1) , is accounted for in Line 2 of the algorithm.
The computational bottleneck after determiningC by matrix inversion is the sorting of rows in Line 2, which takes O(mnlogn) time. Note that F is computed by multiplying C with an incidence matrix, so that it takes only
time in O(mn).
Information centrality can be computed by determining matrix CI defined in the previous section and evaluating (3.2). The total running time is thus O(M(n) +n).
Using the new interpretation as current-flow closeness centrality, we see that it can also be determined from C rather than CI (see Algorithm 3).
Thus sparseness is preserved and only one matrix inversion is required to compute both closeness and betweenness, which corresponds nicely to the fact that shortest-path betweenness and closeness can be computed during the same traversals.
A straight-forward approach for matrix inversion uses Gaussian elimina-tion, leading to a computation time of O(n3). For networks from many ap-plication areas, however, sparse matrix techniques are appropriate. Since Le is symmetric and positive definite, the conjugate gradient method (CGM) can be used with an incomplete Cholesky decomposition as a preconditioner.
This yields a running time of O(mn√
κ), where κ is the condition number of Le times its approximate inverse obtained by applying the preconditioner, see Sect.A.4. A rule of thumb for the condition number isκ∈Θ(n), see [64],
Algorithm 3: Current-flow closeness Input: electrical networkG= (V, E, c)
Output: current-flow closeness cCC: V −→R+
begin cCC ←0 for v ∈V do
C·,v←
0 0T 0Le−1
3.1 ·,v
for w∈V do
increase cCC(v) by Cv,v−2Cw,v increase cCC(w) by Cv,v
for v ∈V do
cCC(v)←1/cCC(v) end
leading to a running time of O(mn1.5), which is faster than the subsequent summation before its improvement to O(mnlogn).
The inverse ofLe can be computed column-by-column as needed in Line3 of the algorithm for closeness centrality. Its memory requirement is inO(m).
For betweenness centrality, O(n2) memory is required in the worst case.
Here it is the current-flow matrix F that is processed row-by-row, implying that columns of Le−1 corresponding to vertices u and w with {u, w} ∈E are needed simultaneously. Therefore, the columnv ∈V needs to be determined only when the first row Fe,· withv ∈e is encountered, and it can be dropped from memory when the last such row has been processed.
To reduce the memory requirements of Algorithm2, we therefore seek an ordering that minimizes the maximum number of columns that have to be kept in memory at the same time. That is, we would like to determine a one-to-one mapping π: V −→ {1, . . . , n}, where
δ(π) := max
1≤j≤n|{u∈V : ∃ w∈V,{u, w} ∈E with π(u)≤j < π(w)}| ≤n is minimum. Unfortunately, this is an N P-hard problem known as vertex separation, see [87], or, equivalently, see [78],minimum pathwidth.
Heuristically, we can find a good ordering π∗ by using algorithms for bandwidth- and envelope-reduction of matrices, since the bandwidth (of the Laplacian matrix ofGordered byπ∗) is an upper bound forδ(π). Algorithm2 is easily modified to use any precomputed ordering. The proven reverse Cuthill-McKee heuristic, see [30], does not increase the asymptotic running
time, while it reduces the memory requirement to O(δ(π∗)n). Note that it can also be employed in the inversion of L.e
3.5 Probabilistic Approximation
In large networks, both running time and space requirements of the algo-rithm for current-flow betweenness are prohibitive. Note that shortest-path closeness can be approximated quickly, see [118].
We show that a similar approach can be used to reduce not only the running time, but also the space requirements of (approximate) current-flow betweenness computations. For large data sets, this is often even more im-portant.
The basic idea is that the betweenness of a vertex, i.e., the throughput over all st-currents, can be approximated using a small fraction of all pairs s 6=t ∈V. A fully polynomial randomized approximation scheme is given in Algorithm 4.
Theorem 3.5.1 There is a randomized algorithm that needs space in O(m) and time in O((1/ε2)m√
κlogn) to approximate current-flow betweenness within an absolute error of ε with high probability.
Proof: LetXv(1), . . . , Xv(k)be independent random variables that returnτs,t(v),
i.e., the scaled expected throughput of k st-currents is equal to the current-flow betweenness. Since 0 ≤τs,t(v)≤1 Hoeffding’s bound, see [71], yields
P
For each selected pair s 6= t ∈ V, the restricted system in Line 4 of Algorithm4 can be solved inO(m√
κ) time and O(m) space using CGM.
Algorithm 4: Randomized approximation scheme for current-flow be-tweenness
Input: electrical networkG= (V, E, c), threshold ε >0, constant ` Output: current-flow betweenness approximation c0CB: V −→R≥0
begin
c0CB ←0 and k←`· d(c∗/ε)2logne for j=1,. . . k do
selects 6=t ∈V uniformly at random and solve LeUe =ρes,t for v ∈V \ {s, t} do
4.1
for e={v, w} ∈E do increase c0CB(v) by c(e)· |Ue(w)−Ue(v)| ·c∗/2k
end
3.6 Experimental Results
We provide empirical evidence that the proposed algorithms for (approxi-mate) computation of current-flow betweenness is practical. It has been im-plemented (like most of the software of this thesis) in Java using the yFiles1 graph data structure and the JMP2 linear algebra package. All experiments were performed on a regular PC with 2.4 GHz clock-speed and 3 GB main memory. Constant ` = 1 for the approximation. See Fig. 3.1.
3.7 Conclusion
Current-flow betweenness and closeness are variants of (shortest-path) be-tweenness and closeness centrality for an alternative model of information spreading. In particular, we introduced current-flow closeness and proved that it is equal to information centrality, the original definition of which is rather unintuitive.
There is one and only one path between each pair of vertices in a tree, and the length of this path equals its resistance. We thus have the following result.
Theorem 3.7.1 The two shortest-path and current-flow measures agree on trees.
1www.yworks.de
2www.math.uib.no/˜bjornoh/jmp/index2.html
0 100 200 300 400 500 600 700
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
computation time in seconds
number of vertices
previous algorithm (conjugate gradient method, without sorting) improved algorithm matrix inversion (conjugate gradient method)
0 0.02 0.04 0.06 0.08 0.1 0.12
0.64 0.32 0.16 0.08 0.04 0.02
maximum error
precision threshold
Figure 3.1: Comparison of total running time for current-flow betweenness on random graphs with average degree 6,8, . . . ,20 and maximum error of ap-proximation on 6 random and 6 AS graphs with approximately 9000 vertices and 20000 edges each
Corollary 3.7.2 Betweenness and closeness can be computed in O(n) time and space on trees.
Proof: A bottom-up followed by a top-down traversal similar to [106].
Finally, we want to remark that there is a straightforward extension of shortest-path betweenness to edges (simply replace the numerators by the number of shortest st-paths that use the edge), see [3]. A similar extension of current-flow betweenness, that can be computed by slight modification of Algorithm 2, is given by
cCB(e) = 1 nB
X
s6=t∈V
τs,t(e) for all e ∈E .
Dynamic Graph Drawing
Spectral methods are naturally suited for dynamic graph layout because, usually, moderate changes of a graph yield moderate changes of the layout.
We discuss some general principles for dynamic graph layout and derive a dynamic spectral layout approach for the animation of small-world models.
4.1 Introduction
The main problem in dynamic graph layout is the balance of layout quality and mental-map preservation [92]. Typically, the problem is addressed by adapting a static layout method, such that it produces similar layouts for successive graphs. While these adaptations are typically ad-hoc [40], oth-ers [18, 11] are based on the formally derived method [19] of integrating difference metrics [22] into the static method. See [20] for an overview of the dynamic graph drawing problem. Spectral layout denotes the use of eigen-vectors of graph-related matrices, such as adjacency or Laplacian matrix as positions, see, e.g., [81]. We argue that spectral methods are naturally suited for dynamic graph layout both from a theoretical and practical point of view, because, usually, moderate changes in the graph yield moderate changes of the layout – this is in particular given for small worlds, see Sect. 4.5, where the initial ring structure provides stability, see Fig. 4.11. Furthermore, up-dates of spectral layouts can be computed efficiently.
This chapter is organized as follows. In Sect.4.2 we define basic notation and recall principles of spectral graph layout. The dynamic graph layout problem is reviewed briefly in Sect. 4.3, and methods for updates between layouts of consecutive graphs are treated in more detail in Sect. 4.4. In Sect. 4.5, our approach for small worlds is introduced, and we conclude with a brief discussion in Sect. 4.6.
4.2 Preliminaries
For ease of exposition we only consider two-dimensional straight-line repre-sentations of undirected, finite, simple graphs G= (V, E) with positive edge weightsω: E →R+, although most techniques and results in this paper eas-ily carry over to other classes of graphs. In straight-line representations, a two-dimensionallayout is determined by a functionp(v) := (x(v), y(v)) onV for the position of v. Most of the time we will reason about one-dimensional layouts p(v) = x(v).
For any graph-related matrix M(G), a spectral layout of G is defined by eigenvectors of M(G). We will only consider layouts derived from the (weighted) Laplacian matrix L(G) of G. Based on the Laplacian, a spectral layout can be defined asp= (u2, u3), whereu2andu3are unit eigenvectors to second and third smallest eigenvalues λ2, λ3 of the corresponding Laplacian matrixL. This has already been used for graph drawing in 1970 by Hall [67].
For sparse graphs of moderate size, a practical method to determine the corresponding eigenvectors is power iteration, also known as von Mises method. Starting from some initial vectorx the iteration is given by
x← Lx
|Lx| (power iteration).
It converges to a unit eigenvector un for the largest eigenvalue λn. Since we are not interested in un we use ˇL := g ·I −L, for some upper bound g ≥ λn (e.g., g := 2 degmax, see Sect. 2.1), which has the same eigenvectors with eigenvalues g = g−λ1 ≥ g −λ2 ≥ . . .≥ g −λn in reversed order. To obtain u2 and u3, respectively, x is orthogonalized with u1 = 1/√
n (and in the case of u3 also with u2) after each iteration step, i.e., the mean value P
v∈V x(v)/n is subtracted from every element of x to ensure x⊥1. Spectral layouts of larger graphs can be computed efficiently using multiscale methods, see [82].
4.3 Dynamic Layout
In our setting, a dynamic graph is a sequence G(0), G(1), . . . of graphs with small edit distance, i.e., G(t+1) is obtained from G(t), t ∈ N, by adding, changing, and deleting only a few vertices and edges.
There are two main scenarios for the animation of a dynamic graph, depending on whether the individual graphs are presented to the layout al-gorithm one at a time, or the entire sequence is known in advance. Layout approaches for theoffline scenario (e.g., [36]) are frequently based on a layout
of the union of all graphs in the sequence. A variant are 2.5D representations in which all graphs are shown at once (e.g., [43]). In the online scenario, the typical approach is to consider only the previous layout (e.g., [40]). A variant in which provisions for likely future changes are made is presented in [32].
Generally, offline methods deliver better results than online methods.
Anyhow, here we concentrate on the online scenario. Moreover, we even con-centrate only on the update step between two consecutive layouts as follows.
4.4 Updates
Assume we are given a sequence of layouts p0, p1. . . for a dynamic graph G(0), G(1), . . .. The step from pt to pt+1 is called logical update, whereas the actual animation of the transition is referred to as the physical update.
While simple, say,linear interpolation of two layouts is most frequently used in graph editors, more sophisticated techniques for morphing are available (see, e.g., [55,56,44]). General morphing strategies do not take into account the method by which origin and target layout are generated. For dynamic spectral layout, at least two additional strategies are reasonable.
4.4.1 Iteration
From now on just consider the step from t = 0 to 1. If the target layout p1
is a spectral layout, the power iteration for its own computation can and should be initialized with p0, that will usually be close to the target layout.
The power iteration then produces intermediate layouts which can be used for the physical update. A way to enhance the smoothness of morphing is needed because of the observation that the first steps of the iteration yield greater movement of the vertices when compared to later steps. Let Lˇ := g·I−L(G(1)). An iteration step then consists of computing the new layout ˇLx/|Lx|ˇ from a given layout x⊥1. Let g = λ1 ≥ λ2 ≥ . . . ≥ λn, and 1/√
n = u1, u2, . . . , un be the eigenvalues and unit eigenvectors of ˇL, respectively. If λ1 > λ2 > λ3 (otherwise proper eigenvectors and eigenvalues have to be chosen in the following) and x = Pn
`=2α`u`, α2 6= 0 the power iteration converges to u2 and
One way to handle this non-linear decay is to use layouts after appropriately spaced numbers of steps, or to use layouts only if the difference to the last used layout exceeds some given threshold in some metric. Both ways will enhance the smoothness of morphing by avoiding the drawing of many small movements at the end of the iteration process. Finally, the overall number of iterations can always be used to choose between layout quality and mental-map preservation, where a small number of iterations favors the latter case.
4.4.2 Interpolation
If both origin and target layout p0 and p1 are spectral layouts, intermediate layouts can also be obtained by computing eigenvectors of some intermediate matrices from L(G(0)) to L(G(1)). We interpolate linearly by
(1−t)L(G(0)) +tL(G(1)), 0≤t≤1 .
Layouts are computed for a sequence of breakpoints (each corresponds to a frame of the animation) 0 ≤ t1 < t2 < · · · < t` ≤ 1 (tj+1−tj constant, or proportional to sin(πj/`), depending on what kind of morphing seems to be appropriate, the latter one slowing down at the beginning and end). For every breakpoint tj+1 the iteration is initialized with the layout oftj, which allows fast convergence and small movements between two succeeding breakpoints.
Deletion and insertion of vertices have to be handled in a different manner, since the matrix dimension changes. See Sect. 4.5 for details.
Figures 4.6and 4.8 show smooth animations of this method. Theoretical justification for smoothness comes along with a theorem by Rellich [104]
applied to the finite dimensional case. Matrix (1−t)L(G(0)) +tL(G(1)) is a perturbed self-adjoint operator L(t) := L(G(0)) +t(L(G(1))−L(G(0))) with corresponding eigenvalues λj(t) and eigenvectorsuj(t), that are holomorphic with respect to t, i.e.
L(t)uj(t) = λj(t)uj(t) , (4.1) where uj(0) are eigenvectors at time t = 0 and uj(1) can be permuted by a permutation π, such that uπ(j)(1) are corresponding eigenvectors at time t = 1. Note that two eigenvectors may only have to be exchanged if their corresponding eigenvalues intersect during the time from 0 to 1. And even then, the power iteration exchanges these eigenvectors sufficiently smoothly for pleasing animations because the corresponding eigenvalues remain within the same range for some time.
Consider λ2 and u2 of a small world with n vertices, where starting from
Consider λ2 and u2 of a small world with n vertices, where starting from