# Efficiently computing current flow betweenness centrality for graphs

Definitions:

Given a graph $$G=(V,E),$$ the current flow betweenness is a node-wise measure that captures the fraction of current through a given node with a unit source (s) sink (t) supply $$b_{st}$$ (1 unit of current inserted at node s, $$b_{st}(s)=1$$ and extracted at node t, $$b_{st}(t)=-1,$$ and $$b_{st}(v)=0$$ for $$v\in V\setminus \{s,t\}$$).

For a fixed s-t pair, the throughput $$\tau$$ of a node $$v$$ is given by:

$$\tau_{st}(v)=\frac{1}{2}\left(-|b_{st}(v)|+\sum_{e\ni v}|I(e)|\right) \tag{1}$$

where $$b_{st}$$ is the supply function defined above for the given $$s,t$$ pair, $$I(e)$$ is the current flowing through edge $$e,$$ and $$e\ni v$$ means all edges incident on vertex $$v$$ (i.e. $$v$$ is part of, irrespective of it being at tail or head of edge).

Now the current flow betweenness centrality of a node $$v$$ is simply a normalized sum over all its throughput for all possible supplied pairs $$s,t,$$ i.e.:

$$c(v)=\frac{1}{(n-1)(n-2)} \sum_{s,t\in V}\tau_{s,t}(v) \tag{2}.$$

My implementation of the current-flow betweenness centrality goes as follows:

• Given a graph $$G,$$ I compute its incidence matrix b, corresponding Laplacian lap, and its inverse in S only once at the begining.
• Then I have a module which takes n ($$n=|V|$$), b, S, conductances, supply nodes s,t and returns the list of currents through edges for the given $$s,t$$ pair as supply.
• Then I have module that computes $$\tau_{st}$$ given in $$(1),$$ in which I use a piecewise function for supply $$b_{st},$$ and use Total[] to compute the sum in $$(1).$$
• Then I have a module that computes $$c$$ given in $$(2),$$ where I use a Table to compute $$\tau$$ of $$v$$ for all possible $$s,t$$ and then again use Total to sum them.
• Finally, to compute $$c$$ for all nodes I create a table that runs over all nodes and calls the module for $$c.$$

Actual implementation with a dummy random graph to showcase:

SeedRandom[123]
n = 15;
m = 20;
G = RandomGraph[{n, m}, VertexLabels -> "Name"]
edges = EdgeList[G];

GDirected =
Graph[Range[n], Map[#[[1]] -> #[[2]] &, edges],
VertexLabels -> "Name"]
conductances = ConstantArray[1., m];
b = -1.*Transpose[IncidenceMatrix[GDirected]];
lap = b\[Transpose].DiagonalMatrix[SparseArray[conductances]].b;
a = SparseArray[ConstantArray[1., {1, n}]];
A = ArrayFlatten[{{lap, a\[Transpose]}, {a, 0.}}];
S = LinearSolve[A];
\[Epsilon] = 1. 10^-8;
s = 1;
t = 2;


Edge current module:

edgecurrents[ncount_, invertedkirch_, incid_, conducarr_, nodei_,
nodej_, threshold_] :=
Module[{n = ncount, solver = invertedkirch, incidmat = incid,
G = conducarr, source = nodei, sink = nodej, eps = threshold},
appliedcurr = 1.;
J = SparseArray[{{source}, {sink}} -> {appliedcurr, -appliedcurr}, \
{n}, 0.];
psi = solver[Join[J, {0.}]][[;; -2]];
edgecurr = G incidmat.psi;
(*define current threshold to take care of small values*)

foundcurrents = Threshold[edgecurr, eps];
Return[foundcurrents, Module];
];


$$\tau$$ module:

tau[edgels_, currls_, source_, sink_, vertex_] :=
Module[{edges = edgels, iedges = currls, s = source, t = sink,
v = vertex},
bst[u_, so_, to_] := Piecewise[{{1., u == so}, {-1., u == to}}, 0.];
If[s == t,
res = 0.,
incidv =
Flatten[Position[
edges, (v \[UndirectedEdge] _ | _ \[UndirectedEdge] v)]];
If[incidv == {},
inoutcurrs = 0.;
,
inoutcurrs = Total[Abs[Part[iedges, incidv]]];
];
res = 0.5*(-Abs[bst[v, s, t]] + inoutcurrs);
];
Return[res, Module];
];


$$c$$ module:

currinbet[vcount_, edgels_, conduc_, vertex_, threshold_] :=
Module[{n = vcount, edges = edgels, conducmat = conduc, v = vertex,
eps = threshold},
taust =
Table[tau[edges, edgecurrents[n, S, b, conducmat, s, t, eps], s,
t, v], {s, n}, {t, n}];
ccb = Total[taust, 2]/((n - 1)*(n - 2));
Return[ccb, Module];
];


Example of currents for $$s=1, t=2:$$

edgecurrents[n, S, b, conductances, s, t, \[Epsilon]]
{0.640145, 0.359855, -0.0198915, -0.200723, -0.039783, -0.640145, \
-0.0994575, -0.0144665, 0., 0.0144665, -0.0198915, -0.0433996, \
0.0578662, -0.0144665, 0.359855, -0.359855, 0.101266, -0.0596745, 0., \
0.}


and computing the current-flow betweenness for all nodes:

vccb = Threshold[
Table[currinbet[n, EdgeList[G], conductances, i, \[Epsilon]], {i, 1,
n}], \[Epsilon]]

{0.182869, 0.403493, 0.268327, 0.052163, 0.253522, 0.240516, \
0.524532, 0.135177, 0., 0.208672, 0.275441, 0., 0., 0.282883, \
0.246786}


The obtained results are cross-checked with the existing Python library Networkx for computing $$c$$ and they are in perfect agreement. But sadly efficiency wise, I am doing terribly.

Improved notebook version after Henrik Schumacher's suggestions can be downloaded here, with a working example.

Questions:

• I (think) have minimized the current through edge calculations since S is simply pre-computed, thanks to Henrik Schumacher's approach here. However, I have the feeling I might be doing some things terribly inefficiently from then onward, as my routine slows down drastically for larger graphs. Is there anywhere I could be doing things much more efficiently?

• Is my module-based approach or use of tables also responsible for part of the slow-down?

• Maybe one line of optimization would be to cast $$(1)$$ and $$(2)$$ into linear-algebraic computations to speed them up, but I currently do not see how to do so.

(Any general feedback for rendering the code more efficient is most welcome of course.)

• It is certainly not the bottle neck, but Return[foundcurrents, Module]; is quite unorthodox and I do not know why you use it (I actually did not know that Return allows a second argument). Just writing foundcurrents (without semicolon) instead of Return[foundcurrents, Module]; would be the Mathematica-way of returning a result of a function... – Henrik Schumacher Nov 5 '19 at 15:17

One potential bottleneck is

incidv = Flatten[Position[edges, (v \[UndirectedEdge] _ | _ \[UndirectedEdge] v)]]


as it involves (i) a search in the rather long list of edges and (ii) pattern matching, which both tend to be rather slow.

A quicker way will be to compute all these lists at once via

vertexedgeincidences = IncidenceMatrix[G]["AdjacencyLists"];


and to access the v-th one like this:

incidv = vertexedgeincidences[[v]]


The numbers

inoutcurrs = Total[Abs[Part[iedges, incidv]]];


can also all be computed at once for all v. This can be done with the help if the incidence matrix

B = IncidenceMatrix[G];


via

B.Abs[iedges]


As a general suggestion: Whenever you find yourself evaluating a Sum or Total of something, try to reprase it into Dot-products of vectors, matrices, etc.

• Hi Henrik, sorry for the late reply, didn't notice how time went by... Really helpful recommendations, incorporating them reduced the absolutetiming by a third so far! I also call the current calculator function half as often as the matrix is symmetric, so now my updates sub-modules are: tau[s_, t_, v_, inoutcurrs_] := With[{bst = Piecewise[{{1., v == s}, {-1., v == t}}, 0.]}, res = 0.5*(-Abs[bst] + inoutcurrs[[v]])]; and currinbet is (...) – user52181 Nov 13 '19 at 13:46
• (...) I edited the post and added a link to the latest version (a notebook with a working example), which is much simpler to show what I have now instead of pasting code again. Do you reckon my main bottleneck lies in the current calculation still? (despite being cut in half), thanks again for all your feedback. – user52181 Nov 13 '19 at 14:51