Conjugate gradient with incomplete Cholesky preconditioner


#1

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.


#2

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.


#3

This would probably be a useful PR to IncompleteSelectedInversion.jl


#4

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.


#5

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


#6

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


#7

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)


#8

Works beautifully, thanks! I guess the PR is ready. Although, I think it does the same amount of work as before :smile: but instead of copying the result, I am copying b into y, but with less allocations.


#9

https://github.com/ettersi/IncompleteSelectedInversion.jl/pull/1