View source: R/exponentiate_matrix.R

exponentiate_matrix | R Documentation |

Calculate the exponential `\exp(T\cdot A)`

of some quadratic real-valued matrix A for one or more scalar scaling factors T.

```
exponentiate_matrix(A, scalings=1, max_absolute_error=1e-3,
min_polynomials=1, max_polynomials=1000)
```

`A` |
A real-valued quadratic matrix of size N x N. |

`scalings` |
Vector of real-valued scalar scaling factors T, for each of which the exponential |

`max_absolute_error` |
Maximum allowed absolute error for the returned approximations. A smaller allowed error implies a greater computational cost as more matrix polynomials need to be included (see below). The returned approximations may have a greater error if the parameter |

`min_polynomials` |
Minimum number of polynomials to include in the approximations (see equation below), even if |

`max_polynomials` |
Maximum allowed number of polynomials to include in the approximations (see equation below). Meant to provide a safety limit for the amount of memory and the computation time required. For example, a value of 1000 means that up to 1000 matrices (powers of A) of size N x N may be computed and stored temporarily in memory. Note that if |

Discrete character evolution Markov models often involve repeated exponentiations of the same transition matrix along each edge of the tree (i.e. with the scaling T being the edge length). Matrix exponentiation can become a serious computational bottleneck for larger trees or large matrices (i.e. spanning multiple discrete states). This function pre-calculates polynomials `A^p/p!`

of the matrix, and then uses linear combinations of the same polynomials for each requested T:

```
\exp(T\cdot A) = \sum_{p=0}^P T^p\frac{A^p}{p!} + ...
```

This function thus becomes very efficient when the number of scaling factors is large (e.g. >10,000).
The number of polynomials included is determined based on the specified `max_absolute_error`

, and based on the largest (by magnitude) scaling factor requested. The function utilizes the balancing algorithm proposed by James et al (2014, Algorithm 3) and the scaling & squaring method (Moler and Van Loan, 2003) to improve the conditioning of the matrix prior to exponentiation.

A 3D numeric matrix of size N x N x S, where N is the number of rows & column of the input matrix A and S is the length of `scalings`

. The [r,c,s]-th element of this matrix is the entry in the r-th row and c-th column of `\exp(scalings[s]\cdot A)`

.

Stilianos Louca

R. James, J. Langou and B. R. Lowery (2014). On matrix balancing and eigenvector computation. arXiv:1401.5766

C. Moler and C. Van Loan (2003). Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review. 45:3-49.

```
# create a random 5 x 5 matrix
A = get_random_mk_transition_matrix(Nstates=5, rate_model="ER")
# calculate exponentials exp(0.1*A) and exp(10*A)
exponentials = exponentiate_matrix(A, scalings=c(0.1,10))
# print 1st exponential: exp(0.1*A)
print(exponentials[,,1])
# print 2nd exponential: exp(10*A)
print(exponentials[,,2])
```

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.