Re: [R] A function for raising a matrix to a power? Ravi Varadhan wrote: > Paul, >  > Your solution based on SVD does not work  Ooops. I really am getting rusty. The idea is based on eigenvalues  which, of course, are not always the same as singular values. Paul even for the matrix in your example > (the reason it worked for e=3, was that it is an odd power and since P is a > permutation matrix.  It will also work for all odd powers, but not for even > powers). >    > However, a spectral decomposition will work for all symmetric matrices (SVD > based exponentiation doesn't work for any matrix).   >  > Here is the function to do exponentiation based on spectral decomposition: >  > expM.sd <- function(X,e){Xsd <- eigen(X); Xsd$vec %*% diag(Xsd$val^e) %*% > t(Xsd$vec)} >  >> P <- matrix(c(.3,.7, .7, .3), ncol=2) >> P >      [,1] [,2] > [1,]  0.3  0.7 > [2,]  0.7  0.3 >> P %*% P >      [,1] [,2] > [1,] 0.58 0.42 > [2,] 0.42 0.58 >> expM(P,2)   >      [,1] [,2] > [1,] 0.42 0.58 > [2,] 0.58 0.42 >> expM.sd(P,2) >      [,1] [,2] > [1,] 0.58 0.42 > [2,] 0.42 0.58 >  > Exponentiation based on spectral decomposition should be generally more > accurate (not sure about the speed).  A caveat is that it will only work for > real, symmetric or complex, hermitian matrices.  It will not work for > asymmetric matrices.  It would work for asymmetric matrix if the > eigenvectors are made to be orthonormal. >  > Ravi. >  > ---------------------------------------------------------------------------- > ------- >  > Ravi Varadhan, Ph.D. >  > Assistant Professor, The Center on Aging and Health >  > Division of Geriatric Medicine and Gerontology  >  > Johns Hopkins University >  > Ph: (410) 502-2619 >  > Fax: (410) 614-9625 >  > Email: rvaradhan@jhmi.edu >  > Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html >  >   >  > ---------------------------------------------------------------------------- > -------- >  >  > -----Original Message----- > From: r-help-bounces@stat.math.ethz.ch > [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Paul Gilbert > Sent: Monday, May 07, 2007 5:16 PM > To: Atte Tenkanen > Cc: r-help@stat.math.ethz.ch > Subject: Re: [R] A function for raising a matrix to a power? >  > You might try this, from 9/22/2006 with subject line Exponentiate a matrix: >  > I am getting a bit rusty on some of these things, but I seem to recall  > that there is a numerical advantage (speed and/or accuracy?) to  > diagonalizing: >  >  > expM <- function(X,e) { v <- La.svd(X); v$u %*% diag(v$d^e) %*% v$vt } >  >  > P <- matrix(c(.3,.7, .7, .3), ncol=2) >  > P %*% P %*% P >       [,1]  [,2] > [1,] 0.468 0.532 > [2,] 0.532 0.468 >  > expM(P,3) >       [,1]  [,2] > [1,] 0.468 0.532 > [2,] 0.532 0.468 >  > I think this also works for non-integer, negative, large, and complex  > exponents: >  >  > expM(P, 1.5) >           [,1]      [,2] > [1,] 0.3735089 0.6264911 > [2,] 0.6264911 0.3735089 >  > expM(P, 1000) >      [,1] [,2] > [1,]  0.5  0.5 > [2,]  0.5  0.5 >  > expM(P, -3) >         [,1]    [,2] > [1,] -7.3125  8.3125 > [2,]  8.3125 -7.3125 >  > expM(P, 3+.5i) >                   [,1]              [,2] > [1,] 0.4713+0.0141531i 0.5287-0.0141531i > [2,] 0.5287-0.0141531i 0.4713+0.0141531i >  > >  > Paul Gilbert >  > Doran, Harold wrote: >  >  > Suppose I have a square matrix P >  > >  > P <- matrix(c(.3,.7, .7, .3), ncol=2) >  > >  > I know that >  > >  > >  >> P * P >  > >  > Returns the element by element product, whereas >  > >  > >  > >  >> P%*%P >  >> >  > >  > Returns the matrix product. >  > >  > Now, P2 also returns the element by element product. But, is there a >  > slick way to write >  > >  > P %*% P %*% P >  > >  > Obviously, P3 does not return the result I expect. >  > >  > Thanks, >  > Harold >  > >  >  > Atte Tenkanen wrote: >> Hi, >> >> Is there a function for raising a matrix to a power? For example if you > like to compute A%*%A%*%A, is there an abbreviation similar to A^3? >> Atte Tenkanen >> >>> A=rbind(c(1,1),c(-1,-2)) >>> A >>      [,1] [,2] >> [1,]    1    1 >> [2,]   -1   -2 >>> A^3 >>      [,1] [,2] >> [1,]    1    1 >> [2,]   -1   -8 >> >> But: >> >>> A%*%A%*%A >>      [,1] [,2] >> [1,]    1    2 >> [2,]   -2   -5 >> >> ______________________________________________ >> R-help@stat.math.ethz.ch mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > ============================================================================ > ======== >  > La version franÃ§aise suit le texte anglais. >  > ---------------------------------------------------------------------------- > -------- >  > This email may contain privileged and/or confidential inform...{{dropped}} >  > ______________________________________________ > R-help@stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. ==================================================================================== La version franÃ§aise suit le texte anglais. ------------------------------------------------------------------------------------ This email may contain privileged and/or confidential inform...{{dropped}} ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.