36

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?

7 Answers 7

33

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

8

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
  • It doesn't have side effects in the global environment - C code is passed a reference to the R objects, so can modify an object in place. This is necessary for certain optimisations, but should never be exposed to the R user.
    – hadley
    Jul 19, 2010 at 1:58
  • @hadley I understand that. But if there is a single reference for two objects (as it seems to be the case above, probably for efficiency) and you let the C code modify the object in place, you are (I think) having side effects in the global environment, right? Jul 20, 2010 at 0:31
  • 2
    Your explanation is basically correct, but you are using suboptimal terminology. It doesn't make sense to talk about modifying the global environment here, because the object might not be in the global environment.
    – hadley
    Jul 20, 2010 at 15:30
7

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.

2

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.

2

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

2
  • 1
    this is not terribly useful, as the original bug was fixed more than two years ago ...
    – Ben Bolker
    Sep 25, 2013 at 19:01
  • plus this is a terrible way to perform matrix exponentiation for large exponents
    – m09
    Oct 14, 2013 at 13:10
2

You can simply use the Eigen values and Eigen vectors to compute the exponential of a matrix ;

# for a given matrix, A of power n

eig_vectors <- eigen(A)$vectors
eig_values <- eigen(A)$values

eig_vectors %*% diag(eig_values)^n %*% solve(eig_vectors)

Alternatively an improved answer from @MichaelChirico. The exponent 0 of a matrix will return its identity matrix instead of NULL.

pow = function(x, n) {
    if (n == 0) {
        I <- diag(length(diag(x)))
        return(I)        
    } 
    Reduce(`%*%`, replicate(n, x, simplify = FALSE))    
}
1
  • 3
    This answer may be correct, but it will be clearer if you add it to the explanation about it. Jan 24, 2021 at 11:03
0

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.

2
  • Calculating A%^%6 also leaves A as (initial A)%^%4. Assigning the result to a new variable, doesn't prevent my initial matrix from being changed.
    – gd047
    Jul 18, 2010 at 9:17
  • sounds like you just have to take the unusual step of assigning the matrix to a new variable first.
    – John
    Jul 18, 2010 at 13:01

Not the answer you're looking for? Browse other questions tagged or ask your own question.