ExponentialAction

Introduction

ExponentialAction.jl is a lightweight package that implements the action of the Matrix exponential using the algorithm of Al-Mohy and Higham [AlMohyHigham2011][Expmv]. For details, see the docstring of expv.

This is the same algorithm used by SciPy's expm_multiply.

While exp(X) is only defined in LinearAlgebra for Diagonal, Symmetric{<:Real}/Hermitian, and StridedMatrix, expv can take an arbitrary matrix type.

Example

julia> using LinearAlgebra, ExponentialAction
julia> A = [1.0 2.0; -2.0 3.0];
julia> B = [1.0, 2.0];
julia> exp(A) * B2-element Vector{Float64}:
 11.445823453626826
 -2.372722541157027
julia> expv(1, A, B)2-element Vector{Float64}:
 11.44582345362683
 -2.3727225411570245
julia> exp(2A) * B2-element Vector{Float64}:
  -81.7559133436684
 -103.56648780897383
julia> expv(2, A, B)2-element Vector{Float64}:
  -81.75591334366831
 -103.56648780897383

Automatic Differentiation (AD)

Special care has been taken to ensure that expv can be differentiated using ForwardDiff.jl, ReverseDiff.jl, ChainRules.jl-compatible packages (e.g. Zygote.jl), and likely others as a result.

This has been achieved by avoiding type constraints, not mutating any arrays, and marking operations only used for control flow as being non-differentiable.

ExponentialUtilities.jl and Expokit.jl both implement an approximation to the action of the matrix exponential using Krylov subspaces. Which package is most efficient or useful depends on the choice of matrix and whether derivatives are needed. If efficiency is important, we recommend choosing from the packages by benchmarking against several of your matrices (if applicable).

Note

The packages use different default tolerances. ExponentialAction.jl is more strict. For a fair comparison, select similar tolerances.