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)[1] != dim(A)[2]) 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)[1],dim(A)[2]))
# 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.
Victor A. Barger 
Very nice!. just what I was looking for.
Keep up the good work :-).
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!
Great, thanks for this. What about defining mpower for negative exponents, e.g.
mpower(M,-p^2)= solve(mpower(M,p^2))
Very nice, thanks for sharing this very effective bit of code. It is strange that this was never added to the base of R.
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)[1]
n <- dim(A)[2]
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)
}
Thanks for sharing! I meet the same problem with R today.
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.