# Solving sparse singular linear systems

**URL:** <https://discourse.julialang.org/t/solving-sparse-singular-linear-systems/33401>\
**Category:** Numerics\
**Created:** [January 15, 2020, 3:49pm UTC](https://discourse.julialang.org/t/solving-sparse-singular-linear-systems/33401 "2020-01-15T15:49:41Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![FHell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fhell/32/2421_2.png) [@FHell](https://discourse.julialang.org/u/FHell)\
**Post date:** [January 15, 2020, 3:49pm UTC](https://discourse.julialang.org/t/solving-sparse-singular-linear-systems/33401/1 "2020-01-15T15:49:41Z")

</div>

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:

```julia
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.

```julia
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.:

```julia
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.:

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

```

or

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

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [January 18, 2020, 3:40pm UTC](https://discourse.julialang.org/t/solving-sparse-singular-linear-systems/33401/2 "2020-01-18T15:40:22Z")

</div>

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](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](https://github.com/JuliaSmoothOptimizers/Krylov.jl).

---

<div class="post-metadata">

**Author:** ![ranocha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ranocha/32/35588_2.png) [@ranocha](https://discourse.julialang.org/u/ranocha)\
**Post date:** [January 19, 2020, 12:18pm UTC](https://discourse.julialang.org/t/solving-sparse-singular-linear-systems/33401/3 "2020-01-19T12:18:12Z")

</div>

The iterative solvers LSQR and LSMR that compute the least norm least square solution are also implemented in [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl), cf. [https://juliamath.github.io/IterativeSolvers.jl/dev/linear\_systems/lsqr/](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/lsqr/).

---

<div class="post-metadata">

**Author:** ![FHell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fhell/32/2421_2.png) [@FHell](https://discourse.julialang.org/u/FHell)\
**Post date:** [January 24, 2020, 1:23pm UTC](https://discourse.julialang.org/t/solving-sparse-singular-linear-systems/33401/4 "2020-01-24T13:23:19Z")

</div>

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