# [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
>    >>
>    >>>>>>> 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:    >
>    >>
>    >>    > 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>
```