R中矩阵乘法的A ^ k? [英] A^k for matrix multiplication in R?
问题描述
假设 A 是一些正方形矩阵.如何在R中轻松对这个矩阵求幂?
Suppose A is some square matrix. How can I easily exponentiate this matrix in R?
我已经尝试了两种方法:带有for-loop hack的Trial 1和Trial 2更加优雅,但是与 A k 简单性相去甚远.
I tried two ways already: Trial 1 with a for-loop hack and Trial 2 a bit more elegantly but it is still a far cry from Ak simplicity.
试验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 ::%^%
Although Reduce
is more elegant, a for-loop solution is faster and seems to be as fast as 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解决方案要慢得多:
On the other hand the matrix.power solution is much, much slower:
system.time(replicate(1000, matrix.power(m1, 4)) )
user system elapsed
0.677 0.013 1.037
@BenBolker是正确的(再次).随着指数的增加,for循环在时间上呈现线性关系,而expm ::%^%函数似乎甚至比log(exponent)更好.
@BenBolker is correct (yet again). The for-loop appears linear in time as the exponent rises whereas the expm::%^% function appears to be even better than 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中矩阵乘法的A ^ k?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!