Re: [R] matrix exponential: M^0

David Firth Thu, 22 Jan 2004 09:31:10 -0800

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 <[EMAIL PROTECTED]>
    on 21 Jan 2004 19:08:38 +0100 writes:

PD> Martyn Plummer <[EMAIL PROTECTED]> 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

______________________________________________ [EMAIL PROTECTED] 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


______________________________________________
[EMAIL PROTECTED] 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

Re: [R] matrix exponential: M^0

Reply via email to