nimble (version 1.4.0)

expm: Matrix Exponential

Description

Compute the the matrix exponential expm(A) by scaling and squaring.

Usage

expm(A, tol = 1e-08)

Value

a matrix that is ans = exp(A).

Arguments

A

Square matrix.

tol

level of accuracy required (default = 1e-8).

Author

Paul van Dam-Bates

Details

This function follows the scaling and squaring algorithm from Sherlock (2021), except that we compute the full matrix exponential. It differs from the standard Taylor scaling and squaring algorithm reviewed by Ruiz et al (2016) and found in common texts, by doing uniformization if the matrix is a generator matrix from a continuous time Markov chain. If using the matrix exponential to create a transition probability matrix in a HMM context just once, this function may be efficient. If dimension is large, we recommend avoiding the matrix exponential and using expAv instead. Note that for computation efficiency, when the columns are non-positive, matrix uniformization is done by A* = A + rho I, where rho = max(abs(diag(A))). When the row sums of the matrix are not zero, then uniformization is not done, and the number of iterations to reach tolerance are approximated based on the a bound of the spectrum, similar to RTMB (Kristensen, 2025).

References

Sherlock, C. (2021). Direct statistical inference for finite Markov jump processes via the matrix exponential. Computational Statistics, 36(4), 2863-2887.

Ruiz, P., Sastre, J., Ibáñez, J., & Defez, E. (2016). High performance computing of the matrix exponential. Journal of computational and applied mathematics, 291, 370-379.

Kristensen K (2025). _RTMB: 'R' Bindings for 'TMB'_. R package version 1.7, commit 6bd7a16403ccb4d3fc13ff7526827540bf27b352, <https://github.com/kaskr/RTMB>.

Examples

A <- rbind(c(-1, 0.25, 0.75), c(0, -2, 2), c(0.25, 0.25, -0.5))
Lambda <- diag(c(0.25, 0.1, 0))
expm((A-Lambda)*2.5)