CODEWITHSUNDEEP
rrcpp
Sosa
Sosa

Reputation: 141

Rcpp Matrix row column permutations

I am trying to mimic the R function that allows to run column and row matrix permutations based on a vector of indices. Like in the following code:

m=matrix(sample(c(0:9),5*5,T),ncol=5,nrow=5)
diag(m)=0
rand=sample(c(1:5))
m[rand,rand]

I tried the following code in c++:

Library(Rcpp)
cppFunction(‘
 NumericMatrix test(NumericMatrix& M, int col, IntegerVector& rand) {
  NumericMatrix M2(col,col);
  for(int a=0;a<col;a++){
    for(int b=a+1;b<col;b++){
       M2(b,a)=M(rand(b),rand(a));
      M2(a,b)=M(rand(a),rand(b));
    }
   }
   return M2;   
}
‘)

But it is very slow: 

microbenchmark::microbenchmark(test(m,5,(rand-1)),m2[rand,rand])

Any ideas how I could speed up the process?

Upvotes: 1

Views: 364

Answers (1)

F. Priv&#233;
F. Priv&#233;

Reputation: 11728

Using a simpler loop:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix test(NumericMatrix& M, int col, IntegerVector& rand) {
  NumericMatrix M2(col,col);
  for(int a=0;a<col;a++){
    for(int b=a+1;b<col;b++){
      M2(b,a)=M(rand(b),rand(a));
      M2(a,b)=M(rand(a),rand(b));
    }
  }
  return M2;   
}

// [[Rcpp::export]]
NumericMatrix test2(const NumericMatrix& M, const IntegerVector& ind) {

  int col = M.ncol();
  NumericMatrix M2(col, col);

  for (int j = 0; j < col; j++)
    for (int i = 0; i < col; i++)
      M2(i, j) = M(ind[i], ind[j]);

  return M2;   
}


/*** R
N <- 500
m <- matrix(sample(c(0:9), N * N, TRUE), ncol = N, nrow = N)
diag(m) <- 0
rand <- sample(N)

all.equal(test(m, ncol(m), rand - 1), m[rand, rand], test2(m, rand - 1))

microbenchmark::microbenchmark(
  test(m, ncol(m), rand - 1),
  m[rand, rand],
  test2(m, rand - 1)
)
*/

For N = 5, the R version is faster, but in terms of nanoseconds.. For example, with N = 500, you get:

Unit: microseconds
                       expr      min       lq     mean   median       uq      max neval
 test(m, ncol(m), rand - 1) 2092.474 2233.020 2843.145 2360.654 2548.050 7412.057   100
              m[rand, rand] 1422.352 1506.117 2064.500 1578.129 1718.345 6700.219   100
         test2(m, rand - 1)  698.595  769.944 1161.747  838.811  928.535 5379.841   100

Upvotes: 2

Related Questions