expmFrechet: 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)

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.

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.

Author

Michael Stadelmann (final polish by Martin Maechler).

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$).

References

see expmCond.

Examples

Run this code

(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

Run the code above in your browser using DataLab