使用Rcpp快速合并矩阵

Question

cbind在R中重复调用时时间相对较长，但对于各种数据类型也很有用。我编写了一段代码，将两个矩阵绑定时比cbind快3倍。但是，在dplyr包中的bind_cols仅比cbind快100倍。可惜它不能接受矩阵作为输入。请问有人能让下面的代码更快吗？另外，我该如何快速绑定稀疏矩阵？以下是我使用的代码：

require( Rcpp )

func <- 'NumericMatrix mmult(NumericMatrix a,NumericMatrix b) {
    //the colnumber of first matrix
    int acoln=a.ncol();
    //the colnumber of second matrix
    int bcoln=b.ncol();
    //build a new matrix, the dim is a.nrow() and acoln+bcoln
    NumericMatrix out(a.nrow(),acoln+bcoln) ;
    for (int j = 0; j < acoln + bcoln; j++) {
        if (j < acoln) {
            out(_,j) = a(_,j);
        } else {
            //put the context in the second matrix to the new matrix
            out(_,j) = b(_,j-acoln);
        }
    }
    return out ;
}'

a <- matrix(rep(1,2000*100),2000)
b <- matrix(rep(2,2000*10),2000)

cppFunction(func)

system.time(for (i in seq(1,800)) {mmult(a,b)})
system.time(for (i in seq(1,800)) {cbind(a,b)})
identical(mmult(a,b),cbind(a,b))

- Weve Wang

Answer 1

借鉴Romain Francois在我之前的Rcpp探险中提出的这个评论的想法，

func1 <- 'NumericMatrix mmult1(NumericMatrix a, NumericMatrix b) {
  int acoln = a.ncol();
  int bcoln = b.ncol();
  NumericMatrix out = no_init_matrix(a.nrow(), acoln + bcoln);
  for (int j = 0; j < acoln + bcoln; j++) {
    if (j < acoln) {
      out(_, j) = a(_, j);
    } else {
      out(_, j) = b(_, j - acoln);
    }
  }
  return out;
}'

cppFunction(func1)
set.seed(42)
a <- matrix(rnorm(1e7), 1e3)
b <- matrix(runif(1e7), 1e3)

identical(mmult(a, b), mmult1(a, b))
#TRUE

library(microbenchmark)
microbenchmark(mmult(a, b), 
               mmult1(a, b), 
               cbind(a, b),
               times = 10)

#Unit: milliseconds
#         expr    min     lq   mean median    uq   max neval
#  mmult(a, b)  69.64  70.52  89.71  72.28 128.8 136.6    10
# mmult1(a, b)  50.84  50.95  69.65  51.43 111.6 114.4    10
#  cbind(a, b) 192.35 194.67 201.13 195.30 196.1 255.9    10

这并不是很重要，但对于这样一个微不足道的变化来说也不错。

- tonytonov