Matrix power in R
Solution 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
Solution 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.
Solution 3
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.
Solution 4
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 :
-
result
, in order to store the output, -
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
.
Solution 5
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
George Dontas
Updated on January 26, 2021Comments
-
George Dontas over 3 years
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?