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.

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) - 1}{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) - 1}{z}\]

Calls phiv_dense on each of the basis vectors to obtain the answer. If A is Diagonal, instead calls the scalar phi on each diagonal element and the return values are also Diagonals

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}}();

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.

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!.