Matrix division \ for sparse and singular matrices

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

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)

You ask why

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

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

Compare with

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

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.

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.

1 Like

True. Actually I do mention it in that post.

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> 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
1 Like