R中的矩阵功率

Question

试图计算R中矩阵的幂,我发现包expm实现了运算符%^%.

所以x%^%k计算矩阵的k次方.

> 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

但令我惊讶的是:

> A
     [,1] [,2] [,3]
[1,]  691 1926  312
[2,] 2331 6502 1056
[3,] 1116 3108  505

不知何故,初始矩阵A已变为A%^%4 !!!

你如何执行矩阵功率操作？

Answer 1

我修复了R-forge源("expm"包),svn rev中的那个bug.53. - > expm R-forge页面 由于某种原因,网页仍显示rev.52,因此以下内容可能尚未解决您的问题(但应在24小时内):

 install.packages("expm", repos="http://R-Forge.R-project.org")

否则,直接获取svn版本并自行安装:

 svn checkout svn://svn.r-forge.r-project.org/svnroot/expm

感谢"gd047",他通过电子邮件提醒我这个问题.请注意,R-forge也有自己的错误跟踪工具.
Martint

Answer 2

这不是一个正确的答案,但可能是一个讨论和理解R的内部工作的好地方.这种错误在我正在使用的另一个包中之前已经悄悄上升.

首先,请注意,首先简单地将矩阵分配给新变量没有帮助:

> 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

我的猜测是R试图通过引用而不是值来进行智能传递.要真正实现这一点,您需要做一些事情来区分A和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

这样做的明确方法是什么？

最后,在包的C代码中,有这样的评论:

• 注意:x会被改变!如果需要,呼叫者必须复制

但我不明白为什么R让C/Fortran代码在全局环境中有副作用.

Answer 3

一个低效的版本（因为首先对矩阵进行对角化更有效）base是：

pow = function(x, n) Reduce(`%*%`, replicate(n, x, simplify = FALSE))

我知道这个问题具体是关于 中的一个老错误expm，但它是目前“矩阵幂 R”的第一个结果之一，所以希望这个小速记对于其他最终在这里只是寻找快速方法的人有用无需安装任何软件包即可运行矩阵幂。