Matrix Exponential

Compute the exponential of a matrix.

expm(x)

Arguments

x

a matrix, typically inheriting from the dMatrix class.

Details

The exponential of a matrix is defined as the infinite Taylor series expm(A) = I + A + A^2/2! + A^3/3! + ... (although this is definitely not the way to compute it). The method for the dgeMatrix class uses Ward's diagonal Pade' approximation with three step preconditioning, a recommendation from Moler & Van Loan (1978) “Nineteen dubious ways...”.

Value

The matrix exponential of x.

References

https://en.wikipedia.org/wiki/Matrix_exponential

Cleve Moler and Charles Van Loan (2003) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review 45, 1, 3–49. doi:10.1137/S00361445024180

for historical reference mostly:
Moler, C. and Van Loan, C. (1978) Nineteen dubious ways to compute the exponential of a matrix. SIAM Review 20, 4, 801–836. doi:10.1137/1020098

Eric W. Weisstein et al. (1999) Matrix Exponential. From MathWorld, https://mathworld.wolfram.com/MatrixExponential.html

Author

This is a translation of the implementation of the corresponding Octave function contributed to the Octave project by A. Scottedward Hodel A.S.Hodel@Eng.Auburn.EDU. A bug in there has been fixed by Martin Maechler.

See also

Package expm, which provides newer (in some cases faster, more accurate) algorithms for computing the matrix exponential via its own (non-generic) function expm(). expm also implements logm(), sqrtm(), etc.

Generic function Schur.

Examples

(m1 <- Matrix(c(1,0,1,1), ncol = 2))
#> 2 x 2 Matrix of class "dtrMatrix"
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    .    1
(e1 <- expm(m1)) ; e <- exp(1)
#> 2 x 2 Matrix of class "dtrMatrix"
#>      [,1]     [,2]    
#> [1,] 2.718282 2.718282
#> [2,]        . 2.718282
stopifnot(all.equal(e1@x, c(e,0,e,e), tolerance = 1e-15))
(m2 <- Matrix(c(-49, -64, 24, 31), ncol = 2))
#> 2 x 2 Matrix of class "dgeMatrix"
#>      [,1] [,2]
#> [1,]  -49   24
#> [2,]  -64   31
(e2 <- expm(m2))
#> 2 x 2 Matrix of class "dgeMatrix"
#>            [,1]      [,2]
#> [1,] -0.7357588 0.5518191
#> [2,] -1.4715176 1.1036382
(m3 <- Matrix(cbind(0,rbind(6*diag(3),0))))# sparse!
#> 4 x 4 sparse Matrix of class "dtCMatrix"
#>             
#> [1,] . 6 . .
#> [2,] . . 6 .
#> [3,] . . . 6
#> [4,] . . . .
(e3 <- expm(m3)) # upper triangular
#> 4 x 4 Matrix of class "dtrMatrix"
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    6   18   36
#> [2,]    .    1    6   18
#> [3,]    .    .    1    6
#> [4,]    .    .    .    1