# Conjugate gradient with incomplete Cholesky preconditioner

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, Krylov.jl and 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:

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 (recommended in IterativeSolvers docs) and the lldl function of 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.

3 Likes

Ended up doing the following:

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.

This would probably be a useful PR to IncompleteSelectedInversion.jl

1 Like

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.

Isnâ€™t this exactly what A_ldiv_B!(y, C.L,b) does? LowerTriangular seems to support this method.

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

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)

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.

4 Likes

Is there an highly optimized Incomplete Cholesky Factorization in Julia?
Something similar to MATLABâ€™s ichol()?

1 Like

Not sure how optimized this is but there is https://github.com/JuliaSmoothOptimizers/LimitedLDLFactorizations.jl which is wrapped in Preconditioners.jl to provide ichol.

1 Like

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. Yet it seems MATLAB is still faster (Also has more options).

Worth benchmarking the Julia version then perhaps attempting to accelerate it with LoopVectorization or Tullio.

1 Like

If your input matrix is posdef, it will return D posdef, effectively computing an incomplete Cholesky.

Does it have all the modes of MATLABâ€™s ichol()?

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.

1 Like

I am not sure this is accurate. Have a look at my implementation IncompleteCholeskyDecompositionThreshold Project.
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() supports this.
2. Pattern Based - IC( l )
MATLABâ€™s ichol() 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.

LDL does 1, 2 and 3. Itâ€™s designed to do 3. 2 is a special case with p = 0.

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?

See the lldl docstring. 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.