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