[R] Exponentiate a matrix

[Duncan Murdoch]

On 9/21/2006 10:40 AM, Doran, Harold wrote:
> Suppose I have a square matrix P
> 
> P <- matrix(c(.3,.7, .7, .3), ncol=2)
> 
> I know that 
> 
>> P * P 
> 
> Returns the element by element product, whereas
> 
>> P%*%P
> 
> Returns the matrix product.
> 
> Now, P^2 also returns the element by element product. But, is there a
> slick way to write
> 
> P %*% P %*% P
> 
> Obviously, P^3 does not return the result I expect.


I don't think there's anything built in, but it's easy to write your own:

"%^%" <- function(mat, pow) {
   stopifnot(length(pow) == 1, all.equal(pow, round(pow)), nrow(mat) == 
ncol(mat))
   pow <- round(pow)
   if (pow < 0) {
     mat <- solve(mat)
     pow <- abs(pow)
   }
   result <- diag(nrow(mat))
   while (pow > 0) {
     result <- result %*% mat
     pow <- pow - 1
   }
   result
}

Now P %^% 3 will give you the matrix cube.

Duncan Murdoch