---
abstract: |
    Matrices associated with graphs, such as the Laplacian, lead to numerous interesting graph problems expressed as linear systems. One field where Laplacian linear systems play a role is network analysis, e.g. for certain centrality measures that indicate if a node (or an edge) is important in the network. One such centrality measure is current-flow closeness.

    To allow network analysis workflows to profit from a fast Laplacian solver, we provide an implementation of the LAMG multigrid solver in the NetworKit package, facilitating the computation of current-flow closeness values or related quantities. Our main contribution consists of two algorithms that accelerate the current-flow computation for one node or a reasonably small node subset significantly. One algorithm is an unbiased estimator using sampling, the other one is based on the Johnson-Lindenstrauss transform. Our inexact algorithms lead to very accurate results in practice. Thanks to them one is now able to compute an estimation of current-flow closeness of one node on networks with tens of millions of nodes and edges within seconds or a few minutes. From a network analytical point of view, our experiments indicate that current-flow closeness can discriminate among different nodes significantly better than traditional shortest-path closeness and is also considerably more resistant to noise – we thus show that two known drawbacks of shortest-path closeness are alleviated by the current-flow variant.\
    **Keywords:** Laplacian linear system, current-flow closeness centrality, algebraic multigrid, commute time, network analysis
author:
- 'Elisabetta Bergamini[^1]'
- Michael Wegner
- 'Dimitar Lukarski[^2]'
- Henning Meyerhenke
bibliography:
- 'references.bib'
title: |
    Estimating Current-Flow Closeness Centrality\
    with a Multigrid Laplacian Solver
---

Introduction
============

Laplacian linear systems $Lp = b$ play a central role in various tasks in algorithmic graph theory and network analysis. Particularly the connection between the graph Laplacian $L$ and electrical networks allows analytical and algorithmic insights [@mavroforakis2015spanning; @DBLP:conf/stacs/BrandesF05]. Based on this connection, several centrality measures, which aim at identifying the most important nodes or edges in a network, have been introduced in the literature. For example, spanning edge centrality [@mavroforakis2015spanning] measures the importance of a node for the connectedness of the graph. Current-flow betweenness [@DBLP:conf/stacs/BrandesF05] indicates the participation of a node in paths between other nodes and current-flow closeness [@DBLP:conf/stacs/BrandesF05] measures the average distance from a node to the other nodes of the network. Differently from shortest-path closeness, here distance is a quantity which takes *all* paths between two nodes into account. In addition to centrality measures, Laplacian linear systems can be used for several other tasks in algorithmic graph theory and algorithmic network analysis, including sparsification [@DBLP:journals/siamcomp/SpielmanS11], graph partitioning [@DBLP:journals/pc/MeyerhenkeMS09], approximate maximum network flow [@DBLP:conf/stoc/ChristianoKMST11] and graph drawing [@DBLP:journals/tvcg/GansnerHN13].

Although algorithms for solving Laplacian linear systems quickly in practice have been proposed, such as the multigrid-based solvers LAMG [@livne2012lean] and CMG [@koutis2011combinatorial], an implementation of these solvers is not available in popular network-analysis frameworks, such as NetworkX [@Hagberg2008], igraph [@igraph], graph-tool [@peixoto_graph-tool_2014] and NetworKit [@Staudt2014]. In this paper, we implement LAMG in NetworKit, bridging the gap between performance-focused frameworks for network analysis and state-of-the-art Laplacian solvers. Also, we formulate the current-flow closeness centrality presented in [@DBLP:conf/stacs/BrandesF05] in terms of linear systems between pairs of nodes and propose two approximation algorithms which work very well in scenarios where the centrality of a single node or a subset of nodes has to be computed, according to our experimental study. To the best of our knowledge, only the NetworkX [@Hagberg2008] package contains an implementation of current-flow closeness. The approach implemented in NetworkX requires to invert the Laplacian of the graph, which takes $\Omega(n^2)$ time, since the inverse of the Laplacian is in general a dense matrix. This does not allow to scale to large networks with millions of nodes and edges, even when we want to compute the closeness of a single node or of a subset of nodes. With this approach, in fact, computing the closeness of one node is just as expensive as computing it for all nodes. On the contrary, our approach can estimate the closeness of a single node very quickly: For example, it requires less than 2 minutes on a network with 50 millions edges.

Thanks to our approach, we can study for the first time the properties of current-flow closeness in large networks with tens of millions of nodes. We compare current-flow closeness with traditional shortest-path closeness and show that the former succeeds in differentiating nodes significantly better than the latter, and is also more resilient to noise. In addition, we study the correlation between centrality measures and degrees in real-world networks, in relation to a recent theoretical result for random geometric graphs [@DBLP:journals/jmlr/LuxburgRH14]. Our experiments show that there is a strong correlation between degrees and current-flow closeness in complex networks, whereas there is basically no correlation is street networks.

Preliminaries {#sec:prelim}
=============

#### Graph basics.

Throughout the paper we consider connected undirected graphs $G = (V, E, w)$ having $n = |V|$ nodes and $m = |E|$ edges. The function $w : E \rightarrow \mathbb{R}$ assigns a weight $w(e)$ to each edge $e = \lbrace u, v \rbrace$ where $u$ and $v$ are the nodes connected by the edge. We introduce the shorter notation $\omega_{uv}$ for the weight $\omega(e)$ of an edge $e = \{u,v\}$. Disconnected graphs can be handled by treating each connected component separately. We denote by $\vec{E}$ the set of bidirected edges of $G$, i.e. $\vec{E} := \{(u,v) : u, v \in V, \{u,v\} \in E \}$.

We say $\pi = (v_1, ..., v_k)$ is a *path* in $G$ if $v_i \in V$ for $i =1,...,k$ and $\{v_i, v_{i+1}\} \in E$ for $i =1, ..., k-1$. We denote the quantity $\sum_{i=1}^{k-1} \omega_{v_i, v_{i+1}}$ as the *length* of path $\pi$. The path(s) of minimum length among all paths between two nodes $u$ and $v$ is (are) called the *shortest path(s)* between $u$ and $v$.

A Laplacian matrix $L$ is defined for an undirected graph $G$ as $L := D-A$, where $A=A(G)$ is the (weighted) adjacency matrix of $G$ and $D=D(G)$ the diagonal matrix storing the (weighted) node degrees: $D_{ii} = \sum_{j = 1}^{n} \omega_{ij}$.

#### Graphs as electrical networks. {#par:elec-nets}

One can regard a graph as an *electrical network* where each edge $\{u,v\}$ corresponds to a resistor with conductance $\omega_{uv}$ (the edge weight) or resistance $1 / \omega_{uv}$. We can interpret the conductance as the ease with which an electrical current can flow through the edge. We can associate a *supply* $b : V \rightarrow \mathbb{R}$ with the electrical network, representing the nodes where current enters and leaves the network. A positive supply $b(v)$ means that current is entering the network from node $v$ and a negative supply means that current is leaving the network. In the following, we will always assume that $\sum_{v \in V} b(v) = 0$ and that $b(s) = +1$ and $b(t) = -1$ for two nodes $s$ and $t$, and that $b(w) = 0\ \forall w \neq s,t$. In the following, we will refer to such a supply as vector $b_{st} \in \mathbb{R}^{n \times 1}$. We could interpret this as $s$ and $t$ being the two poles of a battery: this generates a current $e_{st}: {\vec{E}} \rightarrow \mathbb{R}$ flowing through the network. To each node $v$ we can associate a *potential* $p_{st}(v)$ such that the vector $p_{st} \in \mathbb{R}^{n \times 1}$ satisfies the following linear system: $$\label{eq:system}
Lp_{st} = b_{st}$$ Then, the current flowing through edge $(u,v)$ is defined as $(p_{st}(u) - p_{st}(v))/ \omega_{st}$. Notice that, since $G$ is connected, the rank of the Laplacian is $n-1$ and there are infinitely many vectors $p_{st}$ satisfying Eq. (\[eq:system\]), each of them differing from the other by an additive constant. However, the current is well defined, since it depends on the difference between two potentials.

Related Work {#sec:rel_work}
============

Solving Laplacian linear systems.
---------------------------------

We focus our description on iterative solvers due to their better time complexity on sparse graphs compared to direct solvers. Most iterative solvers reduce the norm of the residual $r = {\left\lVertb - Ax\right\rVert}$ iteratively by altering the current preliminary solution vector $x$ in every iteration. One usually stops when the (relative) residual is below a certain tolerance $\tau$, which yields a vector $x'$ that is a good enough approximation to the actual solution $x$. While there are recent advances in theory to solve special linear systems including Laplacians in nearly-linear time [@Spielman2004; @Kelner2013], those algorithms are not competitive in practice yet [@hoske2015; @boman2015evaluating]. In fact, the Conjugate Gradient (CG) algorithm outperforms the nearly-linear time algorithms in practice even though its asymptotic running time is typically higher.

A popular class of iterative algorithms to solve linear equations quickly in practice is called Algebraic Multigrid (AMG) [@ruge1987algebraic]. The basic idea is to solve the actual linear system by iteratively solving coarser (i.e. smaller) yet similar systems and projecting the solutions of those back to the original system. AMG algorithms can be distinguished by the class of matrices they can handle and the way they construct the coarser systems. Two fast algorithms that are specifically designed for solving Laplacian systems are CMG by Koutis et al. [@koutis2011combinatorial] and LAMG by Livne and Brandt [@livne2012lean].

We decided to use LAMG as linear solver due to its particular design for complex networks. To this end, LAMG alternates between two stages called *elimination* and *aggregation* to construct the coarser systems. The former eliminates low degree nodes in the corresponding graph, the latter partitions nodes into *aggregates* based on a special affinity measure [@livne2012lean]. Both stages reduce the number of nodes and thus define a coarsening mechanism. Based on an extensive evaluation, Livne and Brandt state that running times of LAMG and CMG are comparable but LAMG tends to be more robust in the sense that CMG has large outliers on a small set of systems [@livne2012lean].

Laplacian linear systems for network analysis. {#subsec:LaplaciansInNetworkAnalysis}
----------------------------------------------

The connection between the graph Laplacian and electrical networks (see Section \[par:elec-nets\]) has allowed for the solution of several graph algorithmic problems in terms of Laplacian systems. One of them is a centrality measure called *spanning edge centrality*, which indicates whether an edge is vital for the connectedness of a network [@mavroforakis2015spanning]. The notion of importance for connectedness is also helpful for graph sparsification. A sparsification algorithm takes a dense graph and wants to find a sparser representation (= with fewer edges) with the same vertex set and similar properties [@DBLP:journals/siamcomp/SpielmanS11], e. g. approximately the same cut sizes or eigenvalues. Since processes described by Laplacian linear systems can distinguish sparse from dense graph regions, edges in dense areas are, intuitively speaking, redundant and can be “sparsified” without doing much harm to the cut sizes when the weights of retained edges are properly scaled. This idiosyncrasy allows the use of processes described by Laplacian linear systems also for graph partitioning [@DBLP:journals/pc/MeyerhenkeMS09], approximate maximum network flow [@DBLP:conf/stoc/ChristianoKMST11], and graph drawing [@DBLP:journals/tvcg/GansnerHN13]. Moreover, the connection to electrical flow makes the use in dynamic load balancing of divisible tokens by diffusion [@DiekmannFrommerMonien99efficient] possible.

The interpretation of a graph as an electrical network has also led to the definition of two centrality measures based on current flow, current-flow closeness and current-flow betweenness [@DBLP:conf/stacs/BrandesF05]. Compared to traditional closeness and betweenness centrality, these two measures take *all* paths between two nodes into account and not only *shortest* paths.

Current-flow closeness centrality {#sec:cf_closeness}
=================================

Closeness centrality measures the efficiency of a node in spreading information to the other nodes of the network. Formally, let $d_{\text{SP}}(u,v)$ the shortest-path distance between $u$ and $v$ (i.e. the length of the shortest path(s) between $u$ and $v$). Then, closeness of node $v$ is defined as the inverse of the expected shortest-path distance between $v$ and and a random node $w$: $$\label{eq:closeness}
c_{\text{SP}}(v) := \frac{n-1}{\sum_{w \neq v} d_{\text{SP}}(v,w)}.$$ The smaller the average distance between $v$ and the other nodes, the higher is the closeness of $v$. To better understand the meaning of closeness, let us consider the two graphs in Figure \[fig:comparison\]. Since closeness takes only shortest-path distances into account, the closeness of node $x_1$ in the graph on the left and the score of node $x_2$ in the graph on the right will be exactly the same. However, there is only one path connecting $x_1$ to each of the other nodes. This means that if just a single edge is removed from the graph, $x_1$ will become disconnected from part of the other nodes. For example, let us assume the edges represent streets and $x_1$ is the location of an ambulance. If a congestion occurs, the ambulance in $x_1$ will not be able to reach part of the nodes (or it will take a very long time to reach them). On the other hand, if the ambulance was in $x_2$, a congestion would limit only partially (or not at all) the ability of the ambulance to reach the nodes of the network.

The example above illustrates that traditional closeness is unable to model scenarios where the distance between two nodes does not only depend on the length of the shortest path between them, but also on the number of shortest or relatively short paths between the nodes. For this reason, a variant of closeness named *current-flow closeness centrality* has been introduced [@DBLP:conf/stacs/BrandesF05]. If we see the graph as an electrical network, the *effective resistance* $d_{\text{ER}}(u,v) := p_{uv}(u) - p_{uv}(v)$ can be interpreted as an alternative distance measure between nodes $u$ and $v$. Indeed, if we multiply $d_{\text{ER}}(u,v)$ by the volume of $G$ (i.e. the sum of the weights of the edges in $G$), we get the *commute time* between $u$ and $v$. The commute time between nodes $u$ and $v$ is defined as $H(u,v) + H(v,u)$, where the hitting time $H(x,y)$ is the expected time step in which a random walk in the graph starting in $x$ reaches $y$ for the first time. Thus, the commute time can be seen as the expected time a random walk needs for going from $u$ to $v$ and back again. Since the commute time is based on random walks, it depends on *all* the paths between two nodes. Thus, we can use the effective resistance to define a modified centrality measure [@DBLP:conf/stacs/BrandesF05]: $$\label{eq:closeness2}
c_{\text{ER}}(v) := \frac{n-1}{\displaystyle{\sum_{w \neq v} d_{\text{ER}}(v,w)}} = \frac{n-1}{\displaystyle{\sum_{w \neq v} p_{vw}(v) - p_{vw}(w)}}.$$ By convention, we define $d_{\text{ER}}(v,v) := 0\ \forall v \in V$. To compute $c_{\text{ER}}(v)$ of a node $v$, we can solve $n-1$ linear systems. Alternatively, we could invert the Laplacian matrix $L$ of $G$ (after omitting the row and column corresponding to a node, in order to get a regular matrix $\tilde{L}$), using the property that $p_{vw}(v) - p_{vw}(w) = \tilde{L}^{-1}_{vv} - 2 \tilde{L}^{-1}_{vw} + \tilde{L}^{-1}_{ww}$ [@DBLP:conf/stacs/BrandesF05].

![Shortest-path closeness centrality cannot distinguish between node $x_1$ and node $x_2$.[]{data-label="fig:comparison"}](Comparison.pdf){width="40.00000%"}

Approximating current-flow closeness {#sec:approx}
====================================

As outlined in [@mavroforakis2015spanning], the fastest Laplacian linear solvers with a theoretical time complexity guarantee run in $O(m \log{n} \log{(1/\tau)})$ time. Multigrid methods such as CMG and LAMG are much faster in practice and have an *empirical* running time of $O(m \log{(1/\tau)})$.

Computing current-flow closeness for only one node to the desired tolerance would already require the solution of $n-1$ linear systems, yielding $O(n^2 \log{(1/\tau)})$ time in practice assuming a sparse graph. This is infeasible for large networks with millions of nodes and edges. For this reason, we propose two approximation techniques for computing current-flow closeness for a subset of the nodes in large graphs, and we compare them in our experimental evaluation. The first one is based on a simple sampling approach, which recalls the one used for classical closeness [@DBLP:journals/jgaa/EppsteinW04]. The second one uses the Johnson-Lindenstrauss transform (JLT), which allows to project the system into a lower-dimensional space by using ${\mathcal{O} \bigl (\log n \bigr )}$ random vectors.

Sampling-based approximation. {#sec:sampling}
-----------------------------

The idea is to sample uniformly at random a set $S \subseteq V$ of nodes $S = \{s_1, ..., s_k\}$, which we call *pivots*. To approximate the current-flow closeness of a node $v$, we compute the effective resistance $d_{\text{ER}}(s, v)$ between all nodes $s \in S$ and $v$. Then, the closeness of $v$ can be approximated as $$\tilde{c}_{\text{ER}}(v) := \frac{k}{n} \cdot \frac{n-1}{\sum_{i = 1}^{k} d_{\text{ER}}(v, s_i)}.$$

$\tilde{c}_{\text{ER}}(v)$ is un unbiased estimator for ${c}_{\text{ER}}(v)$ (i.e. $E[\tilde{c}_{\text{ER}}(v)] = {c}_{\text{ER}}(v)$).

We show that $Y \coloneqq \frac{n}{k} \sum_{s_i\in S} d_{\text{ER}}(v, s_i)$ is an unbiased estimator for $s(v) = \sum_{w\in V} d_{\text{ER}}(v,w)$, then the theorem follows directly from the properties of expected value. In the following, we denote the set of $k$-combinations of $V$ with $V_{k}$. $$\begin{split}
E(Y) & = \sum_{S = \{s_1, ..., s_k\} \in V_{k}} \frac{1}{\binom{n}{k}} \frac{n}{k} \sum_{s_i \in S} d_{\text{ER}}(v,s_i) = \\
    & = \frac{1}{\binom{n}{k}} \frac{n}{k} \sum_{w \in V} \binom{n-1}{k-1} d_{\text{ER}}(v,w) \\
    & = \sum_{w\in V} d_{\text{ER}}(v,w).\qquad \qquad \qquad \qquad \qquad \hfill \square
\end{split}$$

With $k$ pivots, the empirical complexity of our approach is ${\mathcal{O} \bigl (k m \log (1/\tau) \bigr )}$ with a multigrid solver. Our experiments in Section \[sec:approx\_comp\] show that a very small $k$ (e.g., $k=10$) is already enough to get a very good approximation.

Projection-based approximation. {#sec:jlt}
-------------------------------

Spielman and Srivastava [@DBLP:journals/siamcomp/SpielmanS11] show how to compute an approximation of effective resistance based on the JLT. Let $B$ be the $m \times n$ incidence matrix where each row corresponds to an edge of $G$ and each node corresponds to a node such that, for edge $e = \{u,v\}$, $B(e,u) = +1$, $B(e,v) = -1$ and $B(e, w) = 0\  \forall w \neq u,v$ (since $G$ is undirected, the direction of edge $e$ can be chosen arbitrarily). Then, they show that the effective resistance between node $u$ and node $v$ can be re-written as $d_{\text{ER}}(u,v) = ||W^{1/2} B L^{\dagger} (e_u - e_v)||_2^2$, where $W$ is the diagonal $m \times m$ matrix such that $W(e,e) = \omega(e)$, $L^{\dagger}$ is the Moore-Penrose pseudoinverse [@DBLP:books/daglib/0086372] of $L$ and $e_u$ is the $n \times 1$ vector such that $e(u) = 1$ and equal to 0 everywhere else. The effective resistances can therefore be seen as pairwise distances between vectors in $\{W^{1/2} B L^{\dagger} e_u \}_{u\in V}$, which allows to apply the JLT: If we project the vectors into a lower-dimensional space spanned by $k = {\mathcal{O} \bigl (\log n \bigr )}$ random vectors, the pairwise distances are approximately preserved. In other words, we can consider the pairwise distances between vectors in $\{QW^{1/2} B L^{\dagger} e_u \}_{u\in V}$, where $Q$ is a random projection matrix of size $k \times m$ with elements in $\{0, +\frac{1}{\sqrt{k}}, -\frac{1}{\sqrt{k}}\}$.

Since we do not want to compute $QW^{1/2} B L^{\dagger}$ directly (it would require to invert $L$), we approximate it by solving $k$ linear systems: for $i = 1,..., k$, the $i$-th row $z_i^T$ of $QW^{1/2} B L^{\dagger}$ can be computed by solving the system $L z_i = \{QW^{1/2} B\}_{\cdot, i}$, see Algorithm \[algo:JLT\] (which we reuse from [@DBLP:journals/siamcomp/SpielmanS11]). Note that the multiplication in Line \[line:multiplication\] requires only ${\mathcal{O} \bigl (2m \log n \bigr )}$ operations, since $B$ is sparse (with $2m$ non-zero entries) and $W$ is diagonal.

Construct random matrix $Q$ Compute $Y = QW^{1/2} B$ \[line:multiplication\] $Z \leftarrow$ empty $k \times n$ matrix

When choosing $k$ in Algorithm \[algo:JLT\] equal to ${\mathcal{O} \bigl (\log n/ \epsilon^2 \bigr )}$ for any $\epsilon > 0$, it was shown [@DBLP:journals/siamcomp/SpielmanS11] that, with probability $\geq 1-1/n$, $$(1-\epsilon){d}_{\text{ER}}(v,w) \leq \tilde{d}_{\text{ER}}(v,w) \leq (1+\epsilon){d}_{\text{ER}}(v,w)$$ for all $(v,w)\in V\times V$.

An approximation of current-flow closeness for node $v$ can therefore be computed as $\tilde{c}_{\text{ER}}(v) :=  (n-1)/\sum_{w \neq v} \tilde{d}_{\text{ER}}(v, w)$. If $(1-\epsilon){d}_{\text{ER}}(v,w) \leq \tilde{d}_{\text{ER}}(v,w) \leq (1+\epsilon){d}_{\text{ER}}(v,w)$ for each $w \neq v$, then also $(1-\epsilon){c}_{\text{ER}}(v) \leq \tilde{c}_{\text{ER}}(v) \leq (1+\epsilon){c}_{\text{ER}}(v)$. This is provably true only with probability $(1-1/n)^{n-1}$. However, our experimental results show that the approximation works well in practice: on all tested instances, $\tilde{c}_{\text{ER}}$ is *always* within a $(1+\epsilon)$-factor from ${c}_{\text{ER}}$ (see Section \[sec:approx\_comp\]).

Experimental Evaluation {#sec:experiments}
=======================

In this section we evaluate the performance of the two approximation algorithms described in Section \[sec:approx\]. First, we want to give some more details on the implementation, the benchmarking setup and the graph instances we used, before elaborating on the results of our evaluation.

Implementation.
---------------

We implemented both approximation algorithms in NetworKit [@Staudt2014], an open-source tool for fast exploratory analysis of massive networks. As became apparent in Section \[subsec:LaplaciansInNetworkAnalysis\], linear systems play a quite important role in network analysis and since both approximation approaches introduced in this paper rely on solving Laplacian systems, we provide the Laplacian solver LAMG by Livne and Brandt [@livne2012lean] in NetworKit with our own new implementation. The original implementation by Livne and Brandt is written in Matlab with some performance-critical parts in C and therefore difficult for us to integrate into large-scale network analysis workflows.

In informal experiments, our C++ implementation of LAMG outperforms their Matlab/C implementation regarding the solve times by a factor of 1.5 on average. In comparison the CMG solver by Koutis et al. [@koutis2011combinatorial] is on average 11% faster than our LAMG implementation on our large test instances. When solving linear systems, in all our experiments we set the relative residual error $\tau$ to $10^{-5}$.

Benchmarking Setup.
-------------------

All experiments were done on a machine equipped with 256 GB RAM and a 2.7 GHz Intel Xeon CPU E5-2680 having 2 sockets with 8 cores each and hyperthreading enabled. The machine runs 64 bit SUSE Linux and we compiled our code with g++-4.8.1 and OpenMP 3.1.

Instances. {#sub:instances}
----------

Tables \[table:small\_instances\] and \[table:instances\] in the appendix show the set of instances we use for our experiments. While Table \[table:small\_instances\] includes rather small complex networks with up to about 150000 edges, Table \[table:instances\] includes larger networks with up to 56 million edges. If a network has more than one connected component, we used the largest connected component (LCC). We ignore self-loops and the direction of edges in case a graph is directed. All the graphs are unweighted.

Approximation algorithms. {#sec:approx_comp}
-------------------------

In this section we compare the two approximation algorithms described in Section \[sec:sampling\] and Section \[sec:jlt\], respectively. We refer to the first one as <span style="font-variant:small-caps;">Sampling</span> and to the second one as <span style="font-variant:small-caps;">Projection</span>. <span style="font-variant:small-caps;">Sampling</span> depends on the number $|S|$ of samples, whereas <span style="font-variant:small-caps;">Projection</span> depends on the dimension $k$ of the $k \times n$ random projection matrix. For simplicity, we call *exact* the approach computing $c_{\text{ER}}$ as in equation \[eq:closeness2\], solving $n-1$ linear systems to the desired tolerance $\tau$.

For our experiments, we select 100 nodes for each of the networks shown in Table \[table:small\_instances\] in the appendix. For each of these nodes, we compute current-flow exactly (to the desired tolerance $\tau$) and the two approximations with different parameters. In particular, we set the number $|S|$ of samples of <span style="font-variant:small-caps;">Sampling</span> to 10, 20, 50, 100, 200, 500, and 1000. When running <span style="font-variant:small-caps;">Projection</span>, we fix $k$ to $\lceil \log n / \epsilon ^2\rceil$ and set $\epsilon$ equal to 0.5, 0.2, 0.1, and 0.05. To measure the accuracy of the algorithms, we use the well-known Spearmann rank correlation coefficient, which measures how close the ranking of nodes determined by the approximation algorithm is close to that of the exact algorithm. We recall that the closer the Spearmann coefficient is to 1, the more correlated are the two rankings, with 0 meaning no correlation and 1 meaning the two ranks are identical.

Figure \[fig:approximations\] reports the accuracy (Spearmann coefficient) and the running times in seconds for each approximation algorithm and for each parameter. We do not report explicitly to which parameter each point in the plot corresponds to, but this can be easily deduced from the running times: a smaller sample size corresponds to a smaller running time for <span style="font-variant:small-caps;">Sampling</span> and a larger $\epsilon$ corresponds to a smaller running time for <span style="font-variant:small-caps;">Projection</span>. Figure \[fig:approximations\] reports, for each approximation algorithm and for each parameter, the average over all networks of Table \[table:small\_instances\] of time and Spearmann coefficient.

The results are quite self-explanatory: the <span style="font-variant:small-caps;">Sampling</span> approach clearly outperforms <span style="font-variant:small-caps;">Projection</span> and its accuracy is extremely high already with only 10 samples. We also compute for each algorithm and parameter the number of rank inversions, i.e. the number of node pairs $\{u,v\}$ for which the approximated closeness of $u$ is smaller than the approximated closeness of $v$, but the exact closeness of $u$ is larger than or equal to the exact one of $v$ (or vice versa). With ten pivots, the average number of rank inversions of <span style="font-variant:small-caps;">Sampling</span> is 12.5; it is always below 10 for higher number of samples. This means that, out of $\binom{100}{2} = 4950$ pairs, less than 10 are inverted, corresponding to $0.2\%$.

![Time vs. Spearmann coefficient for the two approximation algorithms, using different parameters. The points represent the average among the networks of Table \[table:small\_instances\] in the Appendix.[]{data-label="fig:approximations"}](comparison_approximations.pdf){width="50.00000%"}

In addition to accuracy in terms of ranks, we also evaluate the maximum relative error. We define the relative error for a node $v$ as $e(v) = \max\{r(v), 1/r(v)\}$, where $r(v)$ is the ratio between the exact current-flow closeness of $v$ and its approximation. The maximum relative error is then defined as $\max_{v \in V}{e(v)}$. Figure \[fig:rel\_errors\] reports the results. It is interesting to notice that, with respect to this measure, the two algorithms behave quite similarly. Also, notice that the maximum relative error for <span style="font-variant:small-caps;">Projection</span> is always smaller than $\epsilon$ (we recall the values of $\epsilon$ used are 0.05, 0.1, 0.2 and 0.5), although we can only prove that this is true with probability at least $(1-1/n)^{n-1}$.

To summarize, our results show that both algorithms lead to very good accuracy in terms of maximum relative error, whereas the sampling approach better preserves the ranking of nodes, even when the number of samples is very small. For this reason, in our experiments on large graphs, we make use of the sampling approach.

![Time vs. maximum relative error for the two approximation algorithms, using different parameters. The points represent the average among the networks of Table \[table:small\_instances\] in the Appendix.[]{data-label="fig:rel_errors"}](comparison_approximations_rel_error.pdf){width="50.00000%"}

On average (over the instances of Table \[table:small\_instances\]), computing $c_\text{ER}$ on 100 nodes takes more than 20 minutes, whereas using <span style="font-variant:small-caps;">Sampling</span> with 20 pivots takes only 2.87 seconds. Table \[table:add\_results\] in the Appedix shows the detailed running times.

Comparison with shortest-path closeness.
----------------------------------------

As explained in Section \[sec:cf\_closeness\], our intuition is that current-flow closeness should represent the efficiency of a node reaching the other nodes of the network better than shortest-path closeness. To verify this assumption, we first compare the two measures in terms of their capability to discriminate between different nodes. In this experiments, we use the networks of Table \[table:small\_instances\] and compute (exactly) current-flow and shortest-path closeness on 100 randomly chosen nodes. Figure \[fig:rel\_std\_dev\] shows the relative standard deviation for shortest-path and current-flow closeness. The relative standard deviation is defined as the standard deviation divided by the average. It is always significantly higher for current-flow closeness than it is for shortest-path closeness, meaning that there is much more variation in the scores computed by the former.

![Relative standard deviation for shortest-path and current-flow closeness.[]{data-label="fig:rel_std_dev"}](coeff_var.pdf){width="50.00000%"}

Also, similarly to what has been done in [@mavroforakis2015spanning] for spanning edge centrality and edge betweenness centrality, we measure the resilience to noise, in this case for current-flow closeness and shortest-path closeness. The idea is to add edges to the graph and see how well the initial rankings are preserved. Our intuition is that, if we add some edge that creates a shortcut between a node $v$ and some other nodes, the shortest-path closeness of $v$ will be more affected than its current-flow closeness, since the former takes only shortest paths into account. This is confirmed by our experiments, summarized in Figure \[fig:resiliance\]. For each network in Table \[table:small\_instances\], we insert a percentage of the total number of edges varying from $1\%$ to $10\%$. To have a high number of shortcuts involving the sampled nodes, we always add edges between one of the sampled nodes and other nodes of the graph. Figure \[fig:resiliance\] shows, for each percentage of inserted edges, the average among all tested networks of the Spearmann correlation coefficient between the initial ranking and the ranking after the insertions. Figure \[fig:resiliance\] shows that current-flow closeness is more resilient to edge insertions and the difference between the resilience of the two measures increases the more the graph changes.

![Resilience to noise for different percentages of inserted edges. The points represent the average among the networks of Table \[table:small\_instances\] in the Appendix.[]{data-label="fig:resiliance"}](resilience.pdf){width="50.00000%"}

Correlation with degree. {#ref:corr_degrees}
------------------------

In certain random geometric graph models, such as $\epsilon$-graphs, kNN graphs, and Gaussian similarity graphs, it was recently shown [@DBLP:journals/jmlr/LuxburgRH14] that the effective resistance between two nodes $u$ and $v$ converges to $1/\deg{(u)}+1/\deg{(v)}$ when the number of nodes goes to infinity. This result also has implications on current-flow closeness: when the number of nodes goes to infinity in such graphs, $c_{\text{ER}}(v)$ goes to $c_{\text{A}}(v): = (n-1)/\sum_{w\neq v}(1/\deg{(v)} + 1/\deg{(w)} )$. However, the structure of real-world networks (e.g. complex or street networks) is significantly different from random graphs and it is not clear how close current-flow closeness is to this asymptotic value in reality. In this section we therefore study the correlation between current-flow closeness, as well as shortest-path closeness and betweenness, with $c_{\text{A}}$. In our experiments, we consider large networks with up to 56 millions nodes and edges and we use the sampling approach to approximate current-flow closeness (with 20 pivots). Table \[table:times\] shows the running times of our approximation when computing the closeness of a single node. We approximate betweenness using the approach presented in [@DBLP:conf/alenex/GeisbergerSS08], which is already implemented in NetworKit. Since the results are very different for street and complex networks, we study them separately.

Figure \[fig:corr\_degrees\] shows the Spearmann coefficient computed between the three centrality measures and $c_{\text{A}}$ on the complex networks of Table \[table:instances\] in the Appendix. The results show that there is in fact a strong positive correlation. This is weaker for shortest-path closeness, with an average Spearmann coefficient of 0.63 and stronger for betweenness and current-flow closeness, with an average of 0.81 and 0.89, respectively. We obtain similar results on the smaller instances of Table \[table:small\_instances\], where the averages are 0.63, 0.85 and 0.89 for shortest-path closeness, betweenness and current-flow closeness, respectively.

![Correlation with $c_{\text{A}}$ for complex networks.[]{data-label="fig:corr_degrees"}](corr_degrees.pdf){width="50.00000%"}

The results are very different for street networks (Figure \[fig:corr\_degrees\_streets\]). Here the correlation with the degrees is generally very low and sometimes even negative, with an average of -0.02 for closeness, 0.35 for betweenness and 0.07 for current-flow closeness.

While $c_{\text{A}}$ and $c_{\text{ER}}$ are unrelated on street networks, our results show that there is actually a strong correlation between them in complex networks. This behavior is likely due to the different type of degree distributions in the two network classes. While complex networks usually feature a skewed degree distribution with many small, but also some high-degree nodes, the degrees in street networks are closely concentrated.

[8cm]{}[l r]{} Graph & Time approximation \[s\]\
cit-Patents & 125.99\
com-youtube & 5.06\
soc-Epinions1 & 0.28\
com-amazon & 2.56\
hollywood2009 & 107.52\
com-dblp & 1.90\
LiveJournal & 287.27\
Slashdot0902 & 0.41\
roadNet-TX & 6.28\
luxembourg & 0.13\
belgium & 2.39\
netherlands & 5.33\
italy & 10.91\
great-britain & 12.72\
europe & 103.34\

Consequently, in some applications where a very good accuracy is not needed, $c_{\text{A}}$ might be used as an approximation of $c_{\text{ER}}$ in complex networks. The same thing can be said for betweenness, which is only slightly less correlated with $c_{\text{A}}$ than $c_{\text{ER}}$. This is very convenient, since $c_{\text{A}}$ can be computed in ${\mathcal{O} \bigl (m \bigr )}$ time. However, our results in Section \[sec:approx\_comp\] show that our sampling-based approach can compute an extremely accurate approximation in time ${\mathcal{O} \bigl (k m \log (1/\epsilon) \bigr )}$ even when the number $k$ of samples is very small. For this reason, we believe the sampling approach is probably the best option for most applications.

![Correlation with $c_{\text{A}}$ for street networks.[]{data-label="fig:corr_degrees_streets"}](corr_degrees_streets.pdf){width="50.00000%"}

Conclusions
===========

Although many important graph properties can be formulated in terms of Laplacian linear systems, popular network analysis frameworks lack an implementation of state-of-the-art solvers. We bridge this gap by providing a Lean Algebraic Multigrid (LAMG) implementation in NetworKit, our toolkit for large-scale network analysis.

Based on our new LAMG implementation and our algorithms <span style="font-variant:small-caps;">Sampling</span> and <span style="font-variant:small-caps;">Projection</span>, we have computed current-flow closeness centrality and provided the first published results on its behavior on large real-world networks. Our algorithms lead to very accurate results and, thanks to them, we are now able to compute an estimation of current-flow closeness of a subset of nodes on networks with tens of millions of nodes and edges within a few seconds or minutes. In our experiments current-flow closeness alleviates two known problems of shortest-path closeness and can thus be seen as a viable alternative in many scenarios. We have also shown empirically that there is a strong correlation between degrees and both current-flow closeness and betweenness centrality in complex networks, whereas the degree and current-flow closeness are basically unrelated in street networks.

Our approach based on <span style="font-variant:small-caps;">Sampling</span> or <span style="font-variant:small-caps;">Projection</span> is very fast in scenarios where we only need to compute the current-flow closeness for a subset of nodes. However, it might become too expensive if closeness has to be computed for all nodes. In these scenarios an interesting aspect of future work is whether our approximation would still be the best approach or whether inverting the Laplacian of the matrix would be faster in this case.

Instances used for the experiments
==================================

\[table:small\_instances\]

[l r r X r]{} Graph & \#Nodes in LCC & \#Edges in LCC & Description & Ref.\
PGP & 10680 & 24316 & PGP trust network &  [@BaderMSW12dimacs]\
advogato& 5272 & 42816 & Advocato trust network & [@Lasagne]\
Drosophila\_melanogaster & 10424 & 40660 & Interactome &  [@Lasagne]\
Caenorhabditis\_elegans & 4428 & 9659 & Metabolic network &  [@Lasagne]\
CA-HepTh & 8638 & 24806 & Collaboration Network & [@SNAP]\
HC-BIOGRID & 4039 & 10321 & Genetic interaction & [@Lasagne]\
hprd\_pp & 9219 & 36900 & Human proteine interaction &  [@Lasagne]\
Mus\_musculus & 3745 & 5170 & Interactome &  [@Lasagne]\
GoogleNw & 15763 & 148585 & Hyperlinks between web pages & [@Lasagne]\
Homo\_sapiens & 13478 & 61006 & Metabolic network &  [@Lasagne]\
oregon2\_010526 & 11461 & 32730 & AS peering network & [@SNAP]\
as-caida20071105 & 26475 & 53381 & CAIDA AS relationships & [@SNAP]\

\[table:instances\]

[l r r X r]{} Graph & \#Nodes in LCC & \#Edges in LCC & Description & Ref.\
cit-Patents & 3764117 & 16511740 & Citation Network & [@SNAP]\
com-Amazon & 334863 & 925872 & Amazon Product Network & [@SNAP]\
com-DBLP & 317080 & 1049866 & Collaboration Network & [@SNAP]\
com-Youtube & 1134890 & 2987624 & Youtube Social Network & [@SNAP]\
hollywood-2009 & 1069126 & 56306653 & Collaboration Network & [@Davis2011]\
com-LiveJournal & 3997962 & 34681189 & LiveJournal Social Network & [@SNAP]\
Slashdot0902 & 82168 & 504230 & Slashdot Zoo Social Network & [@SNAP]\
soc-Epinions1 & 75877 & 405739 & Epinions Social Network & [@SNAP]\
roadNet-TX & 1351137 & 1879201 & Road Network of Texas & [@SNAP]\
luxembourg.osm & 114599 &119666 & Road Network of Luxembourg &  [@BaderMSW12dimacs]\
belgium.osm & 1441295 & 1549970 & Road Network of Belgium &  [@BaderMSW12dimacs]\
netherlands.osm & 2216688 & 2441238 & Road Network of the Netherlands &  [@BaderMSW12dimacs]\
italy.osm & 6686493 & 7013978 & Road Network of Italy &  [@BaderMSW12dimacs]\
great\_britain.osm & 7733822 & 8156517 & Road Network of Great Britain &  [@BaderMSW12dimacs]\
europe.osm & 50912018 & 54054660 & Road Network of Europe &  [@BaderMSW12dimacs]\

Additional experimental results
===============================

\[table:add\_results\]

[l r r r r]{} Graph & Time exact \[s\] & Time <span style="font-variant:small-caps;">Sampling</span> 20 \[s\] & Spearmann coeff. & Rank Inver.\
PGP & 558.97 & 1.68 & 0.99990 & 0.14%\
advogato & 383.42 & 2.39 & 0.99986 & 0.16%\
Drosophila\_melanogaster & 1077.78 & 3.50 & 0.99986 & 0.12%\
Caenorhabditis\_elegans & 68.50 & 0.64 & 0.99975 & 0.28%\
CA-HepTh & 800.65 & 2.87 & 0.99989 & 0.14%\
HC-BIOGRID & 186.47 & 1.92 & 0.99975 & 0.28%\
hprd\_pp & 988.58 & 4.01 & 0.99990 & 0.10%\
Mus\_musculus & 33.66 & 0.33 & 0.99958 & 0.44%\
GoogleNw & 4612.19 & 8.16 & 0.99987 & 0.14%\
Homo\_sapiens & 1913.06 & 4.90 & 0.99999 & 0.02%\
oregon2\_010526 & 640.97 & 1.49 & 0.99988 & 0.14%\
as-caida20071105 & 3354.62 & 2.62 & 0.99990 & 0.12%\

[^1]: Institute of Theoretical Informatics, Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany

[^2]: Paralution Labs UG & Co. KG, Gaggenau, Germany
