Solving sparse singular linear systems

I have a problem with solving sparse singular linear systems (flow problems on possibly disconnected graphs). I need the minimum norm solution for this. Essentially the sparse version of this:

P = [0.5, -0.5]
M = [1. -1.; -1. 1.]
qr(M, Val(true)) \ P

Is this possible? I don’t need the factorization for anything else.

P = [0.5, -0.5]
M = sparse([1. -1.; -1. 1.])
qr(M) \ P

does work but fails for larger examples and my understanding is that it doesn’t guarantee minimum norm solutions, interestingly it sometimes fails for dense, but works for sparse e.g.:

M2 = zeros(4,4); P2 = rand(4)
M2[1:2,1:2] = [1. -1.;-1. 1.]
M2[3:4,3:4] = [1. -1.;-1. 1.]
MS = sparse(M2)
P2 = M2 * pinv(M2) * P2 # Project on the complement of the kernel (pinv doesn't work for sparse)
qr(M2) \ P2 # LAPACKException(2)
qr(MS) \ P2 # works

Thanks,
Frank

P.S.:

P = [0.5, -0.5]
M = [1. -1.; -1. 1.]
M \ P # SingularException(2)

or

P = [0.5, -0.5]
M = sparse([1. -1.; -1. 1.])
M \ P # SingularException(0)

don’t work (the error messages are not exactly helpful here, is there any way to help with that?)

Your examples are square and symmetric. Will all your use cases have that same structure? As you noticed, sparse QR in Julia doesn’t return the minimum-norm solution. If you’re happy to try an iterative method, have a go with MINRES-QLP: https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/master/src/minres_qlp.jl (symmetric systems only).

Your example M2 is block-diagonal and could be solved as two smaller problems.

If you also have non-symmetric systems (square or rectangular), you can experiment with LSQR and LSMR in https://github.com/JuliaSmoothOptimizers/Krylov.jl.

The iterative solvers LSQR and LSMR that compute the least norm least square solution are also implemented in IterativeSolvers.jl, cf. https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/lsqr/.

Awesome, these are exactly what I was looking for. Many thanks!