Re: R-alpha: matrix exponentiation

Bill Venables (bill@stats.ox.ac.uk)
Wed, 20 Nov 1996 18:00:24 GMT


Date: Wed, 20 Nov 1996 18:00:24 GMT
Message-Id: <199611201800.SAA19316@brunhilde.stats>
From: Bill Venables <bill@stats.ox.ac.uk>
To: thomas@biostat.washington.edu
Subject: Re: R-alpha: matrix exponentiation


 >   On Wed, 20 Nov 1996, Jim Lindsey wrote:
 >
 >   > 3.Here is a matrix exponentiation function for anyone working on
 >   >   stochastic processes.  If anyone knows a better algorithm than brute
 >   >   force, please let me know. Jim
 >
 >   There is a *much* better algorithm, at least for symmetric
 >   matrices.  Get the eigendecomposition and work with that
 >   (this function needs argument checking, of course)

 >   mexp2<-function (mat, n) {
 >	 ei <- eigen(mat)
 >	 t(ei$vectors) %*% diag(exp(ei$values)) %*% ei$vectors
 >   }

You have the wrong term transposed, the diagonal matrix product
can be done more efficiently and `n' is not used.  How about:

     mexp2 <- function(mat) {
	ei <- eigen(mat)
	ei$vectors %*% (exp(ei$values) * t(ei$vectors))
     }

 >   I think that you can do something similar with the singular
 >   value decomposition for general square matrices but I would
 >   have to look it up.

I don't think so; you do need the eigenvalue decomposition.  The
SVD could be used to calculate exp(t(X) %*% X) directly from X.

If you want an exhaustive put-down of nearly all conceivable
methods of calculating matrix exponentials, you might like to
look at:

C. B. Moler and C. F. van Loan. (1978)  Nineteen dubious ways to
compute the matrix exponential.  SIAM Review, 20, 801--836.

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-testers mailing list -- For info or help, send "info" or "help",
To [un]subscribe, send "[un]subscribe"
(in the "body", not the subject !)  To: r-testers-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-