# Matrix division \\ for sparse and singular matrices

**URL:** <https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253>\
**Category:** General Usage\
**Tags:** linearalgebra, sparse\
**Created:** [September 16, 2021, 10:52am UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253 "2021-09-16T10:52:36Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [September 16, 2021, 10:52am UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/1 "2021-09-16T10:52:36Z")

</div>

It seems that the behaviour of matrix division \ is inconsistent if sparse matrices are used:

```julia
using LinearAlgebra, SparseArrays

A1=[1 1 ;0 0.0]

B1=[4; 0.0]

x1=A1\B1 #Error: LAPACKException(2)

x2=sparse(A1)\sparse(B1) #Error: SingularException(2)

x3=sparse(A1)\B1 #Error: SingularException(2)

x4=A1\sparse(B1) #works. Why, if A\B doesn't?

A2=[1 0 ;0 0.0]

B2=[4; 0.0]

x1=A2\B2 #gives 4 and NaN - why if A is still triangular and singular?

x2=sparse(A2)\sparse(B2) #gives 4 and NaN

x3=sparse(A2)\B2 #gives 4 and NaN

x4=A2\sparse(B2) #gives 4 and NaN

A3=[1 1; 1 1]
B3=[4; 4]

x1=A3\B3 #Error: SingularException(2)

x2=sparse(A3)\sparse(B3) #Error: MethodError: no method matching lu! #why a different error?

x3=sparse(A3)\B3 #Error: MethodError: no method matching lu!

x4=A3\sparse(B3) #Error: SingularException(2)

```

Shouldn’t \ work for sparse matrices? Also, why no least-square solutions for a quadratic matrix (as described in [Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra))

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 16, 2021, 11:20am UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/2 "2021-09-16T11:20:53Z")

</div>

You ask why

> [@blob](#):
>
> `x4=A1\sparse(B1)`

works, while all other combinations of imposed sparsity do not. The reason is perhaps related to the following issue: [Request: SPQR backslash operator for sparse right hand sides · Issue #33637 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/33637). It says that Julia does not support sparsity on the right hand side (because the underlying engine - SuitSparse) does not. Therefore in this case the right hand side is kept dense and `A1\x1` is solved. The issue is from 2019 and I have not checked if anything changed, but I guess that it still holds.

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [September 16, 2021, 11:32am UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/3 "2021-09-16T11:32:00Z")

</div>

Yes, that’s an idea, maybe `sparse(B1)` gets converted into just `B1`. But the if it should work for dense matrices, I don’t get why this fails:  
`x1=A1\B1 `

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 16, 2021, 11:37am UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/4 "2021-09-16T11:37:10Z")

</div>

> [@Why does julia systematically fails when doing \\ operation?](https://discourse.julialang.org/t/why-does-julia-systematically-fails-when-doing-operation/67242/9):
>
> It appears that the convention adopted by Julia is that the backslash symbol used in x = A\b with a square A matrix really stands exclusively for the inverse of a matrix (computed perhaps using LU decomposition). If A is singular (or detected close to singular), the call of x= A\b fails. In such situations, SVD is your friend: julia\> G = svd(A) SVD{Float64, Float64, Matrix{Float64}} U factor: 3×3 Matrix{Float64}: -0.332289 0.850246 0.408248 -0.549449 0.177311 -0.816497 -0.766609 -0.4…

Compare with

```julia
julia> x1p = pinv(A1)*B1
2-element Vector{Float64}:
 1.9999999999999993
 1.9999999999999993

```

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [September 16, 2021, 11:54am UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/5 "2021-09-16T11:54:47Z")

</div>

If `x4=A1\sparse(B1)` converts `B1` into a dense matrix, then `A1\sparse(B1)` should be equivalent to `A1\B1`. But in my case, `A1\B1` fails whereas `A1\sparse(B1)` gives a solution. So it seems that the ‘converted B1’ is different somehow - and I don’t know how and why.

I know I can get a solution using other methods, like `pinv`, the point is that I would like to avoid adjusting my solver to every possible matrix that I could get.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 16, 2021, 1:07pm UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/6 "2021-09-16T13:07:41Z")

</div>

> [@Why does julia systematically fails when doing \\ operation?](https://discourse.julialang.org/t/why-does-julia-systematically-fails-when-doing-operation/67242/9):
>
> In such situations, SVD is your friend

QR is a _much_ cheaper way to get least-square solutions of singular/nonsquare systems: `qr(A) \ b`

Also, QR works for sparse matrices too.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 16, 2021, 1:16pm UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/7 "2021-09-16T13:16:06Z")

</div>

True. Actually I do mention it in that post.

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [September 16, 2021, 3:26pm UTC](https://discourse.julialang.org/t/matrix-division-for-sparse-and-singular-matrices/68253/8 "2021-09-16T15:26:30Z")

</div>

I see. You are right. The inconsistency that you mention initially escaped my attention. Indeed, `A\b` and `A\sparse(b)` behave rather inconsistently. For convenience of other readers I include the code again (well, I changed the names of the variables a bit):

```julia
julia> using LinearAlgebra, SparseArrays

julia> A=[1 1 ;0 0.0]
2×2 Matrix{Float64}:
 1.0 1.0
 0.0 0.0

julia> b=[4; 0.0]
2-element Vector{Float64}:
 4.0
 0.0

julia> x1=A\b
ERROR: LAPACKException(2)
Stacktrace:
 [1] chklapackerror(ret::Int64)
   @ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:38
 [2] trtrs!(uplo::Char, trans::Char, diag::Char, A::Matrix{Float64}, B::Vector{Float64})
   @ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:3426
 [3] ldiv!
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:739 [inlined]
 [4] \
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:1673 [inlined]
 [5] \(A::Matrix{Float64}, B::Vector{Float64})
   @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1134
 [6] top-level scope
   @ none:1

julia> x2=sparse(A)\sparse(b)

ERROR: SingularException(2)
Stacktrace:
 [1] _ldiv!(U::UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}, B::Vector{Float64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:706
 [2] ldiv!(A::UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}, B::Vector{Float64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:632
 [3] \(A::UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}, b::SparseVector{Float64, Int64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsevector.jl:1824
 [4] \(A::SparseMatrixCSC{Float64, Int64}, B::SparseVector{Float64, Int64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1562
 [5] top-level scope
   @ none:1

julia> x3=sparse(A)\b
ERROR: SingularException(2)
Stacktrace:
 [1] _ldiv!(U::UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}, B::Vector{Float64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:706
 [2] ldiv!(A::UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}, B::Vector{Float64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:632
 [3] \
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:1673 [inlined]
 [4] \(A::SparseMatrixCSC{Float64, Int64}, B::Vector{Float64})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1562
 [5] top-level scope
   @ none:1

julia> x4=A\sparse(b) 
2-element Vector{Float64}:
 4.0
 0.0

```