Re: [R] A function for raising a matrix to a power?

From: Alberto Vieira Ferreira Monteiro <albmont_at_centroin.com.br>
Date: Sun, 06 May 2007 16:19:49 +0000

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_at_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. Received on Sun 06 May 2007 - 16:29:21 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 Sun 06 May 2007 - 20:33:20 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.