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