R: Fast cartesian product calculation of two numeric matrices

Question

I have two large numeric matrices and want to calculate their cartesian product in R. Is there a way to do it with higher performance and lower memory usage than with my current approach?

EDIT: I added a Rcpp version, which already performs a lot better than my first only-R approach. Since I'm not experienced with Rcpp or RcppArmadillo: Is there a faster/more standardized way to write this Rcpp-function?

m1 <- matrix(sample(0:9, size=110000, replace = TRUE), ncol = 110)
m2 <- matrix(sample(0:9, size=110000, replace = TRUE), ncol = 110)

#Current approach:
m3 <- apply(m1, 1, function(x) x * t(m2))
matrix(m3, ncol = 110, byrow = TRUE)

#EDIT - Rcpp approach
library(Rcpp)
#assuming ncol(m1) == ncol(m2)
cppFunction('IntegerMatrix cartProd(IntegerMatrix m1, IntegerMatrix m2) {
  int nrow1 = m1.nrow(), ncol = m1.ncol(), nrow2 = m2.nrow();
  int orow = 0;
  IntegerMatrix out(nrow1 * nrow2, ncol);

  for (int r1 = 0; r1 < nrow1; r1++) {
    for (int r2 = 0; r2 < nrow2; r2++) {    
      for (int c = 0; c < ncol; c++){
        out(orow, c) = m1(r1, c) * m2(r2, c);
      }
      orow++;
    }
  }
  return out;
}')
m5 <- cartProd(m1, m2)

Show source
| R   | performance   | matrix   | rcpp   | cartesian-product   2017-01-07 19:01 1 Answers

Answers ( 1 )

  1. 2017-01-08 05:01

    The best approach as you have surmised is to use C++ to perform the cartesian product that you desire. Trying to port the code over to Armadillo will yield a minor speed up compared to the pure Rcpp version, which is significantly faster than the written R version. For details on how well each method did, see the benchmark section at the end.


    The first version is almost a direct port into armadillo and actually performs slightly worse than initial pure Rcpp function. The second version uses armadillo's built in submatrix views and each_row() functions to exploit in place evaluation. To achieve parity with the Rcpp version, note the use of the pass-by-reference and the use of a signed integer type yielding const arma::imat&. This avoids a deep copy of the two large integer matrices as the types are matched and a reference is established.

    #include <RcppArmadillo.h>
    // [[Rcpp::depends(RcppArmadillo)]
    
    // --- Version 1
    
    // [[Rcpp::export]]
    arma::imat cartProd_arma(const arma::imat& m1, const arma::imat&  m2) {
      int nrow1 = m1.n_rows, ncol = m1.n_cols, nrow2 = m2.n_rows, orow = 0;
      arma::imat out(nrow1 * nrow2, ncol);
    
      for (int r1 = 0; r1 < nrow1; ++r1) {
        for (int r2 = 0; r2 < nrow2; ++r2) {
          out.row(orow) = m1.row(r1) % m2.row(r2);
          orow++;
        }
      }
      return out;
    }
    
    // --- Version 2
    
    // [[Rcpp::export]]
    arma::imat cartProd_arma2(const arma::imat& m1, const arma::imat&  m2) {
      int nrow1 = m1.n_rows, ncol = m1.n_cols, nrow2 = m2.n_rows, orow = 0;
      arma::imat out(nrow1 * nrow2, ncol);
    
      for (int r1 = 0; r1 < nrow1; ++r1) {
        out.submat(orow, 0, orow + nrow2 - 1, ncol - 1) = m1.row(r1) % m2.each_row();
        orow += nrow2;
      }
      return out;
    }
    

    Quick verification of implementation details matching the initial product

    all.equal( cartProd(m1, m2), cartProd_arma(m1, m2))
    # [1] TRUE
    all.equal( cartProd(m1, m2), cartProd_arma2(m1, m2))
    # [1] TRUE
    

    To generate the benchmarks, I tidied up the initial function slightly by pre-transposing the matrix to avoid multiple transpose calls each time apply was called per row. Furthermore, I've included the function trick showed by @user20650.

    # OP's initial R only solution with slight modifications
    op_R = function(m1, m2){
      m2 <- t(m2)
      m3 <- matrix(apply(m1, 1, function(x) x * m2), ncol = ncol(m1), byrow = TRUE)
    }
    
    # user20650's comment
    so_comment <- function(m1, m2){
      m4 <- matrix(rep(t(m1), each=nrow(m1)) * c(m2), ncol=nrow(m1))
    }
    

    As a result, we have the following microbenchmark

    library("microbenchmark")
    out <- microbenchmark(op_r = op_R(m1, m2), so_comment_r = so_comment(m1, m2), 
                     rcpp = cartProd(m1, m2), arma_v1 = cartProd_arma(m1, m2),
                     arma_v2 = cartProd_arma2(m1, m2),
                     times = 50)
    out
    # Unit: milliseconds
    #         expr       min        lq      mean    median        uq       max neval
    #         op_r 1615.6572 1693.0526 1793.0515 1771.7353 1886.0988 2053.7050    50
    # so_comment_r 2778.0971 2856.6429 2950.5837 2936.7459 3021.4249 3344.4401    50
    #         rcpp  463.6743  482.3118  565.0525  582.1660  614.3714  699.3516    50
    #      arma_v1  597.9004  620.5888  713.4101  726.7572  783.4225  820.3770    50
    #      arma_v2  384.7205  401.9744  490.5118  503.5007  574.6840  622.9876    50
    

    So, from this, we can see that the cartProd_arma2, the submatrix armadillo implementation, is the best function followed closely by cartProd, the pure Rcpp implementation.

◀ Go back