# expm.Higham08: Matrix Exponential [Higham 2008] In expm: Matrix Exponential, Log, 'etc'

## Description

Calculation of matrix exponential e^A with the ‘Scaling & Squaring’ method with balancing.

Implementation of Higham's Algorithm from his book (see references), Chapter 10, Algorithm 10.20.

The balancing option is an extra from Michael Stadelmann's Masters thesis.

## Usage

 `1` ```expm.Higham08(A, balancing = TRUE) ```

## Arguments

 `A` square matrix, may be a `"sparseMatrix"`, currently only if `balancing` is false. `balancing` logical indicating if balancing should happen (before and after scaling and squaring).

## Details

The algorithm comprises the following steps

1. 0.Balancing

2. 1.Scaling

4. 3.Squaring

5. 4.Reverse Balancing

## Value

a matrix of the same dimension as `A`, the matrix exponential of `A`.

## Author(s)

Michael Stadelmann (final polish by Martin Maechler).

## References

Higham, N.~J. (2008). Functions of Matrices: Theory and Computation; Society for Industrial and Applied Mathematics, Philadelphia, PA, USA.

Michael Stadelmann (2009). Matrixfunktionen; Analyse und Implementierung. [in German] Master's thesis and Research Report 2009-12, SAM, ETH Zurich; http://www.sam.math.ethz.ch/reports/2009, or ftp://ftp.sam.math.ethz.ch/pub/sam-reports/reports/reports2009/2009-12.pdf.

For now, the other algorithms `expm`. This will change there will be one function with optional arguments to chose the method !.

`expmCond`, to compute the exponential-condition number.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79``` ```## The *same* examples as in ../expm.Rd {FIXME} -- x <- matrix(c(-49, -64, 24, 31), 2, 2) expm.Higham08(x) ## ---------------------------- ## Test case 1 from Ward (1977) ## ---------------------------- test1 <- t(matrix(c( 4, 2, 0, 1, 4, 1, 1, 1, 4), 3, 3)) expm.Higham08(test1) ## [,1] [,2] [,3] ## [1,] 147.86662244637000 183.76513864636857 71.79703239999643 ## [2,] 127.78108552318250 183.76513864636877 91.88256932318409 ## [3,] 127.78108552318204 163.67960172318047 111.96810624637124 ## -- these agree with ward (1977, p608) ## ---------------------------- ## Test case 2 from Ward (1977) ## ---------------------------- test2 <- t(matrix(c( 29.87942128909879, .7815750847907159, -2.289519314033932, .7815750847907159, 25.72656945571064, 8.680737820540137, -2.289519314033932, 8.680737820540137, 34.39400925519054), 3, 3)) expm.Higham08(test2) expm.Higham08(test2, balancing = FALSE) ## [,1] [,2] [,3] ##[1,] 5496313853692405 -18231880972009100 -30475770808580196 ##[2,] -18231880972009160 60605228702221760 101291842930249376 ##[3,] -30475770808580244 101291842930249200 169294411240850880 ## -- in this case a very similar degree of accuracy. ## ---------------------------- ## Test case 3 from Ward (1977) ## ---------------------------- test3 <- t(matrix(c( -131, 19, 18, -390, 56, 54, -387, 57, 52), 3, 3)) expm.Higham08(test3) expm.Higham08(test3, balancing = FALSE) ## [,1] [,2] [,3] ##[1,] -1.5096441587713636 0.36787943910439874 0.13533528117301735 ##[2,] -5.6325707997970271 1.47151775847745725 0.40600584351567010 ##[3,] -4.9349383260294299 1.10363831731417195 0.54134112675653534 ## -- agrees to 10dp with Ward (1977), p608. ??? (FIXME) ## ---------------------------- ## Test case 4 from Ward (1977) ## ---------------------------- test4 <- structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-10, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0), .Dim = c(10, 10)) E4 <- expm.Higham08(test4) Matrix(zapsmall(E4)) S4 <- as(test4, "sparseMatrix") # some R based expm() methods work for sparse: ES4 <- expm.Higham08(S4, bal=FALSE) stopifnot(all.equal(E4, unname(as.matrix(ES4)))) ## NOTE: Need much larger sparse matrices for sparse arith to be faster! ## ## example of computationally singular matrix ## m <- matrix(c(0,1,0,0), 2,2) eS <- expm.Higham08(m) # "works" (hmm ...) ```

expm documentation built on May 2, 2019, 5:25 p.m.