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)[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.

Advertisements

6 thoughts on “Raising a matrix to a power in R

  1. 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!

  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.

  3. 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)

    }

Leave a Comment

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s