julia - 如何解决两个输入都稀疏的线性系统？

Question

是否有任何等效于 scipy.sparse.linalg.spsolve在 Julia ？这是 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.

我在 Julia 的 LinearAlgebra 中找不到这个和 SparseArrays .有什么我想念的或任何其他选择吗？

谢谢

编辑

例如:

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

在 Julia ，与 \运算符(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

Answer 1

是的，是 \功能。

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

编辑以回答问题编辑:

如果你想要 v此处改为稀疏矩阵 B ，然后我们可以继续使用 QR B的分解(请注意， B 真正稀疏的情况很少见:

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

现在:

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

我们可以测试以确保我们做对了:

julia> mysolve(A, B) ≈ A \ collect(B)
true