A^k for matrix multiplication in R?

Question

Suppose A is some square matrix. How can I easily exponentiate this matrix in R?

I tried two ways already: Trial 1 with a for-loop hack and Trial 2 a bit more elegantly but it is still a far cry from A^k simplicity.

Trial 1

set.seed(10)
t(matrix(rnorm(16),ncol=4,nrow=4)) -> a 
for(i in 1:2){a <- a %*% a}

Trial 2

a <- t(matrix(c(0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0),nrow=4))
i <- diag(4) 
(function(n) {if (n<=1) a else (i+a) %*% Recall(n-1)})(10)

Answer 1

Indeed the expm's package does use exponentiation by squaring.

In pure r, this can be done rather efficiently like so,

"%^%" <- function(mat,power){
    base = mat
    out = diag(nrow(mat))
    while(power > 1){
        if(power %% 2 == 1){
            out = out %*% base
        }
        base = base %*% base
        power  = power %/% 2
    }
    out %*% base
}

Timing this,

m0 <- diag(1, nrow=3,ncol=3)
system.time(replicate(10000, m0%^%4000))#expm's %^% function
user  system elapsed 
0.31    0.00    0.31 
system.time(replicate(10000, m0%^%4000))# my %^% function
user  system elapsed 
0.28    0.00    0.28 

So, as expected, they are the same speed because they use the same algorithm. It looks like the overhead of the looping r code does not make a significant difference.

So, if you don't want to use expm, and need that performance, then you can just use this, if you don't mind looking at imperative code.