expmFrechet {expm}

Frechet Derivative of the Matrix Exponential

Description

Compute the Frechet (actually ‘Fréchet’) derivative of the matrix exponential operator.

Usage

expmFrechet(A, E, method = c("SPS", "blockEnlarge"), expm = TRUE)

Arguments

A

square matrix (n \times n).

E

the “small Error” matrix, used in L(A,E) = f(A + E, A)

method

string specifying the method / algorithm; the default "SPS" is “Scaling + Pade + Squaring” as in the algorithm 6.4 below; otherwise see the ‘Details’ section.

expm

logical indicating if the matrix exponential itself, which is computed anyway, should be returned as well.

Details

Calculation of e^A and the Exponential Frechet-Derivative L(A,E).

When method = "SPS" (by default), the with the Scaling - Padé - Squaring Method is used, in an R-Implementation of Al-Mohy and Higham (2009)'s Algorithm 6.4.

Step 1:

Scaling (of A and E)

Step 2:

Padé-Approximation of e^A and L(A,E)

Step 3:

Squaring (reversing step 1)

method = "blockEnlarge" uses the matrix identity of

f([A E ; 0 A ]) = [f(A) Df(A); 0 f(A)]

for the 2n \times 2n block matrices where f(A) := expm(A) and Df(A) := L(A,E). Note that "blockEnlarge" is much simpler to implement but slower (CPU time is doubled for n = 100).

Value

a list with components

expm

if expm is true, the matrix exponential (n \times n matrix).

Lexpm

the Exponential-Frechet-Derivative L(A,E), a matrix of the same dimension.

Author(s)

Michael Stadelmann (final polish by Martin Maechler).

References

see expmCond.

See Also

expm.Higham08 for the matrix exponential. expmCond for exponential condition number computations which are based on expmFrechet.

Examples

(A <- cbind(1, 2:3, 5:8, c(9,1,5,3)))
E <- matrix(1e-3, 4,4)
(L.AE <- expmFrechet(A, E))
all.equal(L.AE, expmFrechet(A, E, "block"), tolerance = 1e-14) ## TRUE