Frechet Derivative of the Matrix Exponential

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

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.

References

see expmCond.

Author

Michael Stadelmann (final polish by Martin Maechler).

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)))
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    2    5    9
#> [2,]    1    3    6    1
#> [3,]    1    2    7    5
#> [4,]    1    3    8    3
E <- matrix(1e-3, 4,4)
(L.AE <- expmFrechet(A, E))
#> $expm
#>           [,1]     [,2]      [,3]     [,4]
#> [1,] 156521.03 383277.4 1091220.6 632415.2
#> [2,]  97789.55 239462.4  681763.2 395111.7
#> [3,] 138097.98 338165.5  962788.7 557980.6
#> [4,] 135649.36 332170.1  945716.1 548086.0
#> 
#> $Lexpm
#>          [,1]     [,2]     [,3]     [,4]
#> [1,] 732.9072 1571.538 4155.436 2502.036
#> [2,] 472.4668 1017.524 2697.755 1622.058
#> [3,] 650.8293 1396.816 3695.518 2224.455
#> [4,] 640.5896 1375.233 3639.057 2190.267
#> 
all.equal(L.AE, expmFrechet(A, E, "block"), tolerance = 1e-14) ## TRUE
#> [1] TRUE