试图在 R 中计算矩阵的幂，我发现包 expm 实现了运算符 %^%。
所以 x %^% k 计算矩阵的 k 次方。
> A<-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3)
> A %^% 5
[,1] [,2] [,3]
[1,] 6469 18038 2929
[2,] 21837 60902 9889
[3,] 10440 29116 4729
但是，令我惊讶的是：
> A
[,1] [,2] [,3]
[1,] 691 1926 312
[2,] 2331 6502 1056
[3,] 1116 3108 505
不知何故，初始矩阵 A 已更改为 A %^% 4 !!!
如何进行矩阵幂运算？
【问题讨论】：
标签： r matrix linear-algebra
A^5 = (A^4)*A
我想库改变了原始变量 A，因此每个步骤都涉及将结果与原始矩阵 A 相乘。你得到的结果看起来很好，只需将它们分配给一个新的变量。
【讨论】：
虽然源代码在包中不可见，因为它被打包在.dll file中，但我相信包使用的算法是fast exponentiation algorithm，您可以通过查看名为@987654323的函数来学习@ 代替。
你需要两个变量：
result，为了存储输出，mat，作为中间变量。要计算A^6，因为6 = 110（二进制写入），最后是result = A^6和mat = A^4。 A^5 也是如此。
当您尝试计算任何8<n<16 的A^n 时，您可以轻松检查mat = A^8 是否存在。如果是这样，你有你的解释。
封装函数使用初始变量A作为中间变量mat。
【讨论】：
这不是一个正确的答案，但可能是进行此讨论并了解 R 内部工作原理的好地方。这种错误之前在我使用的另一个包中已经悄悄出现。
首先，请注意，简单地将矩阵首先分配给一个新变量并没有帮助：
> A <- B <-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3)
> r1 <- A %^% 5
> A
[,1] [,2] [,3]
[1,] 691 1926 312
[2,] 2331 6502 1056
[3,] 1116 3108 505
> B
[,1] [,2] [,3]
[1,] 691 1926 312
[2,] 2331 6502 1056
[3,] 1116 3108 505
我的猜测是，R 正试图通过引用而不是值进行智能传递。要真正让它发挥作用，您需要做一些事情来区分 A 和 B：
`%m%` <- function(x, k) {
tmp <- x*1
res <- tmp%^%k
res
}
> B <-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3)
> r2 <- B %m% 5
> B
[,1] [,2] [,3]
[1,] 1 2 1
[2,] 3 8 1
[3,] 0 4 1
这样做的显式方法是什么？
最后，在包的 C 代码中，有这样的注释：
但我不明白为什么 R 会让 C/Fortran 代码在全局环境中产生副作用。
【讨论】：
我已在 R-forge 源代码（“expm”包）中修复了该错误， svn 转。 53. --> expm R-forge page 由于某种原因，网页仍然显示 rev.52，所以以下可能还没有 解决您的问题（但应在 24 小时内）：
install.packages("expm", repos="http://R-Forge.R-project.org")
否则直接获取svn版本，自己安装：
svn checkout svn://svn.r-forge.r-project.org/svnroot/expm
感谢“gd047”通过电子邮件提醒我这个问题。
请注意，R-forge 也有自己的错误跟踪工具。
马丁
【讨论】：
非常快速的解决方案不使用任何包正在使用递归： 如果你的矩阵是一个
powA = function(n)
{
if (n==1) return (a)
if (n==2) return (a%*%a)
if (n>2) return ( a%*%powA(n-1))
}
HTH
【讨论】：
在base 中不费吹灰之力的低效版本（因为首先对矩阵进行对角化更有效）是：
pow = function(x, n) Reduce(`%*%`, replicate(n, x, simplify = FALSE))
我知道这个问题专门针对 expm 中的一个旧错误，但它是目前“矩阵幂 R”的第一个结果，所以希望这个简短的速记对最终在这里的其他人有用只是在寻找一种无需安装任何软件包即可运行矩阵电源的快速方法。
【讨论】：
您可以简单地使用特征值和特征向量来计算矩阵的指数；
# for a given matrix, A of power n
eig_vectors <- eigen(A)$vectors
eig_values <- eigen(A)$values
eig_vectors %*% diag(eig_values)^n %*% solve(eig_vectors)
或者来自@MichaelChirico 的改进答案。矩阵的exponent 0 将返回其单位矩阵而不是NULL。
pow = function(x, n) {
if (n == 0) {
I <- diag(length(diag(x)))
return(I)
}
Reduce(`%*%`, replicate(n, x, simplify = FALSE))
}
【讨论】：