Matrix Exponentials

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

Compute a matrix exponential, mutating A when the selected method supports in-place evaluation.

Arguments

• A: square, non-sparse matrix. Mutable inputs are overwritten by most methods; immutable inputs use ExpMethodGeneric and are returned out of place.
• method: exponential implementation, such as ExpMethodNative, ExpMethodDiagonalization, ExpMethodGeneric, or ExpMethodHigham2005.
• cache: optional workspace from alloc_mem(A, method) for repeated calls.

Returns

The matrix exponential. For mutable input and in-place methods, this is the mutated A.

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.

Examples

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 = nothing, expmethod = ExpMethodHigham2005Base())
    -> [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).

Arguments

• z::Number: scalar at which the phi functions are evaluated.
• k::Integer: highest phi-function order to return. It must be nonnegative.

Keywords

• cache: optional (k + 1) by (k + 1) matrix overwritten as workspace.
• expmethod: matrix-exponential algorithm used on the augmented matrix.

Returns

A vector whose entry j + 1 is $\varphi_j(z)$ for j = 0:k.

Examples

phi(0.5, 3)

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

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

Arguments

• A: square numeric matrix.
• k: highest requested phi-function order.

Keywords

• caches: optional PhiPadeCache for repeated dense Float64/ComplexF64 evaluations, or the legacy basis-vector workspace tuple described below.
• expmethod: matrix-exponential implementation for the legacy path.

Returns

A vector out where out[j + 1] approximates $\varphi_j(A)$ for j = 0:k. For diagonal input, every result is diagonal.

Examples

A = [0.0 1.0; -1.0 0.0]
phi(A, 2)

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

Compute matrix phi functions without allocating the output matrices.

Arguments

• out: vector of k + 1 preallocated square matrices. Entry out[j + 1] is overwritten with $\varphi_j(A)$.
• A: square numeric matrix.
• k: highest requested phi-function order.

Keywords

• caches: optional PhiPadeCache or legacy basis-vector workspace.
• expmethod: matrix-exponential implementation for the legacy path.

Returns

The mutated out.

Examples

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.

Arguments

• A: representative strided Float64 or ComplexF64 square matrix.
• p: largest phi-function order needed in subsequent phi! calls.

Returns

A PhiPadeCache sized for matrices matching A and orders through p.

Fields

• As, Apow, Nm, Dm, Dfact, rhs, Naux, Daux, tmp, pow1, and pow2: matrix work buffers for scaling, Padé evaluation, and recovery.
• linsolve: reusable LinearSolve cache for Padé denominator solves.
• absA, rvec1, rvec2, Ncoef, Dcoef, Amat, eta, alpha, tvals, and Cost: real-valued parameter-selection work buffers.
• info: result code from the most recent evaluation (0 on success).

Examples

A = [0.0 1.0; -1.0 0.0]
cache = PhiPadeCache(A, 2)
out = [similar(A) for _ in 0:2]
phi!(out, A, 2; caches = cache)

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

Matrix-exponential method using Higham's 2005 scaling-and-squaring Padé algorithm and generated evaluation kernels.

Arguments

• do_balancing: whether to balance a suitable dense matrix before the Padé evaluation. The no-argument constructor defaults to true; the matrix constructor selects balancing only for strided matrices.

Fields

• do_balancing::Bool: whether to apply matrix balancing.

Returns

An ExpMethodHigham2005 algorithm object for use with exponential!.

Examples

A = [0.0 1.0; -1.0 0.0]
method = ExpMethodHigham2005(A)
exponential!(copy(A), method)

ExpMethodHigham2005Base()

Matrix-exponential method using the Base-compatible Higham scaling-and-squaring implementation.

This is primarily the default reduced-matrix method used by the Krylov APIs. It has no fields or constructor arguments.

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.

Arguments

• k: optional Padé order. ExpMethodGeneric() uses order 13 for common floating-point inputs; ExpMethodGeneric(T) selects an order from T's precision.

Fields

• T: a Val{k} type parameter storing the selected Padé order.

Returns

An ExpMethodGeneric algorithm object for use with exponential!.

Examples

A = BigFloat[0 1; -1 0]
method = ExpMethodGeneric(BigFloat)
exponential!(A, method)

ExpMethodNative()

Matrix-exponential method that delegates to Base.exp.

Returns

An ExpMethodNative algorithm object for use with exponential!.

Examples

A = [0.0 1.0; -1.0 0.0]
exponential!(copy(A), ExpMethodNative())

ExpMethodDiagonalization(enforce_real = true)

Matrix-exponential method based on diagonalization with eigen.

Arguments

• enforce_real: when true (the default), discard small imaginary parts introduced by numerical diagonalization for real input matrices.

Fields

• enforce_real::Bool: whether real input receives a real-valued result.

Returns

An ExpMethodDiagonalization algorithm object for use with exponential!.

Examples

A = [0.0 1.0; -1.0 0.0]
exponential!(copy(A), ExpMethodDiagonalization())

alloc_mem(A, method)

Allocate reusable workspace for exponential!.

alloc_mem is the public allocation hook for matrix-exponential methods. Method implementations may specialize it together with exponential!(A, method, cache) to separate workspace allocation from repeated evaluations. Treat the returned cache as opaque and pass it back only for inputs compatible with the prototype A and the same method configuration.

Arguments

• A: prototype input that determines the workspace's shape, element type, storage type, and execution device.
• method: matrix-exponential method whose workspace should be allocated.

Returns

A method-specific cache accepted as the third argument to exponential!, or nothing when method does not require reusable workspace.

Examples

julia> using ExponentialUtilities, LinearAlgebra

julia> A = [0.0 1.0; -1.0 0.0];

julia> reference = exp(A);

julia> method = ExpMethodHigham2005();

julia> cache = ExponentialUtilities.alloc_mem(A, method);

julia> exponential!(A, method, cache) ≈ reference
true