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

Martin Maechler maechler at stat.math.ethz.ch
Mon Aug 22 09:49:06 CEST 2011

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