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

From: Alberto Monteiro <albmont_at_centroin.com.br>
Date: Wed, 07 May 2008 09:49:56 -0200

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



R-help_at_r-project.org 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. Received on Wed 07 May 2008 - 12:27:41 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 07 May 2008 - 12:30:35 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive