# 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.