KrylovKit.jl
A Julia package collecting a number of Krylov-based algorithms for linear problems, singular value and eigenvalue problems and the application of functions of linear maps or operators to vectors.
Overview
KrylovKit.jl accepts general functions or callable objects as linear maps, and general Julia objects with vector like behavior (see below) as vectors.
The high level interface of KrylovKit is provided by the following functions:
linsolve
: solve linear systemsA*x = b
eigsolve
: find a few eigenvalues and corresponding eigenvectors of an eigenvalue problemA*x = λ x
geneigsolve
: find a few eigenvalues and corresponding vectors of a generalized eigenvalue problemA*x = λ*B*x
svdsolve
: find a few singular values and corresponding left and right singular vectorsA*x = σ * y
andA'*y = σ*x
.exponentiate
: apply the exponential of a linear map to a vectorexpintegrator
: exponential integrator for a linear non-homogeneous ODE, generalization ofexponentiate
Package features and alternatives
This section could also be titled "Why did I create KrylovKit.jl"?
There are already a fair number of packages with Krylov-based or other iterative methods, such as
- IterativeSolvers.jl: part of the JuliaMath organisation, solves linear systems and least square problems, eigenvalue and singular value problems
- Krylov.jl: part of the JuliaSmoothOptimizers organisation, solves linear systems and least square problems, specific for linear operators from LinearOperators.jl.
- KrylovMethods.jl: specific for sparse matrices
- Expokit.jl: application of the matrix exponential to a vector
- ArnoldiMethod.jl: Implicitly restarted Arnoldi method for eigenvalues of a general matrix
- JacobiDavidson.jl: Jacobi-Davidson method for eigenvalues of a general matrix
- ExponentialUtilities.jl: Krylov subspace methods for matrix exponentials and
phiv
exponential integrator products. It has specialized methods for subspace caching, time stepping, and error testing which are essential for use in high order exponential integrators. - OrdinaryDiffEq.jl: contains implementations of high order exponential integrators with adaptive Krylov-subspace calculations for solving semilinear and nonlinear ODEs.
These packages have certainly inspired and influenced the development of KrylovKit.jl. However, KrylovKit.jl distinguishes itself from the previous packages in the following ways:
KrylovKit accepts general functions to represent the linear map or operator that defines the problem, without having to wrap them in a
LinearMap
orLinearOperator
type. Of course, subtypes ofAbstractMatrix
are also supported. If the linear map (always the first argument) is a subtype ofAbstractMatrix
, matrix vector multiplication is used, otherwise it is applied as a function call.KrylovKit does not assume that the vectors involved in the problem are actual subtypes of
AbstractVector
. Any Julia object that behaves as a vector is supported, so in particular higher-dimensional arrays or any custom user type that supports the following functions (withv
andw
two instances of this type andα, β
scalars (i.e.Number
)):Base.:*(α, v)
: multiplyv
with a scalarα
, which can be of a different scalar type; in particular this method is used to create vectors similar tov
but with a different type of underlying scalars.Base.similar(v)
: a way to construct vectors which are exactly similar tov
LinearAlgebra.mul!(w, v, α)
: out of place scalar multiplication; multiply vectorv
with scalarα
and store the result inw
LinearAlgebra.rmul!(v, α)
: in-place scalar multiplication ofv
withα
; in particular withα = false
,v
is the corresponding zero vectorLinearAlgebra.axpy!(α, v, w)
: store inw
the result ofα*v + w
LinearAlgebra.axpby!(α, v, β, w)
: store inw
the result ofα*v + β*w
LinearAlgebra.dot(v,w)
: compute the inner product of two vectorsLinearAlgebra.norm(v)
: compute the 2-norm of a vector
Algorithms in KrylovKit.jl are tested against such a minimal implementation (named
MinimalVec
) in the test suite. This type is only defined in the tests. However, KrylovKit provides two types implementing this interface and slightly more, to make them behave more likeAbstractArrays
(e.g. alsoBase.:+
etc), which can facilitate certain applications:RecursiveVec
can be used for grouping a set of vectors into a single vector like structure (can be used recursively). This is more robust than trying to use nestedVector{<:Vector}
types.InnerProductVec
can be used to redefine the inner product (i.e.dot
) and corresponding norm (norm
) of an already existing vector like object. The latter should help with implementing certain type of preconditioners.
Current functionality
The following algorithms are currently implemented
linsolve
:CG
,GMRES
eigsolve
: a Krylov-Schur algorithm (i.e. with tick restarts) for extremal eigenvalues of normal (i.e. not generalized) eigenvalue problems, corresponding toLanczos
for real symmetric or complex hermitian linear maps, and toArnoldi
for general linear maps.geneigsolve
: an customized implementation of the inverse-free algorithm of Golub and Ye for symmetric / hermitian generalized eigenvalue problems with positive definite matrixB
in the right hand side of the generalized eigenvalue problem $A v = B v λ$. The Matlab implementation was described by Money and Ye and is known asEIGIFP
; in particular it extends the Krylov subspace with a vector corresponding to the step between the current and previous estimate, analogous to the locally optimal preconditioned conjugate gradient method (LOPCG). In particular, with Krylov dimension 2, it becomes equivalent to the latter.svdsolve
: finding largest singular values based on Golub-Kahan-Lanczos bidiagonalization (seeGKL
)exponentiate
: aLanczos
based algorithm for the action of the exponential of a real symmetric or complex hermitian linear map.expintegrator
: exponential integrator for a linear non-homogeneous ODE, computes a linear combination of theϕⱼ
functions which generalizeϕ₀(z) = exp(z)
.
Future functionality?
Here follows a wish list / to-do list for the future. Any help is welcomed and appreciated.
- More algorithms, including biorthogonal methods:
- for
linsolve
: MINRES, BiCG, BiCGStab(l), IDR(s), ... - for
eigsolve
: BiLanczos, Jacobi-Davidson JDQR/JDQZ, subspace iteration (?), ... - for
geneigsolve
: trace minimization, ...
- for
- Support both in-place / mutating and out-of-place functions as linear maps
- Reuse memory for storing vectors when restarting algorithms (related to previous)
- Support non-BLAS scalar types using GeneralLinearAlgebra.jl and GeneralSchur.jl
- Least square problems
- Nonlinear eigenvalue problems
- Preconditioners
- Refined Ritz vectors, Harmonic Ritz values and vectors
- Block versions of the algorithms
- More relevant matrix functions
Partially done:
- Improved efficiency for the specific case where
x
isVector
(i.e. BLAS level 2 operations): any vectorv::AbstractArray
which hasIndexStyle(v) == IndexLinear()
now benefits from a multithreaded (useexport JULIA_NUM_THREADS = x
withx
the number of threads you want to use) implementation that resembles BLAS level 2 style for the vector operations, providedClassicalGramSchmidt()
,ClassicalGramSchmidt2()
orClassicalGramSchmidtIR()
is chosen as orthogonalization routine.