logo标签列表

在R中进行矩阵幂运算

rmatrixlinear-algebra
36

在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 !!!

如何执行矩阵的幂运算？

7个回答
33

我已经在"expm"软件包的R-forge源代码中修复了那个错误，svn版本号为53。--> expm R-forge页面 由于某种原因，网页仍然显示版本号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

8
这不是一个正式的答案，但可能是一个好的讨论场所，以了解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代码中，有这样的注释：

  • NB：x将被更改！如果需要，调用方必须进行复制

但我不明白为什么R允许C / Fortran代码在全局环境中产生副作用。

它不会在全局环境中产生副作用 - C代码被传递一个对R对象的引用，因此可以直接修改对象。这对于某些优化是必要的，但永远不应该暴露给R用户。 - hadley
@hadley 我理解。但是，如果有两个对象的单个引用（正如上面似乎是的，可能是为了提高效率），并且您让C代码直接修改对象，那么您会在全局环境中产生副作用，对吗？ - Eduardo Leoni
2
你的解释基本上是正确的，但你使用了次优术语。在这里谈论修改全局环境并没有意义，因为对象可能不在全局环境中。 - hadley
7
在不费太多力气的情况下，base中一个效率不高的版本（因为首先对矩阵进行对角化更有效）是：
pow = function(x, n) Reduce(`%*%`, replicate(n, x, simplify = FALSE))

我知道这个问题特别涉及到expm中的一个旧漏洞，但此时此刻它是“矩阵幂 R”搜索结果中的首个结果，因此希望这个简短的方法可以对那些没有安装任何软件包但只是想快速运行矩阵幂的人有所帮助。

2
您可以简单地使用特征值和特征向量来计算矩阵的指数；
# 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)

另外，来自@MichaelChirico的改进答案。矩阵的指数0将返回其单位矩阵，而不是NULL。

pow = function(x, n) {
    if (n == 0) {
        I <- diag(length(diag(x)))
        return(I)        
    } 
    Reduce(`%*%`, replicate(n, x, simplify = FALSE))    
}
3
这个答案可能是正确的，但如果你在解释中加入它，会更加清晰。 - E. Zeytinci
2

不使用任何包的非常快速的解决方案是使用递归： 如果您的矩阵是

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

HTH

1
这并不是非常有用，因为原始的错误已经在两年前修复了... - Ben Bolker
此外，这种方法在进行大指数矩阵幂运算时是一种糟糕的方式。 - m09
2
虽然包中的源代码不可见，因为它被打包在一个.dll文件中，但我相信该包使用的算法是快速幂算法，您可以通过查看名为matpowfast的函数来学习这个算法。
您需要两个变量：
  1. result，用于存储输出结果，
  2. mat，作为中间变量。
要计算A^6，由于6 = 110（二进制写法），最终result = A^6，mat = A^4。对于A^5也是一样的。
当您尝试计算任何8<n<16的A^n时，您可以轻松地检查mat = A^8。如果是这样，您就有了解释。
该包函数使用初始变量A作为中间变量mat。
0

A^5 = (A^4)*A

我猜测这个库会改变原始变量A，使得每一步都涉及将到目前为止的结果与原始矩阵A相乘。你得到的结果似乎很好，只需将它们赋值给一个新变量即可。

计算 A%^%6 同时也会使 A 变为 (初始 A)%^%4。将结果赋值给一个新变量并不能防止我的初始矩阵被改变。 - gd047
听起来你只需要采取不寻常的步骤，首先将矩阵分配给一个新变量。 - John
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接