[R] matrix exponential: M^0

[Martyn Plummer]

On Tue, 2004-01-20 at 16:58, Peter Dalgaard wrote:
> Federico Calboli <f.calboli at ucl.ac.uk> writes:
> 
> > Dear All, 
> > 
> > I would like to ask why the zeroeth power of a matrix gives me a matrix
> > of ones rather than the identity matrix:
> 
> Because arithmetic on a matrix works element-wise. M^2 is not equal to
> M %*% M either (but is equal to  M*M).
> 
> (R doesn't have the matrix exponential function. Lifting expm() from
> Octave has been on my (virtual) TODO list for some time now...)

There is nothing in base R, but the Christopher Jackson's msm package
and Jim Lindsey's rmutils both have a matrix exponential function
(MatrixExp and mexp, respectively). Neither are particularly
sophisticated, using either the spectral decomposition or series
expansion.

Calculating the matrix exponential is harder than it looks (I'm sure
Peter knows this). In fact there is a classic paper by Moler and Van
Loan from the 1970s called "Nineteen dubious ways to calculate the
exponential of a matrix", which they updated last year in SIAM. 

There is an interesting looking library of Fortran routines for
calculating matrix exponentials called expokit:

http://www.maths.uq.edu.au/expokit

It has "no commercial use" license.  I'm going to use it in my own
software (JAGS) but it might also be useful in R if you can understand
what the license means.

Martyn