[Expm-developers] Matrix exponential exp(B * t) (for many t)
Martin Maechler
maechler at stat.math.ethz.ch
Mon Aug 22 17:47:01 CEST 2011
>>>>> Ravi Varadhan <rvaradhan at jhmi.edu>
>>>>> on Mon, 22 Aug 2011 14:05:58 +0000 writes:
> Yes, Christophe. It is exactly that implementation by
> Sidje that I have translated into R. Ravi.
Excellent. Thank you, Ravi, and excuse my misunderstanding!
Martin
> -------------------------------------------------------
> 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>
> From: Christophe Dutang [mailto:dutangc at gmail.com] Sent:
> Monday, August 22, 2011 7:00 AM To: Martin Maechler; Ravi
> Varadhan Cc: expm-developers Subject: Re:
> [Expm-developers] Matrix exponential exp(B * t) (for many
> t)
> 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<mailto:maechler at stat.math.ethz.ch>>
>>>>> Martin Maechler
> <maechler at stat.math.ethz.ch<mailto: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<mailto: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<tel:%28410%29%20502-2619> email:
>>> rvaradhan at jhmi.edu<mailto:rvaradhan at jhmi.edu>
>>>
>>> -----Original Message----- From: Martin Maechler
>>> [mailto:maechler at stat.math.ethz.ch<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<mailto: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<tel:%28410%29%20502-2619> email: >
>>> rvaradhan at jhmi.edu<mailto:rvaradhan at jhmi.edu><mailto: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<mailto: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
More information about the Expm-developers
mailing list