LinearAlgebra help

[tushargayan2324]

I am trying to solve a certain question;

Question: There is a 2-d square grid of nodes, with resistors placed randomly to connect those nodes. And eventually, those connections may extend from one end to the other end. Say for now, I am only considering left and right side. I want to know what is the total equivalent resistance of the network, which leads the current from one end to the other.

(For simplicity, consider only one main cluster to extend end to end)

Now, I saw this code for uniform grid case: Resistor mesh - Rosetta Code

using SparseArrays
N = 3
D1 = sparse(I,N-1,N) - spdiagm(N-1,N, 1=>ones(N-1))
D = [ kron(D1, sparse(I,N,N)); kron(sparse(I,N,N), D1) ]
i, j = 1 , 9
b = zeros(N^2); b[i], b[j] = 1, -1
v = (D' * D) \ b
v[i] - v[j]

And made me think I can maybe use this for random network too.

I tried to make an Incidence matrix from a given random lattice (from my previous attempt at percolation problem). And, formed the Laplacian matrix using Incidence matrix.

For sanity check, I tried to make 3x3 grid with all connections. This is what I got for Incidence matrix.

12×9 Matrix{Float64}:
 -1.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0  0.0
 -1.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0  0.0
  0.0  -1.0   0.0   0.0   1.0   0.0   0.0   0.0  0.0
  0.0  -1.0   1.0   0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0  -1.0   0.0   0.0   1.0   0.0   0.0  0.0
  0.0   0.0   0.0  -1.0   0.0   0.0   1.0   0.0  0.0
  0.0   0.0   0.0  -1.0   1.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0  -1.0   0.0   0.0   1.0  0.0
  0.0   0.0   0.0   0.0  -1.0   1.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0  -1.0   0.0   0.0  1.0
  0.0   0.0   0.0   0.0   0.0   0.0  -1.0   1.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0  -1.0  1.0

And the following is the Laplacian matrix

9×9 Matrix{Float64}:
  2.0  -1.0   0.0  -1.0   0.0   0.0   0.0   0.0   0.0
 -1.0   3.0  -1.0   0.0  -1.0   0.0   0.0   0.0   0.0
  0.0  -1.0   2.0   0.0   0.0  -1.0   0.0   0.0   0.0
 -1.0   0.0   0.0   3.0  -1.0   0.0  -1.0   0.0   0.0
  0.0  -1.0   0.0  -1.0   4.0  -1.0   0.0  -1.0   0.0
  0.0   0.0  -1.0   0.0  -1.0   3.0   0.0   0.0  -1.0
  0.0   0.0   0.0  -1.0   0.0   0.0   2.0  -1.0   0.0
  0.0   0.0   0.0   0.0  -1.0   0.0  -1.0   3.0  -1.0
  0.0   0.0   0.0   0.0   0.0  -1.0   0.0  -1.0   2.0

Which is same output as the rosseta code;

9×9 SparseMatrixCSC{Float64, Int64} with 33 stored entries:
  2.0  -1.0    ⋅   -1.0    ⋅     ⋅     ⋅     ⋅     ⋅ 
 -1.0   3.0  -1.0    ⋅   -1.0    ⋅     ⋅     ⋅     ⋅ 
   ⋅   -1.0   2.0    ⋅     ⋅   -1.0    ⋅     ⋅     ⋅ 
 -1.0    ⋅     ⋅    3.0  -1.0    ⋅   -1.0    ⋅     ⋅ 
   ⋅   -1.0    ⋅   -1.0   4.0  -1.0    ⋅   -1.0    ⋅ 
   ⋅     ⋅   -1.0    ⋅   -1.0   3.0    ⋅     ⋅   -1.0
   ⋅     ⋅     ⋅   -1.0    ⋅     ⋅    2.0  -1.0    ⋅ 
   ⋅     ⋅     ⋅     ⋅   -1.0    ⋅   -1.0   3.0  -1.0
   ⋅     ⋅     ⋅     ⋅     ⋅   -1.0    ⋅   -1.0   2.0

But, the problem comes when I try to solve the above linear equation;

M = A' * A
b = zeros(9)
b[1] = 1
b[6] = -1
x = M \ b

The result is somehow, always SingularException(some number) [numbers I remember on top of my head are 0, 1, 13, 9, 2 , etc]

Can someone please tell me whats going on? Because essentially, both the laplacian matrices are same, but this code doesn’t yield the result. Is it because of Sparse matrix that works in the rosetta code?

(Now, when I use the sparse(M), it does give me the result of the completely connected grid problem. But, then why does the normal matrix give error?)

Edit: Typo in last code block

[jishnub]

Could you post the full code (what are A and M?)

[tushargayan2324]

A and M are Incidence and Laplacian matrix respectively. The full code is irrelevant in this case otherwise I would upload it.

[Dan]

I think you forgot to add ones(N,N)/N to the Laplacian in the equation (Based on Wikipedia article)

Resistance distance

[stevengj]

@tushargayan2324

The result is somehow, always SingularException(some number)

That’s because your matrix M is indeed singular — if you multiply it by the vector x = ones(9) # = [1, 1, ..., 1] you’ll see that M*x is indeed zero. Mx = b is still solvable for right-hand sides that are orthogonal to this vector, i.e. which sum to zero, which is indeed the case for the b vector that you have here, but the solution is not unique (you can add any multiple of ones(9) to the solution), and it means that you need to take some care in how you solve the linear system.

This is typical of current/voltage problems in which you don’t set a ground, because the voltage is then not unique — you can add any overall constant — and the constraint that the right-hand side must sum to zero corresponds to the physical constraint that the net current flowing into your circuit must be zero (since there is no ground for it to escape to).

You can solve this in a variety of ways, but the simplest is probably to force one of your voltages to be zero — connect one node (assuming a connected graph) to “ground”. For example, if you set x[end] to be zero, it corresponds to simply deleting the last row and column of your matrix, i.e. solve the non-singular system

x = [M[1:end-1, 1:end-1] \ b[1:end-1]; 0]

(Note that the current injected to/from ground is no longer under your control, so the last row of b must be discarded: effectively, it is assumed to satisfy b[end] = -sum(b[1:end-1]), which is the case in your example.)

@Dan

I think you forgot to add ones(N)/N to the Laplacian in the equation

I think you mean ones(N, N) / N? Yes, another solution is to add any positive constant to all the elements of M, i.e. replace M with M .+ 1 or M .+ 1/N.

The reason this works is that the original Mx = b equation has one solution in which sum(x) == 0 (assuming sum(b) == 0 so that it has a solution at all, as noted above). (Again, as noted above, you can add any constant to solutions x, so you can always choose this constant so that the sum is zero.) For this solution, however, \Phi x = 0 where \Phi is the matrix of all 1’s (using the notation of the Wikipedia article), i.e. ones(N, N). So, for this solution Mx = b \implies (M + \alpha \Phi) x = b for any scalar \alpha. By using any \alpha > 0 the system becomes positive-definite, and hence this x is the only solution.

So, in short, you can do x = (M .+ 1) \ b and it will give the unique solution x whose components sum to zero, assuming sum(b) == 0. (Which is another valid solution, that will differ by an overall additive constant, compared to my method of choosing a “ground” node above. Note that all valid solutions will have the same voltage differences.)

But this is not a great solution technique if M is sparse, because M .+ 1 is non-sparse. One way to get around this fact might be to use the Sherman–Morrison formula, since the matrix of all 1’s is rank 1, but that leads back again to the same issue of M being singular. Easier to choose a ground, in my opinion, at least for direct methods. (For iterative methods you can multiply by M and \Phi quickly and then add the results.)