R中的矩阵幂

我已经修复了R-forge源代码(expm包)中的错误, svn版本。53. --&燃气轮机; expm R-forge page 解决您的问题(但应在24小时内):

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

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

感谢“gd047”通过电子邮件提醒我这个问题。 注意,R-forge也有自己的bug跟踪工具。
马丁特

首先,请注意,仅将矩阵分配给一个新变量并没有帮助:

> 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代码中,有以下注释:

尽管源代码在包中不可见,因为它是打包在 .dll file ,我相信软件包使用的算法是 fast exponentiation algorithm matpowfast 相反。

1. result ,以便存储输出,
2. mat ,作为中间变量。

计算 A^6 6 = 110 (二进制写入),最后, result = A^6 mat = A^4 . 这对我来说是一样的 A^5 .

你可以很容易地检查 mat = A^8 A^n 对于任何 8<n<16 . 如果是,你有你的解释。

package函数使用初始变量 A 作为中间变量 垫子

非常快速的解决方案 正在使用递归:

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

一个低效的版本(因为首先对角化矩阵更有效) base 不费吹灰之力:

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

我知道这个问题是关于 expm ,但这是目前“矩阵幂R”的第一个结果之一,所以希望这个小速记可以对其他人有用,他们最终只是在这里寻找一个快速的方法来运行矩阵幂而不安装任何包。

我假设库对原始变量A进行了变异,因此每一步都需要将结果与原始矩阵A相乘。返回的结果看起来很好,只需将它们赋给一个新变量。

你可以简单地用特征值和特征向量来计算矩阵的指数;

# 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 . 这个 exponent 0 NULL .

pow = function(x, n) {
    if (n == 0) {
        I <- diag(length(diag(x)))
        return(I)        
    } 
    Reduce(`%*%`, replicate(n, x, simplify = FALSE))    
}