[R] matrix exponential: M^0

[David Firth]

Prompted by this thread, I have tidied up a Fortran program I wrote 
with Marina Shapira.  We would be happy for this ("mexp") to become 
part of R, either as a contributed package or as part of the base 
distribution if it's good enough.  I have packaged it and put it at
   http://www.warwick.ac.uk/go/dfirth/software/mexp

The examples in the help file show results on some "difficult" (but 
small) test matrices from the literature.  I would welcome any feedback 
on accuracy in other test cases.

This mexp function provides a choice of methods, Taylor series and Padé 
approximation, both with the usual "squaring and scaling" for increased 
accuracy.

I have not tested it on any platform other than my own (Mac OS X).  
Until I know that it compiles and works on other platforms there seems 
little point in submitting it as a CRAN package.

David


On Thursday, Jan 22, 2004, at 09:33 Europe/London, Martin Maechler 
wrote:

>>>>>> "PD" == Peter Dalgaard <p.dalgaard at biostat.ku.dk>
>>>>>>     on 21 Jan 2004 19:08:38 +0100 writes:
>
>     PD> Martyn Plummer <plummer at iarc.fr> writes:
>>> 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.
>
>     PD> Right (magnificent paper by the way), although I
>     PD> actually hadn't heard about the update. As I remember
>     PD> it, Octave implements what Moler+v.Loan ends up
>     PD> suggesting in the 1978 paper.
>
> The update is actually available online
> from http://epubs.siam.org/sam-bin/dbq/article/41801
> with the extended title "...., 25 Years Later" .
>
> The extension is 8 pages of text + 1.2 pages of references,
> in which (p.42) they say
> ``The matrix exponential is an important computational tool in
>   control theory, so availability of  expm(A)  in early versions
>   of Matlab quite possibily contributed to the system's technical
>   and commercial success.''
>
> and they also tell how expm() is implemented in Matlab
> (scaling, squaring, Padé approximation) which is presumably what
> octave does too?
>
> -----
>
> But -- going back to our original subject --
> Isn't
>
>   expm(A) =  ``e ^ A'' := sum_{n=0}^\Inf  A^n / n!
>
> definitely a bit harder than  A^n for (non-negative) integer n ?
>
> Martin
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html