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

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.

resitanceDistanceMatrix[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[123];
idx = RandomSample[1 ;; VertexCount[g], 5];
resitanceDistanceMatrix[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".

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