Replicate function more efficiently
Suppose I have the following matrices.
x <- matrix(seq(1:4), 2, 2)
y <- matrix(seq(1:4), 2, 2)
I want to do the following.
for(i in 1:5)
{
x <- x %*% y
}
However, this is a simple example. I usually have large matrices for X and Y, and I also have a large number. Thus, using a for-loop can be too time consuming.
Does anyone know whether to use these functions to apply or apply a function.
Thank.
+3
source to share
3 answers
library(expm)
x %*% (y %^% 5)
# [,1] [,2]
#[1,] 5743 12555
#[2,] 8370 18298
Landmarks:
set.seed(42)
x <- matrix(rnorm(1e4), 1e2, 1e2)
y <- matrix(rnorm(1e4), 1e2, 1e2)
fun1 <- function(x, y, j) {
for(i in 1:j)
{
x <- x %*% y
}
x
}
fun2 <- function(x, y, i) {
x %*% (y %^% i)
}
fun3 <- function(x, y, i) {
Reduce("%*%", c(list(x), rep(list(y), i)))
}
library(expm)
all.equal(fun1(x,y,5), fun2(x,y,5))
#[1] TRUE
all.equal(fun1(x,y,5), fun3(x,y,5))
#[1] TRUE
library(microbenchmark)
microbenchmark(fun1(x,y,30),
fun2(x,y,30),
fun3(x,y,30), times=10)
#Unit: milliseconds
# expr min lq median uq max neval
#fun1(x, y, 30) 21.317917 21.908592 22.103380 22.182989 141.933427 10
#fun2(x, y, 30) 5.899368 6.068441 6.235974 6.345301 6.477417 10
#fun3(x, y, 30) 21.385668 21.896274 22.023001 22.086904 22.269527 10
+4
source to share
Just for fun, here's a solution using RcppEigen:
C ++ code:
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using Eigen::Map;
using Eigen::MatrixXd;
typedef Map<MatrixXd> MapMatd;
// [[Rcpp::export]]
NumericMatrix XYpow(NumericMatrix A, NumericMatrix B, const int j) {
const MapMatd X(as<MapMatd>(A)), Y(as<MapMatd>(B));
MatrixXd X1(X);
for (int i = 0; i < j; ++i) X1 = X1 * Y;
return wrap(X1);
}
Then in R:
all.equal(fun2(x,y,5), XYpow(x,y,5)) #[1] TRUE microbenchmark(fun2(x,y,30), XYpow(x,y,30), times=10) #Unit: milliseconds # expr min lq median uq max neval # fun2(x, y, 30) 5.726292 5.768792 5.948027 6.041340 6.276624 10 # XYpow(x, y, 30) 6.926737 7.032061 7.232238 7.512486 7.617502 10
+1
source to share