Solving sparse singular linear systems

Specific Domains Numerics
1

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

2

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.

1 Like
3

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

4

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