# Raising a matrix to a power in R

A surprising deficiency in R is the absence of an operator or built-in function for matrix exponentiation. I had to deal with this today when porting a function from MATLAB to R. Here is my quick-and-dirty solution:

```mpower = function(M,p) {
A = as.matrix(M)

if (dim(A) != dim(A)) stop("not a square matrix")

# M^{-1} = the matrix inverse of M
if (p==-1) return(solve(A))

# M^0 = I
if (p==0) return(diag(1,dim(A),dim(A)))

# M^1 = M
if (p==1) return(A)

if (p < -1) stop("only powers >= -1 allowed")
if (p != as.integer(p)) stop("only integer powers allowed")

R = A
for (i in 2:p) {
R = R %*% A
}
return(R)
}
```

Not the most efficient method, perhaps, but it works. If you know of a better alternative, please add a comment to this blog entry. One possibility might be to port MATLAB’s `expm` function, as defined in `expm1.m` in the MATLAB demos directory, to R.

For a detailed review of methods for computing the exponential of a matrix, see:

Moler, Cleve and Charles Van Loan (2003), “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,” Society for Industrial and Applied Mathematics, 45 (1), 3-49.

## 7 thoughts on “Raising a matrix to a power in R”

1. Ami B. on said:

Very nice!. just what I was looking for.
Keep up the good work :-).

2. Andrew on said:

I have been lost on this one for a few months and hadn’t come up with a simple solution. This isn’t complete but it’s excellent for my purposes. Thank you!

3. Tudor on said:

Great, thanks for this. What about defining mpower for negative exponents, e.g.

mpower(M,-p^2)= solve(mpower(M,p^2))

4. Jasja Dekker on said:

Very nice, thanks for sharing this very effective bit of code. It is strange that this was never added to the base of R.

5. Daniel on said:

Thanks, Victor. That’s very helpful, especially for Markov chains. However, if you had to compute A^5 by hand, first you’d compute A^2, then square the result to get A^4 in two steps, and finally multiply by A once more. I’ve generalized this approach by adding an additional function that, for a nonnegative integer, returns an array of binary digits. This should involve about one half log base 2 of p matrix multiplications rather than p-1. Sorry by the way for the lack of comments — documentation is not my strong suit.

dec2bin = function(n){

p <- n

if(p==0) return(0)

bin.rep <- NULL

k 0){

q <- p%%2^k
bin.rep <- c(sign(q),bin.rep)
p <- p – q
k <- k+1

}

return(bin.rep)

}

eye = function(n){

return(diag(1,n,n))

}

mpower = function(M,p){

A <- as.matrix(M)

m <- dim(A)
n <- dim(A)

if(m != n) stop("Matrix must be square.")

if(p == -1) return(solve(A))

p <- dec2bin(p)

l <- length(p)

B <- eye(n)

for(k in l:1){

if(p[k]) B <- B%*%A
A <- A%*%A

}

return(B)

}

6. Shihai on said:

Thanks for sharing! I meet the same problem with R today.

7. Martin Maechler Sr. on said:

The only comprehensive solution is function Matpow from the R package powerplus. All other functions (expm, etc) have huge limitations, like integer exponents. With Matpow you have no limitations.