[Expm-developers] Matrix exponential exp(B * t) (for many t)

[Martin Maechler]

>>>>> 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]
    >> 
    >>