8

is there any equivalent to scipy.sparse.linalg.spsolve in Julia? Here's the description of the function in Python.

In [59]: ?spsolve                                                                                                                                                                                                  
Signature: spsolve(A, b, permc_spec=None, use_umfpack=True)
Docstring:
Solve the sparse linear system Ax=b, where b may be a vector or a matrix.

I couldn't find this in Julia's LinearAlgebra and SparseArrays. Is there anything I miss or any alternatives?

Thanks

EDIT

For example:

In [71]: A = sparse.csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)                                                                                                                                    

In [72]: B = sparse.csc_matrix([[2, 0], [-1, 0], [2, 0]], dtype=float)                                                                                                                                             

In [73]: spsolve(A, B).data                                                                                                                                                                                        
Out[73]: array([ 1., -3.])

In [74]: spsolve(A, B).toarray()                                                                                                                                                                                   
Out[74]: 
array([[ 0.,  0.],
       [ 1.,  0.],
       [-3.,  0.]])

In Julia, with \ operator

julia> A = Float64.(sparse([3 2 0; 1 -1 0; 0 5 1]))
3×3 SparseMatrixCSC{Float64,Int64} with 6 stored entries:
  [1, 1]  =  3.0
  [2, 1]  =  1.0
  [1, 2]  =  2.0
  [2, 2]  =  -1.0
  [3, 2]  =  5.0
  [3, 3]  =  1.0

julia> B = Float64.(sparse([2 0; -1 0; 2 0]))
3×2 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 1]  =  2.0
  [2, 1]  =  -1.0
  [3, 1]  =  2.0

julia> A \ B
ERROR: MethodError: no method matching ldiv!(::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}, ::SparseMatrixCSC{Float64,Int64})
Closest candidates are:
  ldiv!(::Number, ::AbstractArray) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/generic.jl:236
  ldiv!(::SymTridiagonal, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T; shift) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/tridiag.jl:208
  ldiv!(::LU{T,Tridiagonal{T,V}}, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) where {T, V} at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/lu.jl:588
  ...
Stacktrace:
 [1] \(::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}, ::SparseMatrixCSC{Float64,Int64}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/factorization.jl:99
 [2] \(::SparseMatrixCSC{Float64,Int64}, ::SparseMatrixCSC{Float64,Int64}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/linalg.jl:1430
 [3] top-level scope at REPL[81]:1

1 Answer 1

9

Yes, it's the \ function.

julia> using SparseArrays, LinearAlgebra

julia> A = sprand(Float64, 20, 20, 0.01) + I # just adding the identity matrix so A is non-singular.

julia> typeof(A)
SparseMatrixCSC{Float64,Int64}

julia> v = rand(20);

julia> A \ v
20-element Array{Float64,1}:
 0.5930744938331236
 0.8726507741810358
 0.6846427450637211
 0.3135234897986168
 0.8366321472466727
 0.11338490488638651
 0.3679058951515244
 0.4931583108292607
 0.3057947282994271
 0.27481281228206955
 0.888942874188458
 0.905356044150361
 0.17546911165214607
 0.13636389619386557
 0.9607381212005248
 0.2518153541168824
 0.6237205353883974
 0.6588050295549153
 0.14748809413104935
 0.9806131247053784

Edit in response to question edit:

If you want v here to instead be a sparse matrix B, then we can proceed by using the QR decomposition of B (note that cases where B is truly sparse are rare:

function myspsolve(A, B)
    qrB = qr(B)
    Q, R = qrB.Q, qrB.R
    R = [R; zeros(size(Q, 2) - size(R, 1), size(R, 2))]
    A\Q * R
end

now:

julia> A = Float64.(sparse([3 2 0; 1 -1 0; 0 5 1]))
3×3 SparseMatrixCSC{Float64,Int64} with 6 stored entries:
  [1, 1]  =  3.0
  [2, 1]  =  1.0
  [1, 2]  =  2.0
  [2, 2]  =  -1.0
  [3, 2]  =  5.0
  [3, 3]  =  1.0

julia> B = Float64.(sparse([2 0; -1 0; 2 0]))
3×2 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 1]  =  2.0
  [2, 1]  =  -1.0
  [3, 1]  =  2.0


julia> mysolve(A, B)
3×2 Array{Float64,2}:
  0.0  0.0
  1.0  0.0
 -3.0  0.0

and we can test to make sure we did it right:

julia> mysolve(A, B) ≈ A \ collect(B)
true
5
  • Thanks! See the example my updated question, I couldn't seem to get this working in Julia.. Commented Mar 13, 2020 at 0:33
  • Are you sure your B needs to be sparse? I think your best bet is to convert B to dense because almost no matter what, A\B is going to end up being dense anyways and that’s what Julia supports out of the box. There may be a way to do this while keeping B dense but I don’t know.
    – Mason
    Commented Mar 13, 2020 at 1:52
  • Thanks! That's not the real use case for now, but it's from spsolve docstring example. If that's the case then that's fine :) Commented Mar 13, 2020 at 1:53
  • I have updated my answer to show how to solve this in the case where B is sparse. Can you please edit the question title to something a bit more searchable? How about "how to solve a linear system where both inputs are sparse?"
    – Mason
    Commented Mar 13, 2020 at 2:13
  • 1
    For sure! Thanks a lot! Commented Mar 13, 2020 at 3:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.