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

**URL:** <https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806>\
**Category:** General Usage\
**Created:** [May 27, 2024, 11:34pm UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806 "2024-05-27T23:34:36Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [May 27, 2024, 11:34pm UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806/1 "2024-05-27T23:34:36Z")

</div>

MWE:

```julia
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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [May 27, 2024, 11:49pm UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806/2 "2024-05-27T23:49:07Z")

</div>

Not a simple answer, but see [In place Cholesky solve with zero allocation (simplical vs. supernodal)](https://discourse.julialang.org/t/in-place-cholesky-solve-with-zero-allocation-simplical-vs-supernodal/41115)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2024, 1:47am UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806/3 "2024-05-28T01:47:16Z")

</div>

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!`:

> <https://github.com/JuliaSparse/SparseArrays.jl/issues/319>
>
> It seems like this should be supported:
> \`\`\`
> julia\> using SparseArrays, LinearA…lgebra
> 
> julia\> A = sprand(10,10,0.1); A = A'A + I;
> 
> julia\> F = cholesky(A)
> SuiteSparse.CHOLMOD.Factor{Float64}
> type: LLt
> method: simplicial
> maxnnz: 11
> nnz: 11
> success: true
> 
> julia\> b = rand(10); x = similar(b);
> 
> julia\> ldiv!(x, F, b)
> ERROR: MethodError: no method matching ldiv!(::SuiteSparse.CHOLMOD.Factor{Float64}, ::Vector{Float64})
> \`\`\`
> 
> Similar to #242, except in this case CHOLMOD \*does\* seem to provide an in-place API (whereas currently \[we are calling\](72827cd31f0aa346cfc54158f9c6db1738ecd9e4) the out-of-place \`cholmod\_l\_spsolve\`):
> 
> \<img width="857" alt="image" src="https://user-images.githubusercontent.com/2913679/210141219-6b234456-7501-40da-b5d2-617cbbad4c79.png"\>
> 
> cc @ChrisRackauckas, since this just came up in https://github.com/SciML/LinearSolve.jl/pull/246

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

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [May 28, 2024, 1:55am UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806/4 "2024-05-28T01:55:43Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2024, 1:59am UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806/5 "2024-05-28T01:59:21Z")

</div>

> [@PetrKryslUCSD](#):
>
> You are right, Steven. Unless one does hundreds of solves…

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?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [May 28, 2024, 2:27am UTC](https://discourse.julialang.org/t/is-there-no-ldiv-for-choleski-factorization-of-a-sparse-matrix/114806/6 "2024-05-28T02:27:00Z")

</div>

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.