# Pure-Julia Sparse Cholesky
**URL:** https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293
**Category:** Numerics
**Created:** [August 1, 2025, 3:08pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293 "2025-08-01T15:08:18Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 3:08pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/1 "2025-08-01T15:08:18Z")
I implemented the Cholesky decomposition of sparse positive-definite matrices in CliqueTrees.jl. The implementation is pure-Julia, except for some calls to BLAS and LAPACK kernels. Even the column-ordering algorithm is written in Julia!
```julia-repl
julia> using BenchmarkTools, LinearAlgebra, MatrixMarket, SuiteSparseMatrixCollection
julia> import CliqueTrees, LDLFactorizations, QDLDL
julia> ssmc = ssmc_db(); name = "bcsstk36"; matrix = mmread(joinpath(fetch_ssmc(ssmc[ssmc.name .== name, :]; format="MM")[1], "$(name).mtx"));
julia> @btime cholesky(matrix);
57.618 ms (55 allocations: 120.63 MiB)
julia> @btime CliqueTrees.cholesky(matrix);
67.753 ms (94 allocations: 146.47 MiB)
julia> @btime QDLDL.qdldl(matrix);
209.905 ms (89 allocations: 96.04 MiB)
julia> @btime LDLFactorizations.ldlt(matrix);
274.736 ms (56 allocations: 188.30 MiB)
```
Here is the performance without BLAS or LAPACK.
```julia-repl
julia> @btime CliqueTrees.cholesky(matrix);
97.974 ms (94 allocations: 137.08 MiB)
```
CliqueTrees.jl is faster than LDLFactorizations.jl and QDLDL.jl because it is performing a _supernodal_ decomposition, whereas they are performing _nodal_ decompositions. A supernodal decomposition partitions the input matrix into dense sub-matrices, calling BLAS Level 3 routines on the blocks.
There seems to be some interest in this sort of thing from the community. For example, [this post](https://discourse.julialang.org/t/reality-check-how-useful-exactly-are-suitesparse-and-ldlfactorizations-ldlt/79607) by @PetrKryslUCSD. I am less familiar with LDLT decompositions and quasi-definite matrices, but I could try implementing that next.
The numerical part of the algorithm is contained in [this file](https://github.com/AlgebraicJulia/CliqueTrees.jl/blob/main/src/chol_facts.jl).
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [August 1, 2025, 3:25pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/2 "2025-08-01T15:25:27Z")
I think most of the interest is in [GitHub - NextLinearAlgebra/NextLA.jl](https://github.com/NextLinearAlgebra/NextLA.jl) (espacially if you can adapt it to gpu too) and in LinearSolve.jl for those thing. Quick view on your file nothing should stop you from writing most of it in kernels ( KernelAbstractions.jl) unless they are recursion of course
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 3:35pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/3 "2025-08-01T15:35:11Z")
Hello @yolhan_mannes!
Thank you for the advice. The algorithm is entirely non-recursive. I have never used KernelAbstractions.jl before, but I will look into it!
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [August 1, 2025, 3:36pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/4 "2025-08-01T15:36:18Z")
Its ok you can also just do a PR on NextLA.jl (try adpating to their style before) and add tests and people may help after for the gpu adpatation.
Not sure they want anything to do with sparse though, LinearSolve.jl might be better than
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [August 1, 2025, 6:29pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/5 "2025-08-01T18:29:11Z")
nice! does this work for non native types? being able to do factorizations for fancy types is really nice
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 6:38pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/6 "2025-08-01T18:38:55Z")
Yes! I have a pure-Julia fallback for each BLAS and LAPACK kernel. See [here](https://github.com/AlgebraicJulia/CliqueTrees.jl/blob/728cb43d8fb7830d78f5ecc316ba1b2bf21d8c25/src/chol_fact.jl#L580-L799).
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [August 1, 2025, 6:39pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/7 "2025-08-01T18:39:44Z")
I think most of them are implemented in NextLA.jl already so you won’t have to add them to the PR
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 1, 2025, 7:13pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/8 "2025-08-01T19:13:01Z")
This is really nice.
I wonder if it can be pushed over the default `cholesky()`.
Do you have any plans for the _Incomplete Cholesky Factorization_?
There are some useful variants: Threshold based, Pattern based and Number of Non Zero Elements based. I guess the pattern one makes the most sense Graph wise.
The LDLT variant (For the Incomplete as well) will be greatly appreciated.
@yolhan_mannes , The [`NextLA.jl`](https://github.com/NextLinearAlgebra/NextLA.jl) package seems to be for Dense Linear Algebra. This algorithm is for the sparse case.
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 7:30pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/9 "2025-08-01T19:30:16Z")
Hi @RoyiAvital!
> Do you have any plans for the _Incomplete Cholesky Factorization_?
This is also something that I do not know too much about, but if there is interest then I can try! It looks like incomplete factorizations are used as pre-conditioners?
Ah, I found this: [GitHub - RoyiAvital/IncompleteCholeskyDecomposition: Implementation of the Incomplete Cholesky Decomposition with Thresholding](https://github.com/RoyiAvital/IncompleteCholeskyDecomposition)
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [August 1, 2025, 7:43pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/10 "2025-08-01T19:43:18Z")
Looks really intriguing. I wonder if you have compared performance-wise with [PetrKryslUCSD/Sparspak.jl: Direct solution of large sparse systems of linear algebraic equations in pure Julia](https://github.com/PetrKryslUCSD/Sparspak.jl) [PetrKryslUCSD/Sparspak.jl: Direct solution of large sparse systems of linear algebraic equations in pure Julia](https://github.com/PetrKryslUCSD/Sparspak.jl)?
LU is working, LDLT is implemented, but not tested (yet).
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 1, 2025, 7:43pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/11 "2025-08-01T19:43:56Z")
Indeed.
When solving _very_ large scale sparse systems doing the full Cholesky decomposition (Or other decompositions) might generate infeasible matrix (Too many elements).
Hence the need for _Incomplete Cholesky Factorization_, _Incomplete LDL Factorization_ and _Incomplete LU Factorization_.
Specifically the forms which can limit the number of elements: Pattern based and Number of Non Zero Elements based.
My code is just a redo of a Python code I found and it is basically the Threshold based method which might be nice, yet can not guarantee the number of elements of the output.
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 8:02pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/12 "2025-08-01T20:02:26Z")
Hello @PetrKryslUCSD! Am I doing this right?
```julia
import Sparspak
n = size(matrix, 2)
problem = Sparspak.Problem.Problem(n, n)
@time Sparspak.Problem.insparse!(problem, matrix)
@time Sparspak.Problem.infullrhs!(problem, 1:n)
@time solver = Sparspak.SpkSparseSpdBase._SparseSpdBase(problem)
@time Sparspak.SpkSparseSpdBase._factor!(solver)
```
The output is
```julia
0.197469 seconds (45 allocations: 137.235 MiB, 46.58% gc time)
0.000082 seconds (1 allocation: 32 bytes)
0.004005 seconds (34 allocations: 9.860 MiB)
UndefVarError: `_ldltfactor!` not defined in `Sparspak.SpkSparseSpdBase`
Suggestion: check for spelling errors or missing imports.
```
If this is just a bug then I can dev the repository and fix it.
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [August 1, 2025, 8:07pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/13 "2025-08-01T20:07:44Z")
Sorry, looks like a bug. When I said “not tested” I wasn’t kidding…
```julia-auto
solver = Sparspak.SpkSparseSpdBase._SparseSpdBase(problem)
```
This is probably not the right way. I will have a look.
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 8:17pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/14 "2025-08-01T20:17:35Z")
Probably this is right.
```julia
@time solver = Sparspak.SpkSparseSpdSolver.SparseSpdSolver(problem)
@time Sparspak.SpkSparseSpdSolver.solve!(solver)
```
But it also errors.
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 9:10pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/15 "2025-08-01T21:10:11Z")
I got the factorization part to run without erroring. See [here](https://github.com/samuelsonric/Sparspak.jl).
```julia
import Sparspak
function spk_ldl(matrix)
n = size(matrix, 2)
problem = Sparspak.Problem.Problem(n, n)
Sparspak.Problem.insparse!(problem, matrix)
Sparspak.Problem.infullrhs!(problem, 1:n)
solver = Sparspak.SpkSparseSpdSolver.SparseSpdSolver(problem)
Sparspak.SpkSparseSpdSolver.solve!(solver)
return solver
end
@btime spk_ldl(matrix);
```
The output is
```julia
2.000 s (76606554 allocations: 1.58 GiB)
```
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [August 1, 2025, 9:22pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/16 "2025-08-01T21:22:21Z")
Very cool.
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 10:14pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/17 "2025-08-01T22:14:10Z")
Just released [patch `1.9.1`](https://github.com/JuliaRegistries/General/pull/135795). We now use less memory than SuiteSparse.
```julia-repl
julia> @btime cholesky(matrix)
56.641 ms (55 allocations: 120.63 MiB)
julia> @btime CliqueTrees.cholesky(matrix)
66.369 ms (95 allocations: 109.40 MiB)
```
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [August 1, 2025, 11:12pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/18 "2025-08-01T23:12:06Z")
At the moment I don’t understand what is going on. The LDLT factorization routine has
at this point serious bugs, and does not compile. So I am at a loss to explain why your cold could run.
---
### Author:  [@samuelsonric](https://discourse.julialang.org/u/samuelsonric)
#### Post date: [August 1, 2025, 11:13pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/19 "2025-08-01T23:13:52Z")
I fixed the bugs, I think. See the most recent commit: [GitHub - samuelsonric/Sparspak.jl: Direct solution of large sparse systems of linear algebraic equations in pure Julia](https://github.com/samuelsonric/Sparspak.jl)
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [August 2, 2025, 6:22pm UTC](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293/20 "2025-08-02T18:22:04Z")
Alas, not quite. There is a bug in the factorization (division by zero).
[Next page](https://discourse.julialang.org/t/pure-julia-sparse-cholesky/131293.md?page=2)