Preconditioners

The solvers in Krylov.jl support preconditioners, i.e., transformations that modify a linear system $Ax = b$ into an equivalent form that may yield faster convergence in finite-precision arithmetic. Preconditioning can be used to reduce the condition number of the problem or cluster its eigenvalues or singular values for instance.

The design of preconditioners is highly dependent on the origin of the problem and most preconditioners need to take application-dependent information and structure into account. Specialized preconditioners generally outperform generic preconditioners such as incomplete factorizations.

The construction of a preconditioner necessitates trade-offs because we need to apply it at least once per iteration within a Krylov method. Hence, a preconditioner must be constructed such that it is cheap to apply, while also capturing the characteristics of the original system in some sense.

There exist three variants of preconditioning:

where $P_{\ell}$ and $P_r$ are square and nonsingular.

The left preconditioning preserves the error $x_k - x^{\star}$ whereas the right preconditioning keeps invariant the residual $b - A x_k$. Two-sided preconditioning is the only variant that allows to preserve the hermicity of a linear system.

In Krylov.jl, we call $P_{\ell}^{-1}$ and $P_r^{-1}$ the preconditioners and we assume that we can apply them with the operation $y \leftarrow P^{-1} * x$. It is also common to call $P_{\ell}$ and $P_r$ the preconditioners if the equivalent operation $y \leftarrow P~\backslash~x$ is available. Krylov.jl supports both approaches thanks to the argument ldiv of the Krylov solvers.

How to use preconditioners in Krylov.jl?

Square non-Hermitian linear systems

Methods concerned: CGS, BILQ, QMR, BiCGSTAB, DQGMRES, GMRES, BLOCK-GMRES, FGMRES, DIOM and FOM.

A Krylov method dedicated to non-Hermitian linear systems allows the three variants of preconditioning.

Hermitian linear systems

Methods concerned: SYMMLQ, CG, CG-LANCZOS, CG-LANCZOS-SHIFT, CR, CAR, MINRES, MINRES-QLP and MINARES.

When $A$ is Hermitian, we can only use centered preconditioning $L^{-1}AL^{-H}y = L^{-1}b$ with $x = L^{-H}y$. Centered preconditioning is a special case of two-sided preconditioning with $P_{\ell} = L = P_r^H$ that maintains hermicity. However, there is no need to specify $L$ and one may specify $P_c = LL^H$ or its inverse directly.

Linear least-squares problems

Methods concerned: CGLS, CGLS-LANCZOS-SHIFT, CRLS, LSLQ, LSQR and LSMR.

LSLQ, LSQR and LSMR also handle regularized least-squares problems.

Linear least-norm problems

Methods concerned: CGNE, CRMR, LNLQ, CRAIG and CRAIGMR.

LNLQ, CRAIG and CRAIGMR also handle penalized minimum-norm problems.

Saddle-point and symmetric quasi-definite systems

TriCG and TriMR can take advantage of the structure of Hermitian systems $Kz = d$ with the 2x2 block structure

\[ \begin{bmatrix} \tau E & \phantom{-}A \\ A^H & \nu F \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} b \\ c \end{bmatrix},\]

Generalized saddle-point and unsymmetric partitioned systems

GPMR can take advantage of the structure of general square systems $Kz = d$ with the 2x2 block structure

\[ \begin{bmatrix} \lambda M & A \\ B & \mu N \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} b \\ c \end{bmatrix},\]

Packages that provide preconditioners

• KrylovPreconditioners.jl implements block-Jacobi, IC(0) and ILU(0) preconditioners.
• ILUZero.jl is a Julia implementation of incomplete LU factorization with zero level of fill-in.
• LimitedLDLFactorizations.jl for limited-memory LDLᵀ factorization of symmetric matrices.
• AlgebraicMultigrid.jl provides two algebraic multigrid (AMG) preconditioners.
• RandomizedPreconditioners.jl uses randomized numerical linear algebra to construct approximate inverses of matrices.
• BasicLU.jl uses a sparse LU factorization to compute a maximum volume basis that can be used as a preconditioner for least-norm and least-squares problems.

Examples

using KrylovPreconditioners, Krylov

P⁻¹ = BlockJacobiPreconditioner(A)  # Block-Jacobi preconditioner
x, stats = gmres(A, b, M=P⁻¹)

using Krylov
n, m = size(A)
d = [A[i,i] ≠ 0 ? 1 / abs(A[i,i]) : 1 for i=1:n]  # Jacobi preconditioner
P⁻¹ = diagm(d)
x, stats = symmlq(A, b, M=P⁻¹)

using Krylov
n, m = size(A)
d = [1 / norm(A[:,i]) for i=1:m]  # diagonal preconditioner
P⁻¹ = diagm(d)
x, stats = minres(A, b, M=P⁻¹)

using KrylovPreconditioners, Krylov
Pℓ = ilu(A)
x, stats = gmres(A, b, M=Pℓ, ldiv=true)  # left preconditioning

using LimitedLDLFactorizations, Krylov
P = lldl(A)
P.D .= abs.(P.D)
x, stats = cg(A, b, M=P, ldiv=true)  # centered preconditioning

using ILUZero, Krylov
Pᵣ = ilu0(A)
x, stats = bicgstab(A, b, N=Pᵣ, ldiv=true)  # right preconditioning

using LDLFactorizations, Krylov

M = ldl(E)
N = ldl(F)

# [E   A] [x] = [b]
# [Aᴴ -F] [y]   [c]
x, y, stats = tricg(A, b, c, M=M, N=N, ldiv=true)

using RandomizedPreconditioners, Krylov

 = NystromSketch(A, k, r)

P = NystromPreconditioner(Â, μ)
x, stats = cg(A + μ*I, b; M=P, ldiv=true)

P⁻¹ = NystromPreconditionerInverse(Â, μ)
x, stats = cg(A + μ*I, b; M=P⁻¹)

using SuiteSparse, Krylov
import LinearAlgebra.ldiv!

M = cholesky(E)

# ldiv! is not implemented for the sparse Cholesky factorization (SuiteSparse.CHOLMOD)
ldiv!(y::Vector{T}, F::SuiteSparse.CHOLMOD.Factor{T}, x::Vector{T}) where T = (y .= F \ x)

# [E  A] [x] = [b]
# [Aᴴ 0] [y]   [c]
x, y, stats = trimr(A, b, c, M=M, sp=true, ldiv=true)

using Krylov

C = lu(M)

# [M  A] [x] = [b]
# [B  0] [y]   [c]
x, y, stats = gpmr(A, B, b, c, C=C, gsp=true, ldiv=true)

import BasicLU
using LinearOperators, Krylov

# Least-squares problem
m, n = size(A)
Aᴴ = sparse(A')
basis, B = BasicLU.maxvolbasis(Aᴴ)
opA = LinearOperator(A)
B⁻ᴴ = LinearOperator(Float64, n, n, false, false, (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'T')),
                                                  (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'N')),
                                                  (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'N')))

d, stats = lsmr(opA * B⁻ᴴ, b)  # min ‖AB⁻ᴴd - b‖₂
x = B⁻ᴴ * d                    # recover the solution of min ‖Ax - b‖₂

# Least-norm problem
m, n = size(A)
basis, B = maxvolbasis(A)
opA = LinearOperator(A)
B⁻¹ = LinearOperator(Float64, m, m, false, false, (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'N')),
                                                  (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'T')),
                                                  (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'T')))

x, y, stats = craigmr(B⁻¹ * opA, B⁻¹ * b)  # min ‖x‖₂  s.t.  B⁻¹Ax = B⁻¹b