Matrix Exponentials

E=exponential!(A,[method [cache]])

Computes the matrix exponential with the method specified in method. The contents of A are modified, allowing for fewer allocations. The method parameter specifies the implementation and implementation parameters, e.g. ExpMethodNative, ExpMethodDiagonalization, ExpMethodGeneric, ExpMethodHigham2005. Memory needed can be preallocated and provided in the parameter cache such that the memory can be recycled when calling exponential! several times. The preallocation is done with the command alloc_mem: cache=alloc_mem(A,method). A may not be sparse matrix type, since exp(A) is likely to be dense.

If no method is given, immutable matrices (e.g. StaticArrays' SMatrix) are computed out-of-place with ExpMethodGeneric and the result is returned without modifying A.

Example

julia> A = randn(50, 50);

julia> B = A * 2;

julia> method = ExpMethodHigham2005();

julia> cache = ExponentialUtilities.alloc_mem(A, method); # Main allocation done here

julia> E1 = exponential!(A, method, cache) # Very little allocation here

julia> E2 = exponential!(B, method, cache) # Very little allocation here

phi(z,k[;cache]) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the scalar phi functions for all orders up to k.

The phi functions are defined as

\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - \varphi_k(0)}{z}\]

Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156. Theorem 1).

phi(A,k[;cache]) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the matrix phi functions for all orders up to k. k >= 1.

The phi functions are defined as

\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - \varphi_k(0)}{z}\]

For Float64/ComplexF64 matrices this uses the scaling-and-recovering algorithm of Al-Mohy and Liu (arXiv:2506.01193), which computes all of phi_0(A), ..., phi_k(A) simultaneously in O(k m^3) operations. Other element types fall back to calling phiv_dense on each basis vector (O(m (m+k)^3)). If A is Diagonal, the scalar phi is instead applied to each diagonal element and the return values are also Diagonals.

phi!(out,A,k[;caches]) -> out

Non-allocating version of phi for non-diagonal matrix inputs, writing phi_j(A) into out[j+1].

For dense Float64/ComplexF64 matrices, pass a reusable PhiPadeCache as caches to make repeated evaluations of the same size and order allocation-free (the linear solve runs through an embedded LinearSolve.jl cache):

cache = PhiPadeCache(A, k)
phi!(out, A, k; caches = cache)

Numerical failure (a singular Padé denominator or non-finite result, only possible for pathological inputs such as matrices containing NaN/Inf) does not throw: the outputs are filled with NaN and, when a PhiPadeCache is used, cache.info[] is set to a nonzero return code (0 on success). This lets adaptive integrators reject the step instead of aborting.

For the legacy basis-vector algorithm, caches is instead the tuple (Vector{T}(undef, m), Matrix{T}(undef, m, k+1), Matrix{T}(undef, m+k, m+k)); supplying it forces that code path.

PhiPadeCache(A, p)

Reusable workspace for the in-place phi-function evaluator phi! on strided Float64/ComplexF64 matrices. Allocate once for a given size and p, then pass it as the caches keyword to phi! to reuse across calls. The linear solve runs through an embedded LinearSolve.jl cache (batched matrix right-hand side, GenericLUFactorization); reuse across calls is allocation-free — all buffers, including the LU pivots, live in the workspace and the LinearSolve cache.

After each phi! call, cache.info[] holds a return code: 0 means success; a positive value indicates a numerically singular Padé denominator; -1 means the result was non-finite (NaN/Inf input). Both failures only occur for pathological inputs, and no error is thrown — the outputs are filled with NaN so that adaptive integrators can detect the failed evaluation via cache.info[] != 0 (or an isfinite check) and reject the step rather than abort.

ExpMethodHigham2005(A::AbstractMatrix);
ExpMethodHigham2005(b::Bool=true);

Computes the matrix exponential using the algorithm 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" based on generated code. If a matrix is specified, balancing is determined automatically.

ExpMethodHigham2005Base()

The same as ExpMethodHigham2005 but follows Base.exp closer.

struct ExpMethodGeneric{T}
ExpMethodGeneric()=ExpMethodGeneric{Val{13}}()
ExpMethodGeneric(k::Integer)=ExpMethodGeneric{Val{k}}()
ExpMethodGeneric(::Type{T}) where T = ExpMethodGeneric{Val{pade_order_for_type(T)}}()

Generic exponential implementation of the method ExpMethodHigham2005, for any exp argument x for which the functions LinearAlgebra.opnorm, +, *, ^, and / (including addition with UniformScaling objects) are defined. The type T is used to adjust the number of terms used in the Pade approximants at compile time.

For high-precision types like BigFloat, the Padé order is automatically selected based on the precision to achieve machine-precision accuracy. You can also manually specify the order: ExpMethodGeneric(k) uses a (k,k) Padé approximant. To automatically select based on element type, use ExpMethodGeneric(T) where T is the element type.

See "The Scaling and Squaring Method for the Matrix Exponential Revisited" by Higham, Nicholas J. in 2005 for algorithm details.

ExpMethodNative()

Matrix exponential method corresponding to calling Base.exp.

ExpMethodDiagonalization(enforce_real=true)

Matrix exponential method corresponding to the diagonalization with eigen possibly by removing imaginary part introduced by the numerical approximation.

cache=alloc_mem(A,method)

Pre-allocates memory associated with matrix exponential function method and matrix A. To be used in combination with exponential!.