Preconditioners

Many linear solvers can be accelerated by using what is known as a preconditioner, an approximation to the matrix inverse action which is cheap to evaluate. These can improve the numerical conditioning of the solver process and in turn improve the performance. LinearSolve.jl provides an interface for the definition of preconditioners which works with the wrapped iterative solver packages.

Using Preconditioners

Mathematical Definition

A right preconditioner, $P_r$ transforms the linear system $Au = b$ into the form:

\[AP_r^{-1}(P_r u) = AP_r^{-1}y = b\]

which is solved for $y$, and then $P_r u = y$ is solved for $u$. The left preconditioner, $P_l$, transforms the linear system into the form:

\[P_l^{-1}Au = P_l^{-1}b\]

A two-sided preconditioned system is of the form:

\[P_l^{-1}A P_r^{-1} (P_r u) = P_l^{-1}b\]

Specifying Preconditioners

One way to specify preconditioners uses the Pl and Pr keyword arguments to init or solve: Pl for left and Pr for right preconditioner, respectively. By default, if no preconditioner is given, the preconditioner is assumed to be the identity $I$.

In the following, we will use a left sided diagonal (Jacobi) preconditioner.

import LinearSolve as LS
import LinearAlgebra as LA
n = 4

A = rand(n, n)
b = rand(n)

Pl = LA.Diagonal(A)

prob = LS.LinearProblem(A, b)
sol = LS.solve(prob, LS.KrylovJL_GMRES(), Pl = Pl)
sol.u

4-element Vector{Float64}:
 -0.4787743575456908
  6.566558541442251
 -2.9025985915651593
 -4.1204177656700445

Alternatively, preconditioners can be specified via the precs argument to the constructor of an iterative solver specification. This argument shall deliver a factory method mapping A and a parameter p to a tuple (Pl,Pr) consisting a left and a right preconditioner.

import LinearSolve as LS
import LinearAlgebra as LA
n = 4

A = rand(n, n)
b = rand(n)

prob = LS.LinearProblem(A, b)
sol = LS.solve(prob, LS.KrylovJL_GMRES(precs = (A, p) -> (LA.Diagonal(A), LA.I)))
sol.u

4-element Vector{Float64}:
  0.9308787037046939
 -2.6926956657921934
  0.9323856236325812
  0.04404686672421762

This approach has the advantage that the specification of the preconditioner is possible without the knowledge of a concrete matrix A. It also allows to specify the preconditioner via a callable object and to pass parameters to the constructor of the preconditioner instances. The example below also shows how to reuse the preconditioner once constructed for the subsequent solution of a modified problem.

import LinearSolve as LS
import LinearAlgebra as LA

Base.@kwdef struct WeightedDiagonalPreconBuilder
    w::Float64
end

(builder::WeightedDiagonalPreconBuilder)(A, p) = (builder.w * LA.Diagonal(A), LA.I)

n = 4
A = n * LA.I - rand(n, n)
b = rand(n)

prob = LS.LinearProblem(A, b)
cache = LS.init(prob, LS.KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w = 0.9)))
sol = LS.solve!(cache)
sol.u

B = A .+ 0.1
LS.reinit!(cache, A = B, reuse_precs = true)
sol = LS.solve!(cache, LS.KrylovJL_GMRES(precs = WeightedDiagonalPreconBuilder(w = 0.9)))
sol.u

4-element Vector{Float64}:
 0.18736114498314185
 0.2363276756330867
 0.29624201992601934
 0.3357177240477127

Preconditioner Interface

To define a new preconditioner you define a Julia type which satisfies the following interface:

• Base.eltype(::Preconditioner) (Required only for Krylov.jl)
• LinearAlgebra.ldiv!(::AbstractVector,::Preconditioner,::AbstractVector) and LinearAlgebra.ldiv!(::Preconditioner,::AbstractVector)

Curated List of Pre-Defined Preconditioners

The following preconditioners match the interface of LinearSolve.jl.

• LinearSolve.ComposePreconditioner(prec1,prec2): composes the preconditioners to apply prec1 before prec2.

• LinearSolve.InvPreconditioner(prec): inverts mul! and ldiv! in a preconditioner definition as a lazy inverse.

• LinearAlgera.Diagonal(s::Union{Number,AbstractVector}): the lazy Diagonal matrix type of Base.LinearAlgebra. Used for efficient construction of a diagonal preconditioner.

• Other Base.LinearAlgera types: all define the full Preconditioner interface.

• IncompleteLU.ilu: an implementation of the incomplete LU-factorization preconditioner. This requires A as a SparseMatrixCSC.

• Preconditioners.CholeskyPreconditioner(A, i): An incomplete Cholesky preconditioner with cut-off level i. Requires A as a AbstractMatrix and positive semi-definite.

• AlgebraicMultigrid: Implementations of the algebraic multigrid method. Must be converted to a preconditioner via AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.precmethod(A)). Requires A as a AbstractMatrix. Provides the following methods:

  • AlgebraicMultigrid.ruge_stuben(A)
  • AlgebraicMultigrid.smoothed_aggregation(A)

• AMGCLWrap.jl: algebraic multigrid and relaxation (incomplete LU, SPAI, ...) preconditioners backed by the AMGCL C++ library (CPU, multithreaded via OpenMP). Requires A as a SparseMatrixCSC or a SparseMatrixCSR. AMGPreconBuilder() and RLXPreconBuilder() plug directly into the precs interface, and AMGPrecon(A) / RLXPrecon(A) construct the preconditioners directly. On a 2-D finite-difference Laplacian with ten thousand unknowns, the AMG preconditioner takes KrylovJL_CG from 297 iterations to 10:

  using LinearSolve, AMGCLWrap
  sol = solve(prob, KrylovJL_CG(precs = AMGPreconBuilder()))

• PyAMG via LinearSolvePyAMG.jl: Implementations of the algebraic multigrid method backed by the Python PyAMG library via PythonCall.jl. The Python dependency is installed automatically via CondaPkg.jl. Provides the following solvers through the standard LinearSolve interface:

  • PyAMG() — Ruge–Stüben AMG (default)
  • PyAMG(method = :SmoothedAggregation) — smoothed-aggregation AMG
  • PyAMG(accel = "cg") — AMG preconditioned CG
  • PyAMG_RugeStuben() and PyAMG_SmoothedAggregation() shortcuts

  using LinearSolvePyAMG, LinearSolve, SparseArrays
  A = spdiagm(-1 => -ones(99), 0 => 2ones(100), 1 => -ones(99))
  b = rand(100)
  sol = solve(LinearProblem(A, b), PyAMG(accel = "cg"))

• ILUZero.ILU0Precon(A::SparseMatrixCSC{T,N}, b_type = T): An incomplete LU implementation. Requires A as a SparseMatrixCSC.

• LimitedLDLFactorizations.lldl: A limited-memory LDLᵀ factorization for symmetric matrices. Requires A as a SparseMatrixCSC. Applying F = lldl(A); F.D .= abs.(F.D) before usage as a preconditioner makes the preconditioner symmetric positive definite and thus is required for Krylov methods which are specialized for symmetric linear systems.

• RandomizedPreconditioners.NystromPreconditioner A randomized sketching method for positive semidefinite matrices A. Builds a preconditioner $P ≈ A + μ*I$ for the system $(A + μ*I)x = b$.

• HYPRE.jl A set of solvers with preconditioners which supports distributed computing via MPI. These can be written using the LinearSolve.jl interface choosing algorithms like HYPRE.ILU and HYPRE.BoomerAMG.

• KrylovPreconditioners.jl: Provides GPU-ready preconditioners via KernelAbstractions.jl. At the time of writing the package provides the following methods:

  • Incomplete Cholesky decomposition KrylovPreconditioners.kp_ic0(A)
  • Incomplete LU decomposition KrylovPreconditioners.kp_ilu0(A)
  • Block Jacobi KrylovPreconditioners.kp_block_jacobi(A)

  Reach for these rather than CUSPARSE's ilu02/ic02 directly. Those return the factors packed back into a CuSparseMatrix, which stores the factorization but has no ldiv! to apply it, so it cannot be used as a preconditioner as-is.

GPU algebraic multigrid

AMGX.jl exposes NVIDIA's algebraic multigrid. It is a solver in its own right, but on a GPU sparse matrix it is most useful as the preconditioner of a Krylov method, which is what precs is for:

using LinearSolve, CUDA, AMGX
A = CuSparseMatrixCSR(A_cpu)
alg = KrylovJL_GMRES(precs = (A, p) -> (AMGXPreconditioner(A), I))
sol = solve(LinearProblem(A, b_gpu), alg)

On a 400x400 tridiagonal system that takes GMRES from 14 iterations to 6. A has to be a CuSparseMatrixCSR. config takes a Dict{String, String} of AMGX configuration entries if the default single V-cycle is not what you want.

AMGX.jl needs a system AMGX installation, and points at it with AMGX.set_libAMGX_path.

AMGXPreconditioner(A; config = nothing)

using LinearSolve, CUDA, AMGX
alg = KrylovJL_GMRES(precs = (A, p) -> (AMGXPreconditioner(A), I))
sol = solve(LinearProblem(A_gpu_csr, b_gpu), alg)