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

Question

假设A是某个方阵。如何在 R 中轻松地对该矩阵求幂？

我已经尝试了两种方法:试用 1 使用 for 循环 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) 一起使用会给出 solve.default(M) 中的错误:系统在计算上是奇异的.