Is there no ldiv! for Choleski factorization of a sparse matrix?

General Usage
1

MWE:

using SparseArrays
using LinearAlgebra
a = sparse([1, 2, 3, 4], [1, 2, 3, 4], [1.0, 2.0, 3.0, 4.0], 4, 4)
af = cholesky(a)
b = rand(4)
y = zeros(4)
ldiv!(y, af, b)

Results in error:

ERROR: MethodError: no method matching ldiv!(::SparseArrays.CHOLMOD.Factor{Float64, Int64}, ::Vector{Float64})

Closest candidates are:
  ldiv!(::QRPivoted{T, var"#s994", C} where {var"#s994"<:(StridedMatrix{T} where T), C<:AbstractVector{T}}, ::AbstractVector{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ LinearAlgebra C:\Users\pkonl\VSCode_Julia_portable\assets\julia-1.10.3\share\julia\stdlib\v1.10\LinearAlgebra\src\qr.jl:610
  ldiv!(::QRPivoted, ::AbstractVector)
   @ LinearAlgebra C:\Users\pkonl\VSCode_Julia_portable\assets\julia-1.10.3\share\julia\stdlib\v1.10\LinearAlgebra\src\qr.jl:677
  ldiv!(::LinearAlgebra.TransposeFactorization{T, <:SparseArrays.UMFPACK.UmfpackLU{T}}, ::StridedVecOrMat{T}) where T<:Union{Float64, ComplexF64}
   @ SparseArrays C:\Users\pkonl\VSCode_Julia_portable\assets\julia-1.10.3\share\julia\stdlib\v1.10\SparseArrays\src\solvers\umfpack.jl:935
  ...

Stacktrace:
 [1] ldiv!(Y::Vector{Float64}, A::SparseArrays.CHOLMOD.Factor{Float64, Int64}, B::Vector{Float64})
   @ LinearAlgebra C:\Users\pkonl\VSCode_Julia_portable\assets\julia-1.10.3\share\julia\stdlib\v1.10\LinearAlgebra\src\factorization.jl:168
 [2] top-level scope
   @ REPL[12]:1
2

Not a simple answer, but see In place Cholesky solve with zero allocation (simplical vs. supernodal)

2 Likes
3

Seems like it’s just that no one has gotten around to it. The initial implementation of \ used CHOLMOD’s out-of-place API, but it looks like there is an in-place API we could use to implement ldiv!:

There’s some question of how much of the underlying API should be exposed. e.g. CHOLMOD also lets you re-use a temporary workspace vector, should we expose that?

(I think it hasn’t been a high priority since for people doing large sparse solves the allocation of the solution vector is probably not a big part of the computational cost.)

2 Likes
4

You are right, Steven. Unless one does hundreds of solves…
Personally I am looking for simply avoiding allocations for multiple right-hand sides.

It seems the post linked to above has a decent solution, doesn’t it?

5

If the solution allocations are 1% of the cost for a single solve (even without including the factorization cost), then they should still be 1% of the cost if you do hundreds of solves, no?

6

I went with in-place LU (manque the in-place Cholesky).
The allocations went from 10 GB to 200 MB, and time was reduced nearly in half.

1 Like