performance tuning - 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.)

Answer

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.