Paul, Your solution based on SVD does not work 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: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] 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. ______________________________________________ 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.