Matrix power in R

问题:

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?

回答1:

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

回答2:

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.

回答3:

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.

回答4:

Very quick solution without using any package is using recursivity: if your matrix is a

 powA = function(n)
 {
    if (n==1)  return (a)
    if (n==2)  return (a%*%a)
    if (n>2) return ( a%*%powA(n-1))
 }

HTH

回答5:

A^5 = (A^4)*A

I suppose the library mutates the original variable, A, so that the each step involves multiplying the result-up-till-then with the original matrix, A. The result you get back seem fine, just assign them to a new variable.