我的示例中为什么Rcpp实现比R函数慢得多？

Question

我有一些C++和R的经验，但对Rcpp还是新手。最近，我在以前的一些项目中使用Rcpp取得了巨大成功，因此决定将其应用到一个新项目中。我很惊讶我的Rcpp代码比相应的R函数慢得多。我尝试简化我的R函数来找出原因，但找不到任何线索。非常欢迎您的帮助和意见！

比较R和Rcpp实现的主要R函数：

main <- function(){

  n <- 50000
  Delta <- exp(rnorm(n))
  delta <- exp(matrix(rnorm(n * 5), nrow = n))
  rx <- matrix(rnorm(n * 20), nrow = n)
  print(microbenchmark(c1 <- test(Delta, delta, rx), times = 500))
  print(microbenchmark(c2 <- rcpp_test(Delta, delta, rx), times = 500))

  identical(c1, c2)
  list(c1 = c1, c2 = c2)
}

R 实现：

test <- function(Delta, delta, rx){

  const <- list()
  for(i in 1:ncol(delta)){
    const[[i]] <- rx * (Delta / (1 + delta[, i]))
  }

  const

}

Rcpp实现：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List rcpp_test(NumericVector Delta, 
               NumericMatrix delta, 
               NumericMatrix rx) {

  int n = Delta.length();
  int m = rx.ncol();

  List c; 
  NumericMatrix c1;
  for(int i = 0; i < delta.ncol(); ++i){
    c1 = NumericMatrix(n, m);
    for(int k = 0; k < n; ++k){
      double tmp = Delta[k] / (1 + delta(k, i));
      for(int j = 0; j < c1.ncol(); ++j){
        c1(k, j) = rx(k, j) * tmp; 
      }
    }
    c.push_back(c1);
  }

  return c;

}

我知道使用Rcpp并不能保证提高效率，但是考虑到我这里所展示的简单示例，我不明白为什么Rcpp代码运行速度如此之慢。

Unit: milliseconds
                         expr      min       lq     mean   median       uq      max neval
 c1 <- test(Delta, delta, rx) 13.16935 14.19951 44.08641 30.43126 73.78581 115.9645   500
Unit: milliseconds
                              expr      min       lq     mean  median       uq      max neval
 c2 <- rcpp_test(Delta, delta, rx) 143.1917 158.7481 171.6116 163.413 173.7677 247.5495   500

理想情况下，在我的项目中，rx是矩阵列表。在for循环中，变量i将用于选择要计算的元素。最初，我怀疑将List传递给Rcpp可能会产生较高的开销，因此在此示例中，我假设rx是一个固定的矩阵，用于所有i。看来这并不是速度缓慢的原因。

- Han Zhang

---

在C++中，push_back()在性能方面的代价可能相当大，因此在需要快速执行速度的应用程序中应该避免使用它。最好预先分配所需的内存。 - RHertel

@RHertel 感谢您的评论。但在这个例子中，像Ralf Stubner所做的那样分配List c并没有什么帮助。请参考我对Ralf Stubner答案的回复。 - Han Zhang

请注意，我并没有将我的评论作为对你的问题的答案，因为我没有说使用 push_back() 是你的代码速度慢的主要原因。我所说的是，push_back() 可以消耗相当多的性能，最好在开始时分配内存。此外，我认为指出在循环中动态增长对象是不良的编程风格是有用的。 - RHertel

Answer 1

您的R代码似乎是最优的，即所有真正的工作都在编译代码中完成。对于C++代码，我发现的主要问题是在紧密循环中调用。如果我将其替换为，则C++解决方案几乎与R一样快。如果我将RcppArmadillo加入混合中，我可以获得非常紧凑的语法，但不比纯Rcpp代码更快。对我来说，这表明很难超越编写良好的R代码：

//  [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::export]]
List arma_test(const arma::vec& Delta,
           const arma::mat& delta,
           const arma::mat& rx) {
  int l = delta.n_cols;
  List c(l);

  for (int i = 0; i < l; ++i) {
    c(i) = rx.each_col() % (Delta / (1 + delta.col(i)));
  }

  return c;  
}

// [[Rcpp::export]]
List rcpp_test(NumericVector Delta, 
               NumericMatrix delta, 
               NumericMatrix rx) {

  int n = Delta.length();
  int m = rx.ncol();

  List c(delta.ncol()); 
  NumericMatrix c1;
  for(int i = 0; i < delta.ncol(); ++i){
    c1 = NumericMatrix(n, m);
    for(int k = 0; k < n; ++k){
      double tmp = Delta[k] / (1 + delta(k, i));
      for(int j = 0; j < m; ++j){
        c1(k, j) = rx(k, j) * tmp; 
      }
    }
    c(i) = c1;
  }

  return c;

}

/*** R
test <- function(Delta, delta, rx){

  const <- list()
  for(i in 1:ncol(delta)){
    const[[i]] <- rx * (Delta / (1 + delta[, i]))
  }

  const

}

n <- 50000
Delta <- exp(rnorm(n))
delta <- exp(matrix(rnorm(n * 5), nrow = n))
rx <- matrix(rnorm(n * 20), nrow = n)
bench::mark(test(Delta, delta, rx),
            arma_test(Delta, delta, rx),
            rcpp_test(Delta, delta, rx))
 */

输出：

# A tibble: 3 x 14
  expression     min    mean  median     max `itr/sec` mem_alloc  n_gc n_itr
  <chr>      <bch:t> <bch:t> <bch:t> <bch:t>     <dbl> <bch:byt> <dbl> <int>
1 test(Delt…  84.3ms  85.2ms  84.9ms  86.6ms     11.7     44.9MB     2     4
2 arma_test… 106.5ms 107.7ms 107.7ms 108.9ms      9.28    38.1MB     3     2
3 rcpp_test… 101.9ms 103.2ms 102.2ms 106.6ms      9.69    38.1MB     1     4
# … with 5 more variables: total_time <bch:tm>, result <list>, memory <list>,
#   time <list>, gc <list>

我还明确地将输出列表初始化为所需大小，避免使用 push_back，但这并没有产生很大的差异。对于类似于 Rcpp 的向量数据结构，您应该绝对避免使用 push_back，因为每次扩展向量时都会进行复制。

- Ralf Stubner

---

避免在for循环中调用c1.ncol()可以大大提高效率！这使得我的Rcpp代码比R代码高效50％。我认为这可能取决于平台（我在Mac上测试过），因为在您的基准测试中，Rcpp函数略慢于R函数。 - Han Zhang

这可能与主题无关，但如果可以的话，我也想在这个帖子中问一下。@RHertel和你都建议避免使用push_back。我的问题是，我们知道List c的大小，但不知道其矩阵元素的大小。在这种情况下，我猜将List初始化为List c(delta.ncol())可能无法分配足够的内存。如果我们知道c的每个元素的大小，有没有更有效的方法来初始化c？实际上，c(i)的维度与rx[i]相同。 - Han Zhang

@HanZhang 很有趣。我刚在另一台机器上重新运行了基准测试（速度更快，但也是Linux），相对结果保持不变。但是“YMMV”是每个基准测试都应该添加的内容...;-) 至于为列表分配内存：当列表初始化时，它包含所需数量的指针（指向 SEXP，R的联合数据类型）。实际数据的内存是单独分配的，但这不应该是问题。 - Ralf Stubner

Answer 2

我想补充一下@RalfStubner的精彩回答。

您会注意到，我们在第一个for循环中进行了许多分配（即）。这可能很昂贵，因为我们除了分配内存外，还将每个元素初始化为0。为了提高效率，我们可以将其更改为以下内容：

NumericMatrix c1 = no_init_matrix(n, m)

我还尽可能添加了关键字const。虽然这样做可能使编译器优化某些代码片段存在争议，但我仍然尽可能添加它来提高代码清晰度（即“我不希望此变量发生更改”）。因此，我们有：

// [[Rcpp::export]]
List rcpp_test_modified(const NumericVector Delta, 
                        const NumericMatrix delta, 
                        const NumericMatrix rx) {

    int n = Delta.length();
    int m = rx.ncol();
    int dCol = delta.ncol();

    List c(dCol);

    for(int i = 0; i < dCol; ++i) {
        NumericMatrix c1 = no_init_matrix(n, m);

        for(int k = 0; k < n; ++k) {
            const double tmp = Delta[k] / (1 + delta(k, i));

            for(int j = 0; j < m; ++j) {
                c1(k, j) = rx(k, j) * tmp; 
            }
        }

        c[i] = c1;
    }

    return c;

}

以下是一些基准测试结果（Armadillo方案未列出）：

bench::mark(test(Delta, delta, rx),
            rcpp_test_modified(Delta, delta, rx),
            rcpp_test(Delta, delta, rx))
# A tibble: 3 x 14
  expression     min   mean  median    max `itr/sec` mem_alloc  n_gc n_itr total_time result memory time 
  <chr>      <bch:t> <bch:> <bch:t> <bch:>     <dbl> <bch:byt> <dbl> <int>   <bch:tm> <list> <list> <lis>
1 test(Delt… 12.27ms 17.2ms 14.56ms 29.5ms      58.1    41.1MB    13     8      138ms <list… <Rpro… <bch…
2 rcpp_test…  7.55ms 11.4ms  8.46ms   26ms      87.8    38.1MB    16    21      239ms <list… <Rpro… <bch…
3 rcpp_test… 10.36ms 15.8ms 13.64ms 28.9ms      63.4    38.1MB    10    17      268ms <list… <Rpro… <bch…
# … with 1 more variable: gc <list>

我们发现使用Rcpp版本后，性能提升了50%。

- Joseph Wood

---

你使用 no_init_matrix 的技巧很棒。在我的应用程序中，矩阵中的许多条目将为零，并且我只会给一些列分配值。尽管如此，我将来仍然会使用这个技巧。 - Han Zhang