GitHub

Functions of matrices and linear operators

Applying a function of a matrix or linear operator to a given vector can in some cases also be computed using Krylov methods. One example is the inverse function, which exactly corresponds to what linsolve computes: $A^{-1} * b$. There are other functions $f$ for which $f(A) * b$ can be computed using Krylov techniques, i.e. where $f(A) * b$ can be well approximated in the Krylov subspace spanned by ${b, A * b, A^2 * b, ...}$.

Currently, the only family of functions of a linear map for which such a method is available are the ϕⱼ(z) functions which generalize the exponential function ϕ₀(z) = exp(z) and arise in the context of linear non-homogeneous ODEs. The corresponding Krylov method for computing is an exponential integrator, and is thus available under the name expintegrator. For a linear homogeneous ODE, the solution is a pure exponential, and the special wrapper exponentiate is available:

KrylovKit.exponentiate — Function
function exponentiate(A, t::Number, x; kwargs...)
function exponentiate(A, t::Number, x, algorithm)

Compute $y = exp(t*A) x$, where A is a general linear map, i.e. a AbstractMatrix or just a general function or callable object and x is of any Julia type with vector like behavior.

Arguments:

The linear map A can be an AbstractMatrix (dense or sparse) or a general function or callable object that implements the action of the linear map on a vector. If A is an AbstractMatrix, x is expected to be an AbstractVector, otherwise x can be of any type that behaves as a vector and supports the required methods (see KrylovKit docs).

The time parameter t can be real or complex, and it is better to choose t e.g. imaginary and A hermitian, then to absorb the imaginary unit in an antihermitian A. For the former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to the exponential action of the linear map is obtained. The argument x can be of any type and should be in the domain of A.

Return values:

The return value is always of the form y, info = exponentiate(...) with

  • y: the result of the computation, i.e. y = exp(t*A)*x

  • info: an object of type [ConvergenceInfo], which has the following fields

    • info.converged::Int: 0 or 1 if the solution y was approximated up to the requested tolerance tol.
    • info.residual::Nothing: value nothing, there is no concept of a residual in this case
    • info.normres::Real: a (rough) estimate of the error between the approximate and exact solution
    • info.numops::Int: number of times the linear map was applied, i.e. number of times f was called, or a vector was multiplied with A
    • info.numiter::Int: number of times the Krylov subspace was restarted (see below)
Check for convergence

By default (i.e. if verbosity = 0, see below), no warning is printed if the solution was not found with the requested precision, so be sure to check info.converged == 1.

Keyword arguments:

Keyword arguments and their default values are given by:

  • verbosity::Int = 0: verbosity level, i.e. 0 (no messages), 1 (single message at the end), 2 (information after every iteration), 3 (information per Krylov step)
  • krylovdim = 30: the maximum dimension of the Krylov subspace that will be constructed. Note that the dimension of the vector space is not known or checked, e.g. x₀ should not necessarily support the Base.length function. If you know the actual problem dimension is smaller than the default value, it is useful to reduce the value of krylovdim, though in principle this should be detected.
  • tol = 1e-12: the requested accuracy per unit time, i.e. if you want a certain precision ϵ on the final result, set tol = ϵ/abs(t). If you work in e.g. single precision (Float32), you should definitely change the default value.
  • maxiter::Int = 100: the number of times the Krylov subspace can be rebuilt; see below for further details on the algorithms.
  • issymmetric: if the linear map is symmetric, only meaningful if T<:Real
  • ishermitian: if the linear map is hermitian The default value for the last two depends on the method. If an AbstractMatrix is used, issymmetric and ishermitian are checked for that matrix, otherwise the default values are issymmetric = false and ishermitian = T <: Real && issymmetric.
  • eager::Bool = false: if true, eagerly try to compute the result after every expansion of the Krylov subspace to test for convergence, otherwise wait until the Krylov subspace as dimension krylovdim. This can result in a faster return, for example if the total time for the evolution is quite small, but also has some overhead, as more computations are performed after every expansion step.

Algorithm

This is actually a simple wrapper over more general method expintegrator for for integrating a linear non-homogeneous ODE.

source
KrylovKit.expintegrator — Function
function expintegrator(A, t::Number, u₀, u₁, …; kwargs...)
function expintegrator(A, t::Number, (u₀, u₁, …); kwargs...)
function expintegrator(A, t::Number, (u₀, u₁, …), algorithm)

Compute $y = ϕ₀(t*A)*u₀ + t*ϕ₁(t*A)*u₁ + t^2*ϕ₂(t*A)*u₂ + …$, where A is a general linear map, i.e. a AbstractMatrix or just a general function or callable object and u₀, u₁ are of any Julia type with vector like behavior. Here, $ϕ₀(z) = exp(z)$ and $ϕⱼ₊₁ = (ϕⱼ(z) - 1/j!)/z$. In particular, $y = x(t)$ represents the solution of the ODE $ẋ(t) = A*x(t) + ∑ⱼ t^j/j! uⱼ₊₁$ with $x(0) = u₀$.

Note

When there are only input vectors u₀ and u₁, t can equal Inf, in which the algorithm tries to evolve all the way to the fixed point y = - A \ u₁ + P₀ u₀ with P₀ the projector onto the eigenspace of eigenvalue zero (if any) of A. If A has any eigenvalues with real part larger than zero, however, the solution to the ODE will diverge, i.e. the fixed point is not stable.

Warning

The returned solution might be the solution of the ODE integrated up to a smaller time $t̃ = sign(t) * |t̃|$ with $|t̃| < |t|$, when the required precision could not be attained. Always check info.converged > 0 or info.residual == 0 (see below).

Arguments:

The linear map A can be an AbstractMatrix (dense or sparse) or a general function or callable object that implements the action of the linear map on a vector. If A is an AbstractMatrix, x is expected to be an AbstractVector, otherwise x can be of any type that behaves as a vector and supports the required methods (see KrylovKit docs).

The time parameter t can be real or complex, and it is better to choose t e.g. imaginary and A hermitian, then to absorb the imaginary unit in an antihermitian A. For the former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to the exponential action of the linear map is obtained. The arguments u₀, u₁, … can be of any type and should be in the domain of A.

Return values:

The return value is always of the form y, info = expintegrator(...) with

  • y: the result of the computation, i.e. $y = ϕ₀(t̃*A)*u₀ + t̃*ϕ₁(t̃*A)*u₁ + t̃^2*ϕ₂(t̃*A)*u₂ + …$ with $t̃ = sign(t) * |t̃|$ with $|t̃| <= |t|$, such that the accumulated error in y per unit time is at most equal to the keyword argument tol

  • info: an object of type [ConvergenceInfo], which has the following fields

    • info.converged::Int: 0 or 1 if the solution y was evolved all the way up to the requested time t.
    • info.residual: there is no residual in the conventional sense, however, this value equals the residual time t - t̃, i.e. it is zero if info.converged == 1
    • info.normres::Real: a (rough) estimate of the total error accumulated in the solution, should be smaller than tol * |t̃|
    • info.numops::Int: number of times the linear map was applied, i.e. number of times f was called, or a vector was multiplied with A
    • info.numiter::Int: number of times the Krylov subspace was restarted (see below)

Keyword arguments:

Keyword arguments and their default values are given by:

  • verbosity::Int = 0: verbosity level, i.e. 0 (no messages), 1 (single message at the end), 2 (information after every iteration), 3 (information per Krylov step)
  • krylovdim = 30: the maximum dimension of the Krylov subspace that will be constructed. Note that the dimension of the vector space is not known or checked, e.g. x₀ should not necessarily support the Base.length function. If you know the actual problem dimension is smaller than the default value, it is useful to reduce the value of krylovdim, though in principle this should be detected.
  • tol = 1e-12: the requested accuracy per unit time, i.e. if you want a certain precision ϵ on the final result, set tol = ϵ/abs(t). If you work in e.g. single precision (Float32), you should definitely change the default value.
  • maxiter::Int = 100: the number of times the Krylov subspace can be rebuilt; see below for further details on the algorithms.
  • issymmetric: if the linear map is symmetric, only meaningful if T<:Real
  • ishermitian: if the linear map is hermitian The default value for the last two depends on the method. If an AbstractMatrix is used, issymmetric and ishermitian are checked for that matrix, otherwise the default values are issymmetric = false and ishermitian = T <: Real && issymmetric.
  • eager::Bool = false: if true, eagerly try to compute the result after every expansion of the Krylov subspace to test for convergence, otherwise wait until the Krylov subspace as dimension krylovdim. This can result in a faster return, for example if the total time for the evolution is quite small, but also has some overhead, as more computations are performed after every expansion step.

Algorithm

The last method, without keyword arguments and the different vectors u₀, u₁, … in a tuple, is the one that is finally called, and can also be used directly. Here, one specifies the algorithm explicitly as either Lanczos, for real symmetric or complex hermitian linear maps, or Arnoldi, for general linear maps. Note that these names refer to the process for building the Krylov subspace, and that one can still use complex time steps in combination with e.g. a real symmetric map.

source