比较非连续子矩阵选择的 RcppArmadillo 和 R 运行速度

Question

我的主要目标是使用两组用于行和列的二进制向量来选择不连续的子矩阵。这是我使用 Rcpp、RcppArmadillo 和 RcppEigen 在 C++ 中实现的 MCMC 循环需要执行的众多步骤之一。

三种可能的方法是 (1) 使用 RcppArmadillo，(2) 从 Rcpp 调用我的 R 函数，以及 (3) 直接使用 R 并将结果传递给 C++。虽然最后一个选项对我来说一点都不方便。

然后我比较了这三种场景的性能速度。有趣的是，直接的 R 代码比其他两个要快得多！更令我惊讶的是，当我从 Rcpp 调用确切的 R 函数时，它比我直接从 R 调用它时慢得多。我希望那些运行速度与older post 中的示例中建议的运行速度相对相同.

无论如何，计时结果对我来说似乎有点奇怪。有什么原因吗？我使用带有 El Capitan OS、2.5 Ghz Intel Core i7 的 Macbook Pro。会不会与我的系统、Mac OSX 或我的机器上安装 Rcpp 的方式有关？

提前致谢！

代码如下：

CPP 部分：

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// (1) Using RcppArmadillo functions:
// [[Rcpp::export]]
mat subselect(NumericMatrix X, uvec rows, uvec cols){

  mat XX(X.begin(), X.nrow(),X.ncol(), false);
  mat y = XX.submat(find(rows>0),find(cols>0));
  return (y);
}

// (2) Calling the function from R:
// [[Rcpp::export]]
NumericalMatrix subselect2(NumericMatrix X, NumericVector rows, NumericVector cols){ 

  Environment stats;
  Function submat = stats["submat"];
  NumericMatrix outmat=submat(X,rows,cols);
  return(wrap(outmat));
}

R 部分：

library(microbenchmark)

# (3) My R function:
submat <- function(mat,rvec,cvec){
 return(mat[as.logical(rvec),as.logical(cvec)])
}

# Comparing the performances:

// Generating data:
set.seed(432)
rows <- rbinom(1000,1,0.1)
cols <- rbinom(1000,1,0.1)
amat <- matrix(1:1e06,1000,1000)

//benchmarking:
microbenchmark(subselect(amat,rows,cols),
           subselect2(amat,rows,cols),
           submat(amat,rows,cols))

结果：

                         expr     min       lq     mean    median       uq        max   neval   

  subselect(amat, rows, cols) 893.670 1566.882 2297.991  1675.282 2184.783   8462.142     100
 subselect2(amat, rows, cols) 928.418 1581.553 3554.805  1657.454 2060.837 138801.050     100
     submat(amat, rows, cols)  36.313   55.748   66.782    62.709   73.975    136.970     100

Answer 1

这里有几件事值得一提。首先，您在基准测试设计中犯了一个微妙的错误，该错误对您的 Armadillo 函数 subselect 的性能产生了重大影响。观察：

set.seed(432)
rows <- rbinom(1000, 1, 0.1)
cols <- rbinom(1000, 1, 0.1)

imat <- matrix(1:1e6, 1000, 1000)
nmat <- imat + 0.0

storage.mode(imat)
# [1] "integer"

storage.mode(nmat)
# [1] "double"

microbenchmark(
    "imat" = subselect(imat, rows, cols),
    "nmat" = subselect(nmat, rows, cols)
)
# Unit: microseconds
#  expr      min       lq      mean    median        uq       max neval
#  imat 3088.140 3218.013 4355.2956 3404.4685 4585.1095 21662.540   100
#  nmat  139.298  167.116  223.2271  209.4585  238.6875   533.035   100

---

尽管 R 经常将整数文字（例如 1、2、3、...）视为浮点值，但序列运算符 : 是少数例外之一，

storage.mode(c(1, 2, 3, 4, 5))
# [1] "double"

storage.mode(1:5)
# [1] "integer"

这就是表达式matrix(1:1e6, 1000, 1000) 返回integer 矩阵而不是numeric 矩阵的原因。这是有问题的，因为subselect 期望的是NumericMatrix，而不是IntegerMatrix，并且传递后一种类型会触发深层复制，因此上述基准的差异超过一个数量级。

---

其次，R 函数 submat 和 C++ 函数 subselect 在二进制索引向量分布上的相对性能存在显着差异，这可能是由于底层算法的差异。对于更稀疏的索引（0 比 1 的比例更大），R 函数胜出；而对于更密集的索引，则相反。这似乎也是矩阵大小（或可能只是维度）的函数，如下图所示，其中行和列索引向量是使用 rbinom 生成的，成功参数为 0.0、0.05、0.10、...、0.95 , 1.0 -- 首先使用 1e3 x 1e3 矩阵，然后使用 1e3 x 1e4 矩阵。最后包含此代码。

---

---

基准代码：

library(data.table)
library(microbenchmark)
library(ggplot2)

test_data <- function(nr, nc, p, seed = 123) {
    set.seed(seed)
    list(
        x = matrix(rnorm(nr * nc), nr, nc),
        rv = rbinom(nr, 1, p),
        cv = rbinom(nc, 1, p)
    )
}

tests <- lapply(seq(0, 1, 0.05), function(p) {
    lst <- test_data(1e3, 1e3, p)
    list(
        p = p,
        benchmark = microbenchmark::microbenchmark(
            R = submat(lst[[1]], lst[[2]], lst[[3]]),
            Arma = subselect(lst[[1]], lst[[2]], lst[[3]])
        )
    )
})

gt <- rbindlist(
    Map(function(g) {
        data.table(g[[2]])[
            ,.(Median.us = median(time / 1000)), 
            by = .(Expr = expr)
        ][order(Median.us)][
            ,Relative := Median.us / min(Median.us)
        ][,pSuccess := sprintf("%3.2f", g[[1]])]
    }, tests)
)

ggplot(gt) +
    geom_point(
        aes(
            x = pSuccess, 
            y = Relative, 
            color = Expr
        ),
        size = 2,
        alpha = 0.75
    ) +
    theme_bw() +
    ggtitle("1e3 x 1e3 Matrix")

## change `test_data(1e3, 1e3, p)` to
## `test_data(1e3, 1e4, p)` inside of 
## `tests <- lapply(...) ...` to generate 
## the second plot