return Eigen::VectorXi 的最有效方法，R 的元素超过 2^31-1

Question

我有一个 Eigen::VectorXi 类型的向量 x，其中包含超过 2^31-1 个条目，我想 return 到 R。我可以通过复制来做到这一点x 明智地进入类型为 Rcpp::IntegerVector 的新向量，但这似乎很慢。

我在想：

是否有更有效的解决方法；
为什么在下面的可重现示例中 Rcpp::wrap(x) 不起作用。

test.cpp

#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
SEXP foo(const R_xlen_t size) {
  Eigen::VectorXi x(size);
  for (R_xlen_t i = 0; i < size; i++) {
    x(i) = 1;
  }
  return Rcpp::wrap(x);
}
    
// [[Rcpp::export]]
Rcpp::IntegerVector fooSlow(const R_xlen_t size) {
  Eigen::VectorXi x(size);
  for (R_xlen_t i = 0; i < size; i++) {
    x(i) = 1;
  }
  Rcpp::IntegerVector y(size);
  for (R_xlen_t i = 0; i < size; i++) {
    y(i) = x(i);
  }   
  return y;
}

test.R

Rcpp::sourceCpp("./test.cpp")
a <- foo(2^(31)) # Error in foo(2^(31)) : negative length vectors are not allowed
a <- fooSlow(2^(31)) # runs fine but it's slow

Answer 1

Rcpp::wrap 正在调度 RcppEigen 中实现的特征矩阵和向量的方法。该方法目前似乎不支持 long vectors。（编辑：现在可以；见下文。）

有关负长度的错误由 allocVector3 here. It arises when allocVector3 is called with a negative value for its argument length. My guess is that Rcpp::wrap tries to represent 2^31 as an int, resulting in integer overflow. Maybe this happens here?

抛出

无论如何，您似乎偶然发现了一个错误，因此您可以考虑在 GitHub. (Edit: Never mind - I've just submitted a patch.) (Edit: Patched now, if you'd like to build RcppEigen from sources [commit 5fd125e 或更高版本] 上与 RcppEigen 维护者分享您的示例，以便更新您的 Rcpp::wrap.)

为了回答您的第一个问题，我根据 std::memcpy 将您的两种方法与我自己的方法进行了比较。 std::memcpy 方法支持长向量，仅比 Rcpp::wrap.

稍慢

`std::memcpy` 方法

Eigen::VectorXi x 和 Rcpp::IntegerVector y 下面的 C 数组具有相同的类型 (int) 和长度 (n)，因此它们包含相同的字节数。您可以使用 std::memcpy 将那个字节数从一个内存地址复制到另一个内存地址，而无需 for 循环。困难的部分是知道如何获取地址。 Eigen::VectorXi 有一个成员函数 data 即 returns 底层 int 数组的地址。整数类型的 R 对象使用 R API 中的 INTEGER，做同样的事情。

测试

Rcpp::sourceCpp(code = '
#include <RcppEigen.h>
#include <Rinternals.h>

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Rcpp::IntegerVector f_for(const R_xlen_t n) {
  Eigen::VectorXi x(n);
  for (R_xlen_t i = 0; i < n; ++i) {
    x(i) = i % 10;
  }
  Rcpp::IntegerVector y(n);
  for (R_xlen_t i = 0; i < n; ++i) {
    y(i) = x(i);
  }
  return y;
}

// [[Rcpp::export]]
Rcpp::IntegerVector f_wrap(const R_xlen_t n) {
  Eigen::VectorXi x(n);
  for (R_xlen_t i = 0; i < n; ++i) {
    x(i) = i % 10;
  }
  return Rcpp::wrap(x);
}

// [[Rcpp::export]]
Rcpp::IntegerVector f_memcpy(const R_xlen_t n) {
  Eigen::VectorXi x(n);
  for (R_xlen_t i = 0; i < n; ++i) {
    x(i) = i % 10;
  }
  Rcpp::IntegerVector y(n);
  memcpy(INTEGER(y), x.data(), n * sizeof(int));
  return y;
}
')

n <- 100L
x <- rep_len(0:9, n)
identical(f_for(n),    x) # TRUE
identical(f_wrap(n),   x) # TRUE
identical(f_memcpy(n), x) # TRUE

b <- function(n) microbenchmark::microbenchmark(f_for(n), f_wrap(n), f_memcpy(n), setup = gc(FALSE))

b(2^10)
## Unit: microseconds
##         expr   min     lq     mean median      uq     max neval
##     f_for(n) 6.806 8.5280 15.09497 10.332 11.8900 461.496   100
##    f_wrap(n) 4.469 6.2115 12.60750  8.569  9.7170 435.420   100
##  f_memcpy(n) 4.633 7.0520 13.64193  9.061  9.6965 465.924   100

b(2^20)
## Unit: microseconds
##         expr      min        lq      mean    median       uq      max neval
##     f_for(n) 3094.106 3118.2960 3160.2501 3132.4205 3171.329 3515.996   100
##    f_wrap(n)  864.690  890.0485  912.7006  905.4440  929.593  988.797   100
##  f_memcpy(n)  940.048  971.6590 1001.9805  987.3825 1009.195 1428.235   100

b(2^30)
## Unit: seconds
##         expr      min       lq     mean   median       uq      max neval
##     f_for(n) 3.527164 3.554672 3.575698 3.573021 3.593006 3.693711   100
##    f_wrap(n) 1.119750 1.133130 1.143425 1.139702 1.149030 1.203602   100
##  f_memcpy(n) 1.304877 1.330994 1.343253 1.339099 1.354286 1.422912   100

return Eigen::VectorXi 的最有效方法，R 的元素超过 2^31-1

Most efficient way to return Eigen::VectorXi with more than 2^31-1 elements to R

r

rcpp

rcppeigen

`std::memcpy` 方法

测试

return Eigen::VectorXi 的最有效方法，R 的元素超过 2^31-1

Most efficient way to return Eigen::VectorXi with more than 2^31-1 elements to R

r

rcpp

rcppeigen

std::memcpy 方法

测试

`std::memcpy` 方法