在 R 循环中加速 Rcpp 评估
Speed up Rcpp evaluations within R loop
众所周知,Rcpp 中的实现通常比 R 中的实现快得多。我感兴趣的是是否有良好的实践来加速 具有的 Rcpp 函数的单一评估在 R 循环中进行评估.
考虑以下示例,我在 Rcpp 中使用了一个简单的多元正态生成函数:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
using namespace Rcpp;
// [[Rcpp::export]]
mat mvrnormArma(int n, mat sigma) {
int ncols = sigma.n_cols;
mat Y = randn(n, ncols);
return Y * chol(sigma);
}
假设目标是使用以下两个函数生成 10,000 个 10 维多元正态变量:
PureRcpp = function(n){mvrnormArma(n, diag(10))}
LoopRcpp = function(n){for(ii in 1:n){mvrnormArma(1, diag(10))}}
这里,PureRcpp
当然是更可取且更快的解决方案。但是,在某些应用程序中,可能需要在 R 循环中依赖 mvrnormArma
的单个评估。这是 LoopRcpp
中采用的方法,这肯定是较慢的解决方案。然而,当我对这些进行基准测试并看到第二个解决方案实际上有多慢时,我有点惊讶:
> microbenchmark::microbenchmark(PureRcpp(10000), LoopRcpp(10000))
Unit: milliseconds
expr min lq mean median uq max neval cld
PureRcpp(10000) 2.236624 2.365988 2.578869 2.435268 2.565488 10.79609 100 a
LoopRcpp(10000) 52.590143 53.315655 58.080897 55.406020 62.264711 80.96275 100 b
当我们必须在 R 循环中工作时,这种巨大的减速是否只是我们必须忍受的事情,或者是否有可能减少由于循环导致的开销?我知道我们可以用 C++ 重写所有内容,但目标是尽可能提供快速的“'Rcpp within R loop'”解决方案。
正如罗兰指出的那样,这主要是由于函数调用。但是,您可以通过 optimising/adapting 您的代码来节省一些时间(并获得更准确的比较)。
- 通过引用传递给 Cpp 函数
- 不要在循环中创建对角线
- 在单个调度中使用矢量
- 绘制矢量化随机数
// [[Rcpp::export]]
mat draw_randn(int n, int ncols) {
mat Y = randn(n, ncols);
return(Y);
}
// [[Rcpp::export]]
mat mvrnormArma(mat sigma, mat Y) {
return Y * chol(sigma);
}
// [[Rcpp::export]]
mat mvrnormArma_loop(mat& sigma, rowvec& Y) {
return Y * chol(sigma);
}
并对其进行基准测试。
PureRcpp = function(n) {
Y <- draw_randn(n, 10)
I <- diag(10)
mvrnormArma(I, Y)
}
LoopRcpp = function(n) {
Y <- draw_randn(n, 10)
I <- diag(10)
for(ii in 1:n) {mvrnormArma_loop(I, Y[ii, ])}
}
对我来说减少了大约 10 毫秒。
众所周知,Rcpp 中的实现通常比 R 中的实现快得多。我感兴趣的是是否有良好的实践来加速 具有的 Rcpp 函数的单一评估在 R 循环中进行评估.
考虑以下示例,我在 Rcpp 中使用了一个简单的多元正态生成函数:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
using namespace Rcpp;
// [[Rcpp::export]]
mat mvrnormArma(int n, mat sigma) {
int ncols = sigma.n_cols;
mat Y = randn(n, ncols);
return Y * chol(sigma);
}
假设目标是使用以下两个函数生成 10,000 个 10 维多元正态变量:
PureRcpp = function(n){mvrnormArma(n, diag(10))}
LoopRcpp = function(n){for(ii in 1:n){mvrnormArma(1, diag(10))}}
这里,PureRcpp
当然是更可取且更快的解决方案。但是,在某些应用程序中,可能需要在 R 循环中依赖 mvrnormArma
的单个评估。这是 LoopRcpp
中采用的方法,这肯定是较慢的解决方案。然而,当我对这些进行基准测试并看到第二个解决方案实际上有多慢时,我有点惊讶:
> microbenchmark::microbenchmark(PureRcpp(10000), LoopRcpp(10000))
Unit: milliseconds
expr min lq mean median uq max neval cld
PureRcpp(10000) 2.236624 2.365988 2.578869 2.435268 2.565488 10.79609 100 a
LoopRcpp(10000) 52.590143 53.315655 58.080897 55.406020 62.264711 80.96275 100 b
当我们必须在 R 循环中工作时,这种巨大的减速是否只是我们必须忍受的事情,或者是否有可能减少由于循环导致的开销?我知道我们可以用 C++ 重写所有内容,但目标是尽可能提供快速的“'Rcpp within R loop'”解决方案。
正如罗兰指出的那样,这主要是由于函数调用。但是,您可以通过 optimising/adapting 您的代码来节省一些时间(并获得更准确的比较)。
- 通过引用传递给 Cpp 函数
- 不要在循环中创建对角线
- 在单个调度中使用矢量
- 绘制矢量化随机数
// [[Rcpp::export]]
mat draw_randn(int n, int ncols) {
mat Y = randn(n, ncols);
return(Y);
}
// [[Rcpp::export]]
mat mvrnormArma(mat sigma, mat Y) {
return Y * chol(sigma);
}
// [[Rcpp::export]]
mat mvrnormArma_loop(mat& sigma, rowvec& Y) {
return Y * chol(sigma);
}
并对其进行基准测试。
PureRcpp = function(n) {
Y <- draw_randn(n, 10)
I <- diag(10)
mvrnormArma(I, Y)
}
LoopRcpp = function(n) {
Y <- draw_randn(n, 10)
I <- diag(10)
for(ii in 1:n) {mvrnormArma_loop(I, Y[ii, ])}
}
对我来说减少了大约 10 毫秒。