[R] How do I write a power of matrices?? [was: How do I write a sum of matrixes??]

[Alberto Monteiro]

Jorge Ivan Velez wrote:
> 
> I think the function could be better but try this:
> 
> # Function: M is your matrix and n MUST be an integer>0
> mat.pow<-function(M,n) {
>   result<-M
>     if(n>1){
>        for ( iter in 2:n) result<-M%*%result
>        result
>            }
>     else {result}
>   result
>  }
> 
There are much more efficient ways to compute a power,
just check:
http://en.wikipedia.org/wiki/Exponentiation_by_squaring

Or, translating the pseudo-code to R:

mat.pow <- function(M, n) {
  result <- diag(1, ncol(M))
  while (n > 0) {
    if (n %% 2 == 1) {  
      result <- result %*% M
      n <- n - 1
    }
    M <- M %*% M
    n <- n / 2
  }
  return(result)
}

# The matrix
m <-  rbind(c(1,1,0), c(0,1,1), c(0,0,1))
 
# Goal m^6 = m x m x m x m x m x m
goal=  m %*% m %*% m %*% m %*% m %*% m
 
# matpow
res=mat.pow(m,6)
 
# Check point
all.equal(goal,res)

This algorithm would be fast, unless n is a _very_ big number.

Alberto Monteiro