[Expm-developers] Euclidean (2-) norm and Matrix exponential

[Martin Maechler]

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

>
> Best,
> Ravi.
>
> ________________________________________
> 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]
>
>