[R] matrix exponential: M^0
David Firth
d.firth at warwick.ac.uk
Thu Jan 22 17:12:24 CET 2004
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
More information about the R-help
mailing list