[Expm-developers] Matrix exponential exp(B * t) (for many t)
Christophe Dutang
dutangc at gmail.com
Mon Aug 22 13:00:15 CEST 2011
Hi all,
I just wonder if the method for computing exp(A) %*% v is the one based on
Krylov subspace, http://www.maths.uq.edu.au/expokit/guide.html , which is
implemented in matlab expokit http://www.maths.uq.edu.au/expokit/guide.html
?
A friend of mine who use matlab on a regular basis had told me the expokit
implementation is particularly efficient.
Regards
Christophe
2011/8/22 Martin Maechler <maechler at stat.math.ethz.ch>
> >>>>> Martin Maechler <maechler at stat.math.ethz.ch>
> >>>>> on Sun, 21 Aug 2011 21:33:23 +0200 writes:
>
> > Dear Ravi, On Sat, Aug 20, 2011 at 16:56, Ravi Varadhan
> > <rvaradhan at jhmi.edu> wrote:
> >> Dear Martin,
> >>
> >> Thanks for pointing me to "expm" package. I noted that
> >> it has a rich array of algorithms for exponentiating a
> >> matrix. However, it does not seem to have a way to
> >> efficiently compute exp(A) %*% x0, where x0 is a vector.
>
> > Yes, you are perfectly right and we (the expm-package
> > authors) have been aware that "we need this" but none of
> > us has come around to add it to expm yet.
>
> >> This is the solution of a linear system of ODE: dx/dt =
> >> B, x(0) = x0 (note: A = B * t). Computing exp(B*t) and
> >> then doing the matrix-vector product is inefficient
> >> especially for large matrices.
>
> > yes, and in particular if it should happen --- as
> > typically --- for a whole set of t values.
>
> >> I have gotten interested in this in response to a query
> >> that John Nash had recently posted in r-devel. I have
> >> translated the matlab function by Sidjek to do this. Let
> >> me know if you want this feature in expm().
>
> > Yes, that's really something we should add to the expm
> > package.. but may be not to the expm() function per se.
> > I think a main use case should be to pass a whole vector t
> > values (and a fixed x0 vector), right? I'd definitely
> > look at what you sent us - at the beginning of this
> > thread, to see how it fits into this framework. Thank you
> > very much, Ravi, for "insisting" ;-b Martin
>
> Hmm. Today I had a careful look at the Ravi-pexp.R that you
> had sent me, and I really don't see how your pexp() would solve the
> above problem....
>
> but maybe I misunderstood and you meant something else?
>
> Martin
>
>
> >> ________________________________________ From: Ravi
> >> Varadhan Sent: Friday, August 19, 2011 5:31 PM To:
> >> 'Martin Maechler' Cc: 'Doug and Martin Subject: RE:
> >> Euclidean (2-) norm and Matrix exponential
> >>
> >> Dear Martin,
> >>
> >> Thanks for the clarifications and the new pointers. I
> >> concur on all the points that you have made. I will take
> >> a look at the documentation of these functions that you
> >> mentioned.
> >>
> >> It is interesting, the point that that you make about the
> >> lack of usefulness of 2-norm. Is it really the case this
> >> only serves theoretical purposes and that there is no
> >> *practical* need for this norm?
> >>
> >> In any case, I have a neat little function that uses the
> >> power method + sequence acceleration to estimate 2-norm
> >> in a somewhat efficient manner.
> >>
> >> Best regards, Ravi.
> >>
> >> -------------------------------------------------------
> >> Ravi Varadhan, Ph.D. Assistant Professor, Division of
> >> Geriatric Medicine and Gerontology School of Medicine
> >> Johns Hopkins University
> >>
> >> Ph. (410) 502-2619 email: rvaradhan at jhmi.edu
> >>
> >> -----Original Message----- From: Martin Maechler
> >> [mailto:maechler at stat.math.ethz.ch] Sent: Friday, August
> >> 19, 2011 5:23 PM To: Ravi Varadhan Cc: 'Doug and Martin
> >> Subject: Re: Euclidean (2-) norm and Matrix exponential
> >>
> >>>>>>> Ravi Varadhan <rvaradhan at jhmi.edu> on Fri, 19
> >>>>>>> Aug 2011 17:52:24 +0000 writes:
> >>
> >> > Dear Martin & Doug, I have a simple R function for
> >> > computing matrix exponential that is based on a >
> >> Pade-Polynomial approximation of Sidjek. I translated it
> >> > to R from his matlab code. It appears to be simpler
> >> than > your code in "Matrix", but my limited testing
> >> shows that > it produces identical results, and is a
> >> bit faster.
> >>
> >> Dear Ravi, you must have overlooked the fact that Doug
> >> and I (and two others) have been authoring a package
> >> 'expm' (on R-forge and CRAN) which contains newer and
> >> partly considerably better algorithms for expm() {and
> >> also some for logm() and sqrtm()}.
> >>
> >> and the topic appeared on R-devel / R-help several times
> >> in the last couple of years... but you are "well
> >> excused": You don't see anything at all about that when
> >> you read help(expm) from Matrix ==> I have added a note
> >> (and \seealso{..} \link's) to the Matrix help page. So
> >> the next Matrix version will point to these .... thanks
> >> to you!
> >>
> >> > I also noted that there is no function in base R and
> >> in > Matrix for computing the Euclidean norm or 2-norm
> >> of a > matrix. This is trivially easy, but
> >> surprisingly there is > nothing available in R (at
> >> least as far as I know). I > have also written a
> >> simple function for this as well as > computing the
> >> condition number of a matrix that is simpler > than
> >> condest() in Matrix.
> >>
> >> Well, but did you see kappa() and rcond() in addition to
> >> condest() ?
> >>
> >> So kappa() --- the oldest of *all* the functions
> >> mentioned --- has been in S I think even before S-plus!
> >> --- already does the 2-norm condition number. You are
> >> right that the 2-norm itself is not available,
> >> explicitly, but then it is AFAIR considerably slower than
> >> the other norms, and "everybody should now" how to
> >> trivially compute it, if they really want. But as said,
> >> I think it is +- clear to numerical analysis experts that
> >> you work with the 1- and possibly Inf-norm instead, as
> >> both are considerably simpler.
> >>
> >> OTOH, you make a good point that this should be better
> >> documented.. ... though I wonder why you didn't follow
> >> the "seealso" links on the condest or (Matrix) rcond help
> >> page, pointing to the base rcond + kappa help page...
> >>
> >> We'd be grateful if you can propose documentation
> >> improvements here.
> >>
> >> Best regards, Martin
> >>
> >> > Perhaps, you might want to consider these for "base"
> >> or > even "Matrix". In any case, I would be
> >> interested in your > thoughts.
> >>
> >> > Best regards, Ravi.
> >>
> >> >
> >> -------------------------------------------------------
> >> > Ravi Varadhan, Ph.D. Assistant Professor, Division of
> >> > Geriatric Medicine and Gerontology School of
> >> Medicine > Johns Hopkins University
> >>
> >> > Ph. (410) 502-2619 email: >
> >> rvaradhan at jhmi.edu<mailto:rvaradhan at jhmi.edu>
> >>
> >> > xpexpm.R, pexpm.R [Click mouse-2 to save to a file]
> >>
> >>
> _______________________________________________
> Expm-developers mailing list
> Expm-developers at lists.r-forge.r-project.org
>
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/expm-developers
>
--
Christophe DUTANG
Ph. D. student at ISFA, Lyon, France
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/expm-developers/attachments/20110822/69bca338/attachment.html>
More information about the Expm-developers
mailing list