Matrix power in R

user_170792 0
Viewed 42235

Trying to compute the power of a matrix in R, I found that package expm implements the operator %^%.

So x %^% k computes the k-th power of a matrix. 

> A<-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3)

> A %^% 5
      [,1]  [,2] [,3]
[1,]  6469 18038 2929
[2,] 21837 60902 9889
[3,] 10440 29116 4729

but, to my surprise:

> A
     [,1] [,2] [,3]
[1,]  691 1926  312
[2,] 2331 6502 1056
[3,] 1116 3108  505

somehow the initial matrix A has changed to A %^% 4 !!!

How do you perform the matrix power operation?

linear-algebra r matrix
7 Answers

An inefficient version (since it's more efficient to first diagonalize your matrix) in base without much effort is:

pow = function(x, n) Reduce(`%*%`, replicate(n, x, simplify = FALSE))

I know this question is specifically about an old bug in expm, but it's one of the first results for "matrix power R" at the moment, so hopefully this little shorthand can be useful for someone else who ends up here just looking for a quick way to run matrix powers without installing any packages.

user_35769840

You can simply use the Eigen values and Eigen vectors to compute the exponential of a matrix ;

# for a given matrix, A of power n

eig_vectors <- eigen(A)$vectors
eig_values <- eigen(A)$values

eig_vectors %*% diag(eig_values)^n %*% solve(eig_vectors)

Alternatively an improved answer from @MichaelChirico. The exponent 0 of a matrix will return its identity matrix instead of NULL.

pow = function(x, n) {
    if (n == 0) {
        I <- diag(length(diag(x)))
        return(I)        
    } 
    Reduce(`%*%`, replicate(n, x, simplify = FALSE))    
}
user_122058330
Related