Matrix power in R

Matrix power in R

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

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 dont understand why R lets C/Fortran code have side effects in the global environment.

Matrix power in R

An inefficient version (since its 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 its 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.

Leave a Reply

Your email address will not be published.