andreasnoack/ExponentialUtilities.jl

Utility functions for exponential integrators for the SciML scientific machine learning ecosystem

ExponentialUtilities

ExponentialUtilities is a package of utility functions used by the exponential integrators in OrdinaryDiffEq. It is a lightweight pure Julia package with no external dependencies, so it can also be used independently.

Matrix-phi-vector product

The main functionality of ExponentialUtilities is the computation of matrix-phi-vector products. The phi functions are defined as

ϕ_0(z) = exp(z)
ϕ_(k+1)(z) = (ϕ_k(z) - 1) / z

In exponential algorithms, products in the form of ϕ_m(tA)b is frequently encountered. Instead of computing the matrix function first and then computing the matrix-vector product, the common alternative is to construct a Krylov subspace K_m(A,b) and then approximate the matrix-phi-vector product.

expv and phiv

expv(t,A,b;kwargs) -> exp(tA)b
phiv(t,A,b,k;kwargs) -> [ϕ_0(tA)b ϕ_1(tA)b ... ϕ_k(tA)b][, errest]

For phiv, all ϕ_m(tA)b products up to order k is returned as a matrix. This is because it's more economical to produce all the results at once than individually. A second output is returned if errest=true in kwargs. The error estimate is given for the second-to-last product, using the last product as an estimator. If correct=true, then ϕ_0 through ϕ_(k-1) are updated using the last Arnoldi vector. The correction algorithm is described in [1].

You can adjust how the Krylov subspace is constructed by setting various keyword arguments. See the Arnoldi iteration section for more details.

expv_timestep and phiv_timestep

Unlike expv and phiv, the timestepping methods divide t into smaller time steps and compute the product step-by-step. By doing this in smaller chunks, the methods allow for finer error control as well as adaptation. The timestepping algorithm is described in [1], which is based upon the numerical package Expokit [2].

exp_timestep(ts,A,b;kwargs) -> U

Evaluates the matrix exponentiation-vector product u = exp(tA)b using time stepping.

phiv_timestep(ts,A,[b_0 b_1 ... b_p];kwargs) -> U

Evaluates the linear combination of phi-vector products u = ϕ_0(tA)b_0 + tϕ_1(tA)b_1 + ... + t^pϕ_p(tA)b_p using time stepping.

In both cases, ts is an array of time snapshots for u, with U[:,j] ≈ u(ts[j]). ts can also be just one value, in which case only the end result is returned and U is a vector.

Apart from keyword arguments that affect the computation of Krylov subspaces (see the Arnoldi iteration section), you can also adjust the timestepping behavior using the arguments. By setting adaptive=true, the time step and Krylov subsapce size adaptation scheme of Niesen & Wright is used and the relative tolerance of which can be set using the keyword parameter tol. The delta and gamma parameter of the adaptation scheme can also be adjusted. The tau parameter adjusts the size of the timestep (and for adaptive=true, the initial timestep). By default, it is calculated using a heuristic formula by Niesen & Wright.

Support for matrix-free operators

You can use any object as the "matrix" A as long as it implements the following linear operator interface:

• Base.eltype(A)
• Base.size(A, dim)
• LinearAlgebra.mul!(y, A, x) (for computing y = A * x in place).
• LinearAlgebra.opnorm(A, p=Inf). If this is not implemented or the default implementation can be slow, you can manually pass in the operator norm (a rough estimate is fine) using the keyword argument opnorm.
• LinearAlgebra.ishermitian(A). If this is not implemented or the default implementation can be slow, you can manually pass in the value using the keyword argument ishermitian.

Matrix-phi function phi

phi(z,k[;cache]) -> [ϕ_0(z),ϕ_1(z),...,ϕ_k(z)]

Compute ϕ function directly. z can both be a scalar or a AbstractMatrix (note that unlike the previous functions, you need to use a concrete matrix). This is used by the caching versions of the ExpRK integrators to precompute the operators.

Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used [2].

Arnoldi iteration arnoldi

arnoldi(A,b[;m,tol,opnorm,iop,cache]) -> Ks

Performs Arnoldi iterations to obtain the Krylov subspace K_m(A,b). The result is a KrylovSubspace that can be used by phiv via the alternative interface

phiv(t,Ks,k;kwargs) -> [ϕ_0(tA)b ϕ_1(tA)b ... ϕ_k(tA)b][, errest]

The reason for having this alternative interface is that we may want to compute ϕ_m(tA)b for different values of t. In this case, we can compute Ks just once (which is expensive) and follow up with several phiv calls using Ks (which is not as expensive).

For arnoldi, if A is hermitian, then the more efficient Lanczos algorithm is used instead. For cases when A is almost hermitian or when accuracy is not important, the incomplete orthogonalization procedure (IOP) can be used by setting the IOP length iop in kwargs.

For the other keyword arguments, m determines the dimension of the Krylov subspace and tol is the relative tolerance used to determine the "happy-breakdown" condition. You can also set custom operator norm in opnorm, e.g. efficient norm estimation functions instead of the default LinearAlgebra.opnorm. Only opnorm(A, Inf) needs to be defined.

_exp!

_exp(A)

A pure Julia implementation of a non-allocating matrix exponential using the Destructive matrix exponential using algorithm from Higham, 2008. Mostly generic, though the coefficients are geared towards 64-bit floating point calculations, and the use of BLAS requires a StridedMatrix.

exp_generic

exp(x, vk=Val{10}())

A pure Julia generic implementation of the exponential function using the scaling and squaring method, working on any x for which the functions LinearAlgebra.opnorm, +, *, ^, and / (including addition with UniformScaling objects) are defined. Use the argument vk to adjust the number of terms used in the Pade approximants at compile time.

Advanced features

"In-place" versions for phi, arnoldi, expv, phiv, expv_timestep and phiv_timestep are available as phi!, arnoldi!, expv!, phiv!, expv_timestep! and phiv_timestep!. You can refer to the docstrings for more information.

In addition, you may provide pre-allocated caches to the functions to further improve efficiency. In particular, dedicated cache types for expv! and phiv! are available as ExpvCache and PhivCache. Both of them can be dynamically resized, which is crucial for the adaptive phiv_timestep! method.

References

[1] Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.

[2] Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-15.

[3] Koskela, A. (2015). Approximating the matrix exponential of an advection-diffusion operator using the incomplete orthogonalization method. In Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353). Springer, Cham.

[4] HIGHAM, N. J. (2005). "THE SCALING AND SQUARING METHOD FOR THE MATRIX EXPONENTIAL REVISITED." SIAM J. MATRIX ANAL. APPL.Vol. 26, No. 4, pp. 1179–1193