Menu
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
630
The Sparse Solvers library in the Accelerate framework handles the solution of systems of equations where the coefficient matrix is sparse. That is, most of the entries in the matrix are zero. The Sparse Solvers library provides a sparse counterpart to the dense factorizations and linear solvers that LAPACK provides.
A system of linear equations is called sparse if only relatively small number of its matrix elements are nonzero. It is wasteful to use general methods of linear algebra for such problems because when we use all the elements of matrix (zero and nonzero) we perform operations.
Sparse linear algebra is at the heart of a most partial differential equation solvers and hence they are extremely common in the computational sciences. Finance problems, structural mechanics, data mining, operations research … the list of problems based on sparse linear algebra is extensive.
The solution of the sparse linear system is usually the most computationally demanding of the three steps. Solution methods include direct factorization and preconditioned iterative methods, methods which can vary dramatically in required storage and computational cost for different problems.
The linear systems in AEM can reach millions of equations and unknowns, which makes solving them efficiently a nontrivial issue. The methods for solving sparse linear systems are generally split in two categories: direct solvers and iterative solvers.
When solving linear systems, you have two methods at your disposal, and which one you choose depends on the problem: If the coefficient of any variable is 1, which means you can easily solve for it in terms of the other variable, then substitution is a very good bet.
When two or more linear equations are grouped together, they form a system of linear equations. In this section, we will focus our work on systems of two linear equations in two unknowns. We will solve larger systems of equations later in this chapter. An example of a system of two linear equations is shown below.
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With