# Conjugate gradient with incomplete Cholesky preconditioner
**URL:** https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809
**Category:** Numerics
**Created:** [December 17, 2017, 2:10am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809 "2017-12-17T02:10:18Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 17, 2017, 2:10am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/1 "2017-12-17T02:10:18Z")
I have been trying my luck with using the conjugate gradient method to solve a sparse symmetric positive definite system of equations. I found [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl), [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) and [KrylovMethods.jl](https://github.com/lruthotto/KrylovMethods.jl), of which the first seems to be the most active and well documented. IterativeSolvers.jl also supports supplying a preconditioner assuming that it respects the following interface:
```julia
These preconditioners should support the operations
- `A_ldiv_B!(y, P, x)` computes `P \ x` in-place of `y`;
- `A_ldiv_B!(P, x)` computes `P \ x` in-place of `x`;
- and `P \ x`
```
The docstring of `A_ldiv_B!` highlights that `P` must be a factorization object, not a matrix. All of this is great.
Being the most popular general-purpose preconditioner (AFAIK), I thought of trying an incomplete Cholesky factor. A simple search got me to the `cldlt` function of [IncompleteSelectedInversion.jl](https://github.com/ettersi/IncompleteSelectedInversion.jl) (recommended in IterativeSolvers docs) and the `lldl` function of [LLDL.jl](https://github.com/JuliaSmoothOptimizers/LLDL.jl). These implement an incomplete LDL’ factorization instead, but I am not complaining! The former package returns a sparse lower triangular matrix as a `SparseMatrixCSC{Float64,Int64}` with the strict lower triangle hosting the strict lower triangle of `L`, and the diagonal is the diagonal of `D`. The latter returns a tuple of:
1. A sparse strictly lower triangular matrix as a `SparseMatrixCSC{Float64,Int64}` hosting the strict lower triangle of `L`,
2. A vector of the diagonal elements of `D`, and
3. The diagonal shift value used (not important for me).
Unfortunately none of these outputs can be just plugged in the `A_ldiv_B!` function because they are not factors, so I cannot use them as preconditioners in the `cg` function of IterativeSolvers.
Now, I can apply the preconditioner myself in some linear operator and pass that instead to `cg`, then post-process the output, but I think such an important procedure should be more streamlined than this. So I guess my question is how to create my own `Factorization` subtype instance using the information provided by the incomplete LDL’ factorization functions? Calling the constructor `Base.SparseArrays.CHOLMOD.Factor{Float64}` naively doesn’t just work, and reading a bit of cholmod.jl didn’t really help.
Edit: I could obviously subtype `Factorization` and implement the above functions myself in a PR or something. But I am asking if such a thing already exists.
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 17, 2017, 9:46am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/2 "2017-12-17T09:46:28Z")
Ended up doing the following:
```julia
module IncompleteCholeskyPreconditioner
using IncompleteSelectedInversion
import Base.LinAlg: A_ldiv_B!, \
struct CholeskyPreconditioner{T, S <: AbstractSparseMatrix{T}}
L::LowerTriangular{T, S}
end
function CholeskyPreconditioner(A, c)
L = cldlt(A,c)
@inbounds for j in 1:size(L, 2)
d = sqrt(L[j,j])
L[j,j] = d
for i in Base.Iterators.drop(nzrange(L,j), 1)
L.nzval[i] *= d
end
end
return CholeskyPreconditioner(LowerTriangular(L))
end
function A_ldiv_B!(y::AbstractVector{T}, C::CholeskyPreconditioner{T, S}, b::AbstractVector{T}) where {T, S}
y .= C.L \ b
return y
end
(\)(C::CholeskyPreconditioner{T, S}, b::AbstractVector{T}) where {T, S <: AbstractSparseMatrix{T}} = C.L \ b
export CholeskyPreconditioner
end
```
If you think there is a better package or way to do this, a suggestion is more than welcome.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 17, 2017, 10:48am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/3 "2017-12-17T10:48:27Z")
This would probably be a useful PR to IncompleteSelectedInversion.jl
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 17, 2017, 10:55am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/4 "2017-12-17T10:55:37Z")
Yes, I will submit it as a PR once I am convinced there is no other obvious more efficient way of doing it. In particular, I am concerned with the line `y .= C.L \ b`. The function `\` allocates internally for the output vector and then it is copied to `y`. I am not sure if there is a version that takes `y` as input and fills it up with the result straight away.
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [December 17, 2017, 11:12am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/5 "2017-12-17T11:12:23Z")
Isn’t this exactly what A\_ldiv\_B!(y, C.L,b) does? LowerTriangular seems to support this method.
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 17, 2017, 11:15am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/6 "2017-12-17T11:15:45Z")
Seems not! Only tried it in v0.6.1.
`MethodError: no method matching A_ldiv_B!(::Array{Float64,1}, ::LowerTriangular{Float64,SparseMatrixCSC{Float64,Int64}}, ::Array{Float64,1})`
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [December 17, 2017, 11:23am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/7 "2017-12-17T11:23:26Z")
Indeed! I’m not sure whether master needs a patch for that, but as a workaround you can copy b to y and use the two-arguments version
A\_ldiv\_B!(L::LowerTriangular{T,\<:SparseMatrixCSC{T}}, B::StridedVecOrMat) where {T} = fwdTriSolve!(L.data, B)
(base/sparse/linalg.jl)
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 17, 2017, 11:33am UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/8 "2017-12-17T11:33:56Z")
Works beautifully, thanks! I guess the PR is ready. Although, I think it does the same amount of work as before 😄 but instead of copying the result, I am copying `b` into `y`, but with less allocations.
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [December 17, 2017, 12:17pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/9 "2017-12-17T12:17:02Z")
[https://github.com/ettersi/IncompleteSelectedInversion.jl/pull/1](https://github.com/ettersi/IncompleteSelectedInversion.jl/pull/1)
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 19, 2020, 6:33pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/10 "2020-09-19T18:33:12Z")
Is there an highly optimized Incomplete Cholesky Factorization in Julia?
Something similar to MATLAB’s [`ichol()`](https://www.mathworks.com/help/matlab/ref/ichol.html)?
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [September 19, 2020, 6:35pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/11 "2020-09-19T18:35:52Z")
Not sure how optimized this is but there is [https://github.com/JuliaSmoothOptimizers/LimitedLDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LimitedLDLFactorizations.jl) which is wrapped in Preconditioners.jl to provide ichol.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 19, 2020, 6:46pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/12 "2020-09-19T18:46:49Z")
This is my question as well. Since it only assumes Symmetric Matrix I wonder if one can also assume Positive Definite Matrix can lead for farther optimizations.
I tried myself by implementing a Python Code in `C` in my [`IncompleteCholeskyDecompositionThreshold` Project](https://github.com/RoyiAvital/IncompleteCholeskyDecompositionThreshold). Yet it seems MATLAB is still faster (Also has more options).
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [September 19, 2020, 6:55pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/13 "2020-09-19T18:55:09Z")
Worth benchmarking the Julia version then perhaps attempting to accelerate it with LoopVectorization or Tullio.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [August 27, 2021, 1:40pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/14 "2021-08-27T13:40:44Z")
If your input matrix is posdef, it will return D posdef, effectively computing an incomplete Cholesky.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 27, 2021, 1:48pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/15 "2021-08-27T13:48:56Z")
Does it have all the modes of [MATLAB’s `ichol()`](https://www.mathworks.com/help/matlab/ref/ichol.html)?
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [August 27, 2021, 2:57pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/16 "2021-08-27T14:57:40Z")
Matlab’s ichol is exclusively drop-tolerance-based. LimitedLDLFactorizations is first and foremost a limited-memory factorization, i.e., it lets you predefine how much fill-in you’re willing to tolerate. In addition to that, I added a simple drop tolerance long ago. Please open issues (or better yet, pull requests) for feature requests.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 6, 2021, 12:39pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/17 "2021-09-06T12:39:03Z")
I am not sure this is accurate. Have a look at my implementation [`IncompleteCholeskyDecompositionThreshold` Project](https://github.com/RoyiAvital/IncompleteCholeskyDecomposition).
From what I understand we can have 3 types of incomplete (At least those are 3 options):
1. Threshold Based - IC( \tau ) .
MATLAB’s [`ichol()`](https://www.mathworks.com/help/matlab/ref/ichol.html) supports this.
2. Pattern Based - IC( l )
MATLAB’s [`ichol()`](https://www.mathworks.com/help/matlab/ref/ichol.html) supports this with `l = 0`, namely only where non zero elements exists. This basically also guarantees upper memory boundary.
3. Number of Non Zero Elements - IC( p )
MATLAB doesn’t support this. This allows exact memory allocation control.
So what I see is that MATLAB has 2 out of 3 and has a mode to control the amount of memory.
It also has a mode for `Modified Cholesky` which compensate the diagonal components for the inaccurate decomposition.
If I get it correctly, `LimitedLDLFactorizations.jl` is doing something like (2), hence the memory allocation is a multiplication of the number of non zeros of the input matrix. Am I right?
I can see there is also support for Thresholding so if one allocates enough memory then this become (1) as well.
So we’re missing (3) and the `Modified ICT`.
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [September 6, 2021, 12:52pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/18 "2021-09-06T12:52:55Z")
LDL does 1, 2 and 3. It’s designed to do 3. 2 is a special case with p = 0.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [September 6, 2021, 1:02pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/19 "2021-09-06T13:02:54Z")
OK, So If it does (3) why can’t I pre allocate memory and hand it to the function?
What about the diagonal compensation? I saw something like that on code. Could you clarify it?
---
### Author:  [@dpo](https://discourse.julialang.org/u/dpo)
#### Post date: [September 6, 2021, 1:46pm UTC](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809/20 "2021-09-06T13:46:49Z")
See [the lldl docstring](https://juliasmoothoptimizers.github.io/LimitedLDLFactorizations.jl/dev/reference/#LimitedLDLFactorizations.lldl-Union%7BTuple%7BSparseArrays.SparseMatrixCSC%7BTv,%20Ti%7D%7D,%20Tuple%7BTi%7D,%20Tuple%7BTv%7D%7D%20where%20%7BTv%3C:Number,%20Ti%3C:Integer%7D). You control the amount of memory allocated for the factors via the integer p. p = 0 means that no fill-in is allowed in the factors, i.e., at most nnz(A) elements are reserved for them. p \> 0 means that nnz(A) + p \* n elements are reserved.
There’s no option to pass in your own memory, but feel free to open an issue (or better yet, a pull request) if that’s a crucial feature for you.
[Next page](https://discourse.julialang.org/t/conjugate-gradient-with-incomplete-cholesky-preconditioner/7809.md?page=2)