Julia: Request: SPQR backslash operator for sparse right hand sides

[TomDeWeer]

I would like to invert a sparse matrix using LU factorization, via the following code:

using LinearAlgebra
using SuiteSparse
using SparseArrays
function invertMatrix(A)
    Af = lu(A)
    I = spdiagm(0 => ones(size(A)[1]))
    Ainv =  Af \ I
end 

The factorization Af is is a SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64} object and supports the backslash operator via the SPQR module, according to https://julialinearalgebra.github.io/SuiteSparse.jl/latest/umfpack.html. This backslash operator only supports right hand sides that inherit from StridedVecOrMat (see https://julialinearalgebra.github.io/SuiteSparse.jl/latest/spqr.html). Supplying the sparse identity matrix I constructed above thus fails with a method error. I therefore request that the SPQR backslash operator is extended to work with sparse right hand sides. The return type is also a dense array, so ideally this also becomes sparse array.

Kind regards,
Tom

PS: I know I can just make my identity matrix dense, but this causes excessive memory requirements to save all the zeros. Writing a for loop to go around this is my current solution, but is also suboptimal.

[TomDeWeer]

I did some research and SuiteSparse has no solver that can incorporate sparse right hand sides without converting it to dense. Since Julia relies on SuiteSparse, this request is therefore improperly placed. Closing this issue.

[ViralBShah]

I believe that is accurate - SuiteSparse doesn't yet support sparse RHS (although it may be possible for KLU to support it at some point). One can do various tricks to speed up sparse RHS solutions, but I believe you can't get the benefits of the multifrontal or supernodal methods usually in those cases.

The best way to do all this would be to implement the Gilbert Peierls algorithm natively in Julia: https://ecommons.cornell.edu/bitstream/handle/1813/6623/86-783.pdf

[ViralBShah]

Could be a Google Summer of Code project potentially for someone with the right level of background and interest.

[DrTimothyAldenDavis]

CSparse has a cs_spsolve function that can do both x=L\b and x=U\b, both using Gilbert-Peierl's method, and where b is sparse. It assumes a workspace of size n, and takes O(n) time to initialize it. If the workspace were passed in pre-initialized, and then the parts used were cleared before returning, then it could be used to create a sparse forward/backsolve to compute a single column of the inverse. That function would be a good place to start since it's so simple.

[ViralBShah]

@DrTimothyAldenDavis - So great to see you here!