R cbind mrix veloce utilizzando Rcpp

Question

cbind in R richiede tempo relativamente a chiamate ripetute, ma è anche potente per vari tipi di dati. Ho scritto codice che è 3 volte più veloce di cbind quando si vincolano due matrici. Ma nel pacchetto dplyr è semplicemente 100 volte più veloce di cbind. È solo un peccato che non possa prendere come input la matrice. Qualcuno può rendere il codice sotto più veloce. Inoltre, come faccio a legare velocemente una matrice sparsa? Ecco il codice che ho usato:

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)) 

Answer 1

Prendendo a prestito un'idea da questo comment da Romain Francois su uno dei miei precedenti avventure 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 

Non un affare enorme, ma non male per un tale cambiamento banale.