我正在使用 Matrix 包创建一个包含大量零的大型 (~14000x14000) 稀疏矩阵。有谁知道计算这个矩阵的幂的最佳方法？
我试过 A_pow2 = A%^%2 但我得到错误： A %^% 2 中的错误：不是矩阵。这是一个返回相同错误的简单示例：
A = matrix(3,2,2)
A = Matrix(A,sparse=TRUE)
Apow2 = A%^%2
【问题讨论】：
Matrix pkg 的%^% 运算符吗？您是否尝试过简单地使用^ 代替%^%？
%^% 运算符。只是将它定位在expm pkg
标签： r sparse-matrix
（感谢@Roland 的 cmets 编辑）
自定义函数或许可以解决您的问题。根据?expm::`%^%`的文档
计算矩阵的 k 次方。而 x^k 计算元素明智 幂，x %^% k 对应于 k - 1 次矩阵乘法，x %*% x %*% ... %*% x。
我们可以编写一个新的中缀运算符来执行 k-1 次乘法。不确定它的扩展性如何，但它适用于较小的示例。
> library(Matrix)
> library(expm)
> A = matrix(3,2,2)
> B = Matrix(A,sparse=TRUE)
>
> # changed lapply to rep list
> `%^^%` = function(x, k) Reduce(`%*%`, rep(list(x), k))
> # per Roland for loop approach will be better on memory
> `%^^%` = function(x, k) {for (i in 1:(k - 1)) x <- x %*% x; x}
>
> as.matrix(B%^^%2)
[,1] [,2]
[1,] 18 18
[2,] 18 18
> A%^%2
[,1] [,2]
[1,] 18 18
[2,] 18 18
【讨论】：
`%^^%` = function(x, k) Reduce(`%*%`, rep(list(x), k))
for (i in 1:(k - 1)) x <- x %*% x（可能会通过 JIT 编译进行优化）。