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

[Alberto Vieira Ferreira Monteiro]

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