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.