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:

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 = 
       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., \

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, \

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.


  • 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.)

| improve this question | | | | |
  • $\begingroup$ 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... $\endgroup$ – 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];



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.

| improve this answer | | | | |
  • $\begingroup$ 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 (...) $\endgroup$ – user929304 Nov 13 '19 at 13:46
  • $\begingroup$ (...) 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. $\endgroup$ – user929304 Nov 13 '19 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.