一个^ k用于R中的矩阵乘法？

Question

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

我已经尝试了两种方法:试用1使用for-loop hack和试验2更优雅但它仍然与A ^k简单相去甚远.

试验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)

Answer 1

如果A是可对角线的,您可以使用特征值分解:

matrix.power <- function(A, n) {   # only works for diagonalizable matrices
   e <- eigen(A)
   M <- e$vectors   # matrix for changing basis
   d <- e$values    # eigen values
   return(M %*% diag(d^n) %*% solve(M))
}

当A不可对角化时,矩阵M(特征向量矩阵)是单数.因此,使用它A = matrix(c(0,1,0,0),2,2)会给Error in solve.default(M) : system is computationally singular.

Answer 2

该expm套餐有一个%^%运营商:

library("sos")
findFn("{matrix power}")
install.packages("expm")
library("expm")
?matpow
set.seed(10);t(matrix(rnorm(16),ncol=4,nrow=4))->a
a%^%8

Answer 3

虽然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(指数).

> 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