Re: [R] A function for raising a matrix to a power? Atte, Your matrix A is not symmetric, that is why exponentiation using spectral decomposition, expM.sd, does not give you the correct answer.   Convert A to a symmetric matrix:  A <- (A + t(A))/2 then the results will all match. 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: Atte Tenkanen [mailto:attenka@utu.fi]  Sent: Wednesday, May 09, 2007 10:07 AM To: Paul Gilbert Cc: Ravi Varadhan; r-help@stat.math.ethz.ch Subject: Re: [R] A function for raising a matrix to a power? Hello, Thanks for many replys. I tested all the functions presented.  I'm a beginner in linear algebra, but now I have a serious question ;-) Here is a matrix A, which determinant is 3, so it is nonsingular.  Then there are similar computer runs done with each function proposed. I have calculated powers of A between A^274 - A^277. In the first version (see the outputs lower), when n=277,  there comes "zero in the upper right corner".  Can you say, if the first function is more reliable than others?  What can you say about the accuracy of calculations? -Atte #-------------------------------------------------------------# # Matrix A: A=rbind(c(-3,-4,-2),c(3,5,1),c(2,-1,4)) A det(A) # 1st version: "%^%"<-function(A,n){   if(n==1) A else {B<-A; for(i in (2:n)){A<-A%*%B}}; A   } for(i in 274:277){print(i);print(A%^%i)} # 2nd version: "%^%" <- function(A, n) if(n == 1) A else A %*% (A %^% (n-1)) for(i in 274:277){print(i);print(A%^%i)} # 3rd version: mp <- function(mat,pow){ ans <- mat for ( i in 1:(pow-1)){ ans <- mat%*%ans } return(ans) } for(i in 274:277){print(i);print(mp(A,i))} # 4th version: library(Malmig) mtx.exp for(i in 274:277){print(i);print(mtx.exp(A,i))} # 5th version: matrix.power <- function(mat, n) {   # test if mat is a square matrix   # treat n < 0 and n = 0 -- this is left as an exercise   # trap non-integer n and return an error   if (n == 1) return(mat)   result <- diag(1, ncol(mat))   while (n > 0) {     if (n %% 2 != 0) {       result <- result %*% mat       n <- n - 1     }     mat <- mat %*% mat     n <- n / 2   }   return(result) } for(i in 274:277){print(i);print(matrix.power(A,i))} # 6th version: expM.sd <- function(X,e){Xsd <- eigen(X); Xsd$vec %*% diag(Xsd$val^e) %*%t(Xsd$vec)} for(i in 274:277){print(i);print(expM.sd(A,i))} #-----------------OUTPUTS---------------------------# > A=rbind(c(-3,-4,-2),c(3,5,1),c(2,-1,4)) >  > A      [,1] [,2] [,3] [1,]   -3   -4   -2 [2,]    3    5    1 [3,]    2   -1    4 > det(A) [1] 3 >  > # 1st version: > "%^%"<-function(A,n){ +   if(n==1) A else {B<-A; for(i in (2:n)){A<-A%*%B}}; A +   } >  > for(i in 274:277){print(i);print(A%^%i)} [1] 274                [,1]           [,2]           [,3] [1,]  1.615642e+131  3.231283e+131 -3.940201e+115 [2,] -5.385472e+130 -1.077094e+131  4.925251e+114 [3,] -3.769831e+131 -7.539661e+131  7.880401e+115 [1] 275                [,1]           [,2]           [,3] [1,]  4.846925e+131  9.693850e+131 -1.182060e+116 [2,] -1.615642e+131 -3.231283e+131   0.000000e+00 [3,] -1.130949e+132 -2.261898e+132  4.728241e+116 [1] 276                [,1]           [,2]           [,3] [1,]  1.454078e+132  2.908155e+132 -1.576080e+116 [2,] -4.846925e+131 -9.693850e+131   0.000000e+00 [3,] -3.392848e+132 -6.785695e+132  1.576080e+117 [1] 277                [,1]           [,2]           [,3] [1,]  4.362233e+132  8.724465e+132   0.000000e+00 [2,] -1.454078e+132 -2.908155e+132 -1.576080e+116 [3,] -1.017854e+133 -2.035709e+133  3.782593e+117 >  >  > # 2nd version: > "%^%" <- function(A, n) if(n == 1) A else A %*% (A %^% (n-1)) >  > for(i in 274:277){print(i);print(A%^%i)} [1] 274                [,1]           [,2]           [,3] [1,]  1.615642e+131  3.231283e+131 -3.101125e+114 [2,] -5.385472e+130 -1.077094e+131  1.533306e+114 [3,] -3.769831e+131 -7.539661e+131  5.664274e+114 [1] 275                [,1]           [,2]           [,3] [1,]  4.846925e+131  9.693850e+131 -8.158398e+114 [2,] -1.615642e+131 -3.231283e+131  4.027430e+114 [3,] -1.130949e+132 -2.261898e+132  1.492154e+115 [1] 276                [,1]           [,2]           [,3] [1,]  1.454078e+132  2.908155e+132 -2.147761e+115 [2,] -4.846925e+131 -9.693850e+131  1.058350e+115 [3,] -3.392848e+132 -6.785695e+132  3.934193e+115 [1] 277                [,1]           [,2]           [,3] [1,]  4.362233e+132  8.724465e+132 -5.658504e+115 [2,] -1.454078e+132 -2.908155e+132  2.782660e+115 [3,] -1.017854e+133 -2.035709e+133  1.038290e+116 >  >  > # 3rd version: >  > mp <- function(mat,pow){ + ans <- mat + for ( i in 1:(pow-1)){ + ans <- mat%*%ans + } + return(ans) + } >  > for(i in 274:277){print(i);print(mp(A,i))} [1] 274                [,1]           [,2]           [,3] [1,]  1.615642e+131  3.231283e+131 -3.101125e+114 [2,] -5.385472e+130 -1.077094e+131  1.533306e+114 [3,] -3.769831e+131 -7.539661e+131  5.664274e+114 [1] 275                [,1]           [,2]           [,3] [1,]  4.846925e+131  9.693850e+131 -8.158398e+114 [2,] -1.615642e+131 -3.231283e+131  4.027430e+114 [3,] -1.130949e+132 -2.261898e+132  1.492154e+115 [1] 276                [,1]           [,2]           [,3] [1,]  1.454078e+132  2.908155e+132 -2.147761e+115 [2,] -4.846925e+131 -9.693850e+131  1.058350e+115 [3,] -3.392848e+132 -6.785695e+132  3.934193e+115 [1] 277                [,1]           [,2]           [,3] [1,]  4.362233e+132  8.724465e+132 -5.658504e+115 [2,] -1.454078e+132 -2.908155e+132  2.782660e+115 [3,] -1.017854e+133 -2.035709e+133  1.038290e+116 >  >  > # 4th version >  > library(Malmig) >  > mtx.exp function (X, n)  {     if (n != round(n)) {         n <- round(n)         warning("rounding exponent `n' to", n)     }     phi <- diag(nrow = nrow(X))     pot <- X     while (n > 0) {         if (n%%2)              phi <- phi %*% pot         n <- n%/%2         pot <- pot %*% pot     }     return(phi) } > for(i in 274:277){print(i);print(mtx.exp(A,i))} [1] 274                [,1]           [,2]           [,3] [1,]  1.615642e+131  3.231283e+131 -2.156709e+114 [2,] -5.385472e+130 -1.077094e+131  1.218501e+114 [3,] -3.769831e+131 -7.539661e+131  3.460636e+114 [1] 275                [,1]           [,2]           [,3] [1,]  4.846925e+131  9.693850e+131 -5.325150e+114 [2,] -1.615642e+131 -3.231283e+131  3.083014e+114 [3,] -1.130949e+132 -2.261898e+132  8.310626e+114 [1] 276                [,1]           [,2]           [,3] [1,]  1.454078e+132  2.908155e+132 -1.297786e+115 [2,] -4.846925e+131 -9.693850e+131  7.750249e+114 [3,] -3.392848e+132 -6.785695e+132  1.950919e+115 [1] 277                [,1]           [,2]           [,3] [1,]  4.362233e+132  8.724465e+132 -3.108580e+115 [2,] -1.454078e+132 -2.908155e+132  1.932685e+115 [3,] -1.017854e+133 -2.035709e+133  4.433079e+115 >  > # 5th version >  > matrix.power <- function(mat, n) + { +   # test if mat is a square matrix +   # treat n < 0 and n = 0 -- this is left as an exercise +   # trap non-integer n and return an error +   if (n == 1) return(mat) +   result <- diag(1, ncol(mat)) +   while (n > 0) { +     if (n %% 2 != 0) { +       result <- result %*% mat +       n <- n - 1 +     } +     mat <- mat %*% mat +     n <- n / 2 +   } +   return(result) + } >  > for(i in 274:277){print(i);print(matrix.power(A,i))} [1] 274                [,1]           [,2]           [,3] [1,]  1.615642e+131  3.231283e+131 -2.156709e+114 [2,] -5.385472e+130 -1.077094e+131  1.218501e+114 [3,] -3.769831e+131 -7.539661e+131  3.460636e+114 [1] 275                [,1]           [,2]           [,3] [1,]  4.846925e+131  9.693850e+131 -5.325150e+114 [2,] -1.615642e+131 -3.231283e+131  3.083014e+114 [3,] -1.130949e+132 -2.261898e+132  8.310626e+114 [1] 276                [,1]           [,2]           [,3] [1,]  1.454078e+132  2.908155e+132 -1.297786e+115 [2,] -4.846925e+131 -9.693850e+131  7.750249e+114 [3,] -3.392848e+132 -6.785695e+132  1.950919e+115 [1] 277                [,1]           [,2]           [,3] [1,]  4.362233e+132  8.724465e+132 -3.108580e+115 [2,] -1.454078e+132 -2.908155e+132  1.932685e+115 [3,] -1.017854e+133 -2.035709e+133  4.433079e+115 >  >  > # 6th versio >  > expM.sd <- function(X,e){Xsd <- eigen(X); Xsd$vec %*% diag(Xsd$val^e) %*%t(Xsd$vec)} >  > for(i in 274:277){print(i);print(expM.sd(A,i))} [1] 274                [,1]           [,2]           [,3] [1,]  8.215127e+129 -2.738376e+129 -1.916863e+130 [2,] -2.738376e+129  9.127919e+128  6.389543e+129 [3,] -1.916863e+130  6.389543e+129  4.472680e+130 [1] 275                [,1]           [,2]           [,3] [1,]  2.464538e+130 -8.215127e+129 -5.750589e+130 [2,] -8.215127e+129  2.738376e+129  1.916863e+130 [3,] -5.750589e+130  1.916863e+130  1.341804e+131 [1] 276                [,1]           [,2]           [,3] [1,]  7.393614e+130 -2.464538e+130 -1.725177e+131 [2,] -2.464538e+130  8.215127e+129  5.750589e+130 [3,] -1.725177e+131  5.750589e+130  4.025412e+131 [1] 277                [,1]           [,2]           [,3] [1,]  2.218084e+131 -7.393614e+130 -5.175530e+131 [2,] -7.393614e+130  2.464538e+130  1.725177e+131 [3,] -5.175530e+131  1.725177e+131  1.207624e+132 >  >  >  >  > 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 info...{{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.