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

Suppose A is some square matrix. How can I easily exponentiate this matrix in R?

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

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 A^k simplicity.

我已经尝试了两种方法：试用1使用for-loop黑客，试用2更优雅，但它仍然与Ak简单相去甚远。

Trial 1

试验1

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

Trial 2

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

6 个解决方案

#1

Although Reduce is more elegant, a for-loop solution is faster and seems to be as fast as expm::%^%

尽管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

On the other hand the matrix.power solution is much, much slower:

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

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

@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).

@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

#2

The expm package has an %^% operator:

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

#3

If A is diagonizable, you could use eigenvalue decomposition:

如果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))
}

When A is not diagonalizable, the matrix M (matrix of eigenvectors) is singular. Thus, using it with A = matrix(c(0,1,0,0),2,2) would give Error in solve.default(M) : system is computationally singular.

当A不可对角化时，矩阵M（特征向量矩阵）是单数。因此，将它与A =矩阵（c（0,1,0,0），2,2）一起使用会在solve.default（M）中给出Error：系统在计算上是单数的。

#4

A shorter solution with eigenvalue decomposition:

具有特征值分解的较短解决方案：

"%^%" <- function(S, power)
         with(eigen(S), vectors %*% (values^power * t(vectors)))

#5

Indeed the expm's package does use exponentiation by squaring.

实际上，expm的包确实通过平方来使用取幂。

In pure r, this can be done rather efficiently like so,

在纯粹的r中，这可以像这样有效地完成，

"%^%" <- function(mat,power){
    base = mat
    out = diag(nrow(mat))
    while(power > 1){
        if(power %% 2 == 1){
            out = out %*% base
        }
        base = base %*% base
        power  = power %/% 2
    }
    out %*% base
}

Timing this,

这个时间，

m0 <- diag(1, nrow=3,ncol=3)
system.time(replicate(10000, m0%^%4000))#expm's %^% function
user  system elapsed 
0.31    0.00    0.31 
system.time(replicate(10000, m0%^%4000))# my %^% function
user  system elapsed 
0.28    0.00    0.28

So, as expected, they are the same speed because they use the same algorithm. It looks like the overhead of the looping r code does not make a significant difference.

因此，正如预期的那样，它们的速度相同，因为它们使用相同的算法。看起来循环r代码的开销并没有产生显着差异。

So, if you don't want to use expm, and need that performance, then you can just use this, if you don't mind looking at imperative code.

因此，如果您不想使用expm并且需要该性能，那么您可以使用它，如果您不介意查看命令式代码。

#6

-1

Here's a solution that's not terribly efficient but works and is easy to understand/implement, works in base, and works almost instantly despite being way less efficient than, e.g., the expm solution:

这是一个效率不高但工作且易于理解/实现的解决方案，在基础上工作，并且几乎立即工作，尽管效率低于例如expm解决方案：

eval(parse(text = paste(rep("A", k), collapse = '%*%')))

e.g.,

例如。，

set.seed(032149)
# force A to be a Markov matrix so it converges
A = matrix(abs(rnorm(100)), 10, 10)
A = A/rowSums(A)
eval(parse(text = paste(rep("A", 1000), collapse = '%*%')))[1:5, 1:5]
#            [,1]       [,2]      [,3]      [,4]       [,5]
# [1,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [2,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [3,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [4,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [5,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486

Again, if efficiency is paramount, the other answers are far superior. But this requires less coding overhead and no packages for situations when you just want to compute some matrix powers in your scratch work.

同样，如果效率至关重要，那么其他答案就更为优越。但是，当您只想在临时工作中计算某些矩阵功率时，这需要较少的编码开销并且不需要打包。

#1