Re: [R] A function for raising a matrix to a power? Ron E. VanNimwegen wrote: > > I was looking for a similar operator, because R uses scalar products > when raising a matrix to a power with "^".  There might be something > more elegant, but this little loop function will do what you need for a > matrix "mat" raised to a power "pow": > > mp <- function(mat,pow){ > ans <- mat > for ( i in 1:(pow-1)){ > ans <- mat%*%ans > } > return(ans) > } > This function is extremely inefficient for high powers of pow [an unhappy choice of variables, because pow is the power function in many languages] A better method would keep track of the powers of two, and would optimize according to it. For example, in order to compute mat^42, we just have to compute mat^2, mat^4, mat^8, mat^16, mat^32, and then mat^40 and mat^42, with "only" 7 multiplications, instead of 41. See the technique in... http://en.wikipedia.org/wiki/Exponentiation_by_squaring 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) } # Test case: mat <- matrix(c(1,1,1,0), 2, 2) mat42 <- matrix.power(mat, 42) det(mat) # -1 det(mat42) # not 1, because of truncation errors, but close enough :-) # I hate floating point arithmetics :-/ # Another test case: mat <- matrix(c(1,1,0,1), 2, 2) mat314159 <- matrix.power(mat, 314159) # mat314159 is matrix(c(1,314159,0,1), 2, 2) Alberto Monteiro ______________________________________________ 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.