快速计算两个数值矩阵的笛卡尔积。

Question

我有两个大型数值矩阵，想要在R中计算它们的笛卡尔积。是否有比我当前的方法更高效、内存使用更低的方法？

编辑：我添加了一个Rcpp版本，它的性能已经比我的第一个纯R方法好很多了。由于我对Rcpp或RcppArmadillo并不熟悉：是否有更快/更标准化的方法来编写这个Rcpp函数？

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)

- Daniel_Kuehn

---

这会略微提高一些速度，但是... m4 = matrix(rep(t(m1), each=nrow(m1)) * c(m2), ncol=nrow(m1)) ; all.equal(m3, m4) - user20650

Rcpp是一个选项吗？ - Ben Bolker

@BenBolker：是的，Rcpp也是一个选项。 - Daniel_Kuehn

Answer 1

正如你所推测的那样，最好的方法是使用 C++ 来执行所需的笛卡尔积。尝试将代码移植到 Armadillo 上会比纯 Rcpp 版本略微提速，而后者比编写的 R 版本快得多。有关每种方法的表现如何的详细信息，请参见结尾的基准测试部分。

---

第一个版本几乎是直接移植到armadillo中，并且实际上比最初的纯Rcpp函数性能略差。第二个版本使用了armadillo内置的子矩阵视图和each_row()函数，以利用原地评估。为了实现与Rcpp版本的相同性能，请注意使用传递引用并使用带符号整数类型产生const arma :: imat&。这避免了两个大型integer矩阵的深层复制，因为类型匹配并且建立了一个引用。

#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;
}

---

快速验证实现细节是否与初始产品匹配。

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

---

为了生成基准测试，我稍微整理了初始函数，通过预转置矩阵来避免每次对每行调用apply时进行多次转置调用。此外，我还包括了@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))
}

作为结果，我们有以下微基准测试。

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

因此，我们可以看出，子矩阵的Armadillo实现cartProd_arma2是最好的函数，其次是纯Rcpp实现的cartProd。

- coatless

---

非常好。有趣的是，最好的C++实现仍然不比原始的R实现快超过4倍... - Ben Bolker