Matrix power in R

This is not a proper answer, but may be a good place to have this discussion and understand the inner workings of R. This sort of bug has crept up before in another package I was using.

First, note that simply assigning the matrix to a new variable first does not help:

> A <- B <-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3)
> r1 <- A %^% 5
> A
     [,1] [,2] [,3]
[1,]  691 1926  312
[2,] 2331 6502 1056
[3,] 1116 3108  505
> B
     [,1] [,2] [,3]
[1,]  691 1926  312
[2,] 2331 6502 1056
[3,] 1116 3108  505

My guess is that R is trying to be smart passing by reference instead of values. To actually get this to work you need to do something to differentiate A from B:

`%m%` <- function(x, k) {
    tmp <- x*1
    res <- tmp%^%k
    res
}
> B <-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3)
> r2 <- B %m% 5
> B
     [,1] [,2] [,3]
[1,]    1    2    1
[2,]    3    8    1
[3,]    0    4    1

What is the explicit way of doing this?

Finally, in the C code for the package, there is this comment:

• NB: x will be altered! The caller must make a copy if needed

But I don't understand why R lets C/Fortran code have side effects in the global environment.

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.

Although the source-code is not visible in the package since it is packed in a .dll file, I believe the algorithm used by the package is the fast exponentiation algorithm, which you can study by looking at the function called matpowfast instead.

You need two variables :

1. result, in order to store the output,
2. mat, as an intermediate variable.

To compute A^6, since 6 = 110 (binary writing), in the end, result = A^6 and mat = A^4. This is the same for A^5.

You could easily check if mat = A^8 when you try to compute A^n for any 8<n<16. If so, you have your explanation.

The package function uses the initial variable A as the intermediate variable mat.

I have fixed that bug in the R-forge sources (of "expm" package), svn rev. 53. --> expm R-forge page For some reason the web page still shows rev.52, so the following may not yet solve your problem (but should within 24 hours):

 install.packages("expm", repos="http://R-Forge.R-project.org")

Otherwise, get the svn version directly, and install yourself:

 svn checkout svn://svn.r-forge.r-project.org/svnroot/expm

Thanks to "gd047" who alerted me to the problem by e-mail. Note that R-forge also has its own bug tracking facilities.
Martint