r - R中矩阵乘法的A ^ k？

假设A是一些方阵。如何在R中轻松对这个矩阵求幂？

我已经尝试了两种方法：带有for-loop hack的Trial 1和Trial 2更优雅，但是与Ak的简单性相去甚远。

试验1

set.seed(10)
t(matrix(rnorm(16),ncol=4,nrow=4)) -> a
for(i in 1:2){a <- a %*% a}

试用2

a <- t(matrix(c(0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0),nrow=4))
i <- diag(4)
(function(n) {if (n<=1) a else (i+a) %*% Recall(n-1)})(10)

最佳答案

尽管Reduce更为优雅，但是for循环解决方案更快，并且似乎与expm ::％^％一样快。

m1 <- matrix(1:9, 3)
m2 <- matrix(1:9, 3)
m3 <- matrix(1:9, 3)
system.time(replicate(1000, Reduce("%*%" , list(m1,m1,m1) ) ) )
#   user  system elapsed
#  0.026   0.000   0.037
mlist <- list(m1,m2,m3)
m0 <- diag(1, nrow=3,ncol=3)
system.time(replicate(1000, for (i in 1:3 ) {m0 <- m0 %*% m1 } ) )
#   user  system elapsed
#  0.013   0.000   0.014

library(expm)  # and I think this may be imported with pkg:Matrix
system.time(replicate(1000, m0%^%3))
# user  system elapsed
#0.011   0.000   0.017

另一方面，matrix.power解决方案要慢得多：

system.time(replicate(1000, matrix.power(m1, 4)) )
   user  system elapsed
  0.677   0.013   1.037

@BenBolker是正确的（再次）。随着指数的增加，for循环在时间上呈线性，而expm ::％^％函数似乎比log（exponent）更好。

> m0 <- diag(1, nrow=3,ncol=3)
> system.time(replicate(1000, for (i in 1:400 ) {m0 <- m0 %*% m1 } ) )
   user  system elapsed
  0.678   0.037   0.708
> system.time(replicate(1000, m0%^%400))
   user  system elapsed
  0.006   0.000   0.006

关于r - R中矩阵乘法的A ^ k？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/9459421/

r中

r - R中矩阵乘法的A ^ k？