I'm using the Matrix package to create a large (~14000x14000) sparse matrix with a lot of zeros. Does anyone know the best way to calculate the power of this matrix?
I tried A_pow2 = A%^%2 but I get the error: Error in A %^% 2 : not a matrix. Here's a simple example that returns the same error:
A = matrix(3,2,2)
A = Matrix(A,sparse=TRUE)
Apow2 = A%^%2
(edited thanks to @Roland's comments)
A custom function might be able to solve your issue. Per documentation of ?expm::`%^%`
Compute the k-th power of a matrix. Whereas x^k computes element wise powers, x %^% k corresponds to k - 1 matrix multiplications, x %*% x %*% ... %*% x.
We can write a new infix operator to perform the multiplication k-1 times. Not sure how well it will scale, but it works in smaller examples.
> 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