Solving "Resistance between two nodes on a grid" problem in Mathematica

In the context of resistor networks and finding the (equivalent) resistance between two arbitrary nodes, I am trying to learn how to write a generic approach in Mathematica, generic as in an approach that also lends itself to large spatially randomly distributed graphs as well (not just lattices), where then one has to deal with sparse matrices. Before getting there, I've tried simply recreating a piece of algorithm written in Julia for solving an example on a square grid, with all resistances set to 1.

Here's the grid where each edge depicts a resistor between its incident nodes (all resistance values are assumed to be $$1 \Omega$$) and two arbitrary nodes ($$A$$ at {2,2} and $$B$$ at {7,8}) are highlighted, question is to find the resistance between them. In the Julia's code snippet, the approach of injecting a current and measuring the voltages at the two nodes is adopted, as shown below: (source)

N = 10
D1 = speye(N-1,N) - spdiagm(ones(N-1),1,N-1,N)
D = [ kron(D1, speye(N)); kron(speye(N), D1) ]
i, j = N*1 + 2, N*7+7
b = zeros(N^2); b[i], b[j] = 1, -1
v = (D' * D) \ b
v[i] - v[j]

Output: 1.6089912417307288

I've tried to recreate exactly the same approach in Mathematica, here's what I have done:

n = 10;
grid = GridGraph[{n, n}];
i = n*1 + 2;
j = n*7 + 7;
b = ConstantArray[0, {n*n, 1}];
b[[i]] = {1};
b[[j]] = {-1};
incidenceMat = IncidenceMatrix[grid];
matrixA = incidenceMat.Transpose[incidenceMat];
v = LinearSolve[matrixA, b]

I feel very silly, but I must be missing something probably very obvious as LinearSolve does not manage to find a solution (for the chosen nodes the answer is know to be $$1.608991...$$, which is obtained by taking the potential difference between A and B since the current is set to 1).

Questions

• Have I mis-interpreted something in my replication of the algorithm sample written in Julia?

• It would be very interesting and useful if someone could comment on how extensible these approaches are to more general systems (2d, 3d and not only for lattices). For instance, which approaches would be more suitable to adopt in Mathematica for larger resistor networks (in terms of efficiency, as one would have to deal with very sparse matrices probably).

As a side-note, on the same Rosetta article, there are two alternative code snippets provided for Mathematica (which follows Maxima's approach, essentially similar to the one written Julia). In case someone is interested I include them here: (source for both)

gridresistor[p_, q_, ai_, aj_, bi_, bj_] :=
Block[{A, B, k, c, V}, A = ConstantArray[0, {p*q, p*q}];
Do[k = (i - 1) q + j;
If[{i, j} == {ai, aj}, A[[k, k]] = 1, c = 0;
If[1 <= i + 1 <= p && 1 <= j <= q, c++; A[[k, k + q]] = -1];
If[1 <= i - 1 <= p && 1 <= j <= q, c++; A[[k, k - q]] = -1];
If[1 <= i <= p && 1 <= j + 1 <= q, c++; A[[k, k + 1]] = -1];
If[1 <= i <= p && 1 <= j - 1 <= q, c++; A[[k, k - 1]] = -1];
A[[k, k]] = c], {i, p}, {j, q}];
B = SparseArray[(k = (bi - 1) q + bj) -> 1, p*q];
LinearSolve[A, B][[k]]];
N[gridresistor[10, 10, 2, 2, 8, 7], 40]

Alternatively:

graphresistor[g_, a_, b_] :=
LinearSolve[
SparseArray[{{a, a} -> 1, {i_, i_} :> Length@AdjacencyList[g, i],
Alternatives @@ Join[#, Reverse /@ #] &[
List @@@ EdgeList[VertexDelete[g, a]]] -> -1}, {VertexCount[
g], VertexCount[g]}], SparseArray[b -> 1, VertexCount[g]]][[b]];
N[graphresistor[GridGraph[{10, 10}], 12, 77], 40]

In addition to Carl Woll's post:

Computing the pseudoinverse of a the graph Laplacian matrix (a.k.a. the KirchhoffMatrix) is very expensive and in general leads to a dense matrix that, if the graph is too large, cannot be stored in RAM. In the case that you have to compute only a comparatively small block of the resistance distance matrix, you can employ sparse methods as follows:

Generating a graph with 160000 vertices.

g = GridGraph[{400, 400}, GraphLayout -> None];
L = N@KirchhoffMatrix[g];

The idea here is that I know in advance that $$\mathbf{A}$$ is symmetric and positive semidefinite and that $$\operatorname{ker}(\mathbf{L}) = \operatorname{im}(\mathbf{L})^\perp = \mathbb{R} \, \mathbf{1}.$$ (The latter holds only if the graph is connected.)

Fix a vector $$\mathbf{b}$$ and denote the Kirchhoff matrix by $$\mathbf{L}$$ and its pseudoinverse by $$\mathbf{L}^\dagger$$. Denote the orthogonal projection of $$\mathbf{b}$$ onto $$\operatorname{im}(\mathbf{L})$$ by $$\mathbf{y}$$, so that we have $$\mathbf{b} = \mathbf{y} + \mathbf{1} \, \lambda$$ with some $$\lambda \in \mathbb{R}$$. The orthogonal projector onto $$\operatorname{im}(\mathbf{L})$$ is given by $$\mathbf{L} \, \mathbf{L}^\dagger$$ so that we have $$\mathbf{y} = \mathbf{L} \, \mathbf{L}^\dagger \, \mathbf{b} = \mathbf{L} \, \mathbf{x}$$. Thus: $$\mathbf{b} = \mathbf{L} \, \mathbf{x} + \mathbf{1} \, \lambda.$$ We have $$\operatorname{ker}(\mathbf{L})^\perp = \operatorname{ima}(\mathbf{L}^\dagger)$$, hence $$\mathbf{1}^\intercal \, \mathbf{x} = \mathbf{1}^\intercal \, \mathbf{L}^\dagger \mathbf{b} = 0,$$ hence $$\mathbf{1}^\intercal \, \mathbf{x} = 0.$$

That is, it suffices to solve the linear saddle point system $$\begin{pmatrix} \mathbf{L} & \mathbf{1} \\ \mathbf{1}^\intercal &0 \end{pmatrix} \begin{pmatrix} \mathbf{x} \\ \lambda \end{pmatrix} = \begin{pmatrix} \mathbf{b} \\ \mathbf{0} \end{pmatrix}.$$ The good things are that the saddle point matrix is (i) invertible and (ii) usually quite sparse. So we may employ LinearSolve to solve this linear system.

The following builds the saddle point matrix A and computes an $$LU$$-factorization S of it (You may read S basically as the inverse of A).

A = With[{a = SparseArray[ConstantArray[1., {1, VertexCount[g]}]]},
ArrayFlatten[{{L, a\[Transpose]}, {a, 0.}}]
];
S = LinearSolve[A]; // AbsoluteTiming

Applying the pseudoinverse of L to a vector b is now equivalent to

b = RandomReal[{-1, 1}, VertexCount[g]];
x = S[Join[b, {0.}]][[1 ;; -2]];

We may exploit that via the following helper function; internally, it computes only few columns of the pseudoinverse and returns the corresponding resistance graph matrix.

resistanceDistanceMatrix[S_LinearSolveFunction, idx_List] :=
Module[{n, basis, Γ},
n = S[[1, 1]];
basis = SparseArray[
Transpose[{idx, Range[Length[idx]]}] -> 1.,
{n, Length[idx]}
];
Γ = S[basis][[idx]];
(* stealing from Carl Woll *)
Outer[Plus, Diagonal[Γ], Diagonal[Γ]] - Γ - Transpose[Γ]
];

Let's compute the resistance distance matrix for 5 random vertices:

SeedRandom;
idx = RandomSample[1 ;; VertexCount[g], 5];
resistanceDistanceMatrix[S, idx] // MatrixForm

$$\left( \begin{array}{ccccc} 0. & 2.65527 & 2.10199 & 2.20544 & 2.76988 \\ 2.65527 & 0. & 2.98857 & 2.85428 & 2.3503 \\ 2.10199 & 2.98857 & 0. & 2.63996 & 3.05817 \\ 2.20544 & 2.85428 & 2.63996 & 0. & 3.04984 \\ 2.76988 & 2.3503 & 3.05817 & 3.04984 & 0. \\ \end{array} \right)$$

This requires $$k$$ linear solves for $$k (k-1) /2$$ distances, so it is even more efficient than the method you posted (which needs one linear solve per distance).

The most expensive part of the code is to generate the LinearSolveFunction S. Thus, I designed the code so that S can be reused.

Under the hood, a sparse LU-factorization is computed via UMFPACK. Since the graph g is planar, this is guaranteed to be very quick compared to computing the whole pseudoinverse.

For nonplanar graphs, things become complicated. Often, using LU-factorization will work in reasonable time. But that is not guaranteed. If you have for example a cubical grid in 3D, LU-factorization will take much longer than a 2D-problem of similar size even if you measure size by the number of nonzero entries. In such cases, iterative linear solvers with suitable preconditioners may perform much better. One such method (with somewhat built-in preconditioner) is the (geometric or algebraic) multigrid method. You can find an implementation of such a solver along with a brief explanation of its working here. For a timing comparison of linear solves on a cubical grid topology see here. The drawback of this method is that you have to create a nested hierarchy of graphs on your own (e.g. by edge collapse). You may find more info on the topic by googling for "multigrid"+"graph".

• Sorry for my late comment, many thanks for this wonderful analysis and extremely useful discussions. You mentioned that this approach works for connected graphs only. To better understand its domain of validity, while the connected requirement is met, would it work e.g. for i) cases where there may be large loops present? ii) random graphs with spatially distributed nodes in 3D? iii) what about networks with random resistance values?(inhomogeneity) As a side-question, given two nodes, do you reckon one can use this method to simply detect edges that are carrying non-zero current between them?
– user52181
Mar 13 '19 at 15:30
• What I meant is the computation of the pseudoinverse: If you have several connected components you have to add a row and a column for each component to the matrix: The entries have to be all 1 for vertices of the given coponent and zero otherwise. However, the resistance distance formula from en.wikipedia.org/wiki/Resistance_distance#Definition is quite likely not correct in this case (the resitance between vertices in different connected components of the graph should be infinite). Mar 13 '19 at 15:38
• At i) I don't see at the moment why loops should be a problem. At ii) The method should work but may become slower. It really depends on the topology of the graph at hand, and thus indirectly also on the graph probability distribution. At iii) There is a notion of a graph Laplacians for weighted graphs. Probably one has to employ the inverses of the edge resitances as weights, but I am not sure about that. At your last question: I don't understand it. Maybe you can explain what you mean? Mar 13 '19 at 15:44
• See also here: match.pmf.kg.ac.rs/electronic_versions/Match50/… I had a short gilmpse into the paper and I found the statement that for computing resistance distance, every generalized inverse will do (and will return the same result). So I guess the method should also apply to networks with varying edge resistances. Mar 13 '19 at 15:51
• IIRC, instead of L = N@KirchhoffMatrix[g], you have to use the matrix L from this post in the case of nonconstant edge resistances. In that post, I defined the vector conductances which is a function on edges; if you set conductances to the inverse of the edge resistances, everything should work out correctly. Mar 19 '19 at 10:46

Based on rcampion2012's answer to Efficient Implementation of Resistance Distance for graphs?, you could use:

resistanceGraph[g_] := With[{Γ = PseudoInverse[N @ KirchhoffMatrix[g]]},
Outer[Plus, Diagonal[Γ], Diagonal[Γ]] - Γ - Transpose[Γ]
]

Then, you can find the resistance using:

r = resistanceGraph[GridGraph[{10, 10}]];
r[[12, 68]]

1.60899

• Nice answer, but why did you take the [[12,68]] part and what are the right terms to take in the general case of n x n system? Mar 13 '19 at 11:04
• @Loki There are no right choices here, one is usually interested to know the equivalent resistance between two arbitrary points of the network, in this example the 12th and 68th nodes are chosen (they correspond to the nodes I had highlighted in the question, and if you collapse the grid into a vector you'd get these indices).
– user52181
Mar 13 '19 at 11:35
• Thanks, I hadn't noticed the small red dots in the graph. Nice answer, my +1 Mar 13 '19 at 11:40