Как решить линейную систему, где оба входа редки? - PullRequest
Новая
       22

Как решить линейную систему, где оба входа редки?

6 голосов
/ 13 марта 2020

есть ли в Джулии эквивалент scipy.sparse.linalg.spsolve? Вот описание функции в 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.

Я не смог найти ее в 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> 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 Ответ

8 голосов
/ 13 марта 2020

Да, это функция \. 

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