为了进行C++/Rcpp编程,我实现了一个(样本)标准偏差函数:

#include <Rcpp.h>
#include <vector>
#include <cmath>
#include <numeric>

// [[Rcpp::export]]
double cppSD(Rcpp::NumericVector rinVec)
{
  std::vector<double> inVec(rinVec.begin(),rinVec.end());
  int n = inVec.size();
  double sum = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  double mean = sum / inVec.size();

  for(std::vector<double>::iterator iter = inVec.begin();
      iter != inVec.end(); ++iter){
        double temp;
        temp= (*iter - mean)*(*iter - mean);
        *iter = temp;
      }

  double sd = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  return std::sqrt( sd / (n-1) );
}

考虑到可以在 vector 上调用它,我还决定从Armadillo库中测试stddev函数:
#include <RcppArmadillo.h>

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

using namespace Rcpp;

// [[Rcpp::export]]

double armaSD(arma::colvec inVec)
{
  return arma::stddev(inVec);
}

然后,我针对基本R函数sd()对这两个函数进行了基准测试,以获取一些大小不同的 vector :
Rcpp::sourceCpp('G:/CPP/armaSD.cpp')
Rcpp::sourceCpp('G:/CPP/cppSD.cpp')
require(microbenchmark)
##
## sample size = 1,000: armaSD() < cppSD() < sd()
X <- rexp(1000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
#      expr    min     lq median     uq    max neval
# armaSD(X)  4.181  4.562  4.942  5.322 12.924   100
#     sd(X) 17.865 19.766 20.526 21.287 86.285   100
#  cppSD(X)  4.561  4.941  5.321  5.701 29.269   100
##
## sample size = 10,000: armaSD() < cppSD() < sd()
X <- rexp(10000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
#      expr    min     lq  median      uq     max neval
# armaSD(X) 24.707 25.847 26.4175 29.6490  52.455   100
#     sd(X) 51.315 54.356 55.8760 61.1980 100.730   100
#  cppSD(X) 26.608 28.128 28.8885 31.7395 114.413   100
##
## sample size = 25,000: armaSD() < cppSD() < sd()
X <- rexp(25000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
#      expr     min       lq  median      uq     max neval
# armaSD(X)  66.900  67.6600  68.040  76.403 155.845   100
#     sd(X) 108.332 111.5625 122.016 125.817 169.910   100
#  cppSD(X)  70.320  71.0805  74.692  80.203 102.250   100
##
## sample size = 50,000: cppSD() < sd() < armaSD()
X <- rexp(50000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
#      expr     min       lq   median      uq     max neval
# armaSD(X) 249.733 267.4085 297.8175 337.729 642.388   100
#     sd(X) 203.740 229.3975 240.2300 260.186 303.709   100
#  cppSD(X) 162.308 185.1140 239.6600 256.575 290.405   100
##
## sample size = 75,000: sd() < cppSD() < armaSD()
X <- rexp(75000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
#      expr     min       lq   median       uq     max neval
# armaSD(X) 445.110 479.8900 502.5070 520.5625 642.388   100
#     sd(X) 310.931 334.8780 354.0735 379.7310 429.146   100
#  cppSD(X) 346.661 380.8715 400.6370 424.0140 501.747   100

对于较小的样本,我的C++函数cppSD()stats::sd()快,但由于stats::sd()是矢量化的,所以对于较大尺寸的 vector ,它的速度却慢得多,这让我并不感到惊讶。但是,我没想到arma::stddev()函数会导致这种性能下降,因为它似乎也以 vector 方式运行。我使用arma::stdev()的方式是否存在问题,或者仅仅是以可以更有效地处理较大 vector 的方式(以我假设的stats::sd()编写)C?任何输入将不胜感激。

更新:

尽管我的问题最初是关于正确使用arma::stddev的问题,而不是关于尝试找到可能的最有效方法来计算样本标准差的问题,但有趣的是Rcpp::sd sugar函数的性能如此出色。为了使事情变得更有趣,我将下面的arma::stddevRcpp::sd函数与我从JJ Allaire的Rcpp Gallery文章中的两篇改编的RcppParallel版本进行了基准测试-herehere:
library(microbenchmark)
set.seed(123)
x <- rnorm(5.5e06)
##
Res <- microbenchmark(
  armaSD(x),
  par_sd(x),
  sd_sugar(x),
  times=500L,
  control=list(warmup=25))
##
R> print(Res)
Unit: milliseconds
        expr       min        lq      mean    median        uq       max neval
   armaSD(x) 24.486943 24.960966 26.994684 25.255584 25.874139 123.55804   500
   par_sd(x)  8.130751  8.322682  9.136323  8.429887  8.624072  22.77712   500
 sd_sugar(x) 13.713366 13.984638 14.628911 14.156142 14.401138  32.81684   500

这是在运行64位linux,i5-4200U CPU @ 1.60GHz处理器的笔记本电脑上;但是我猜想par_sdsugar_sd之间的差异在Windows机器上不会那么重要。

以及RcppParallel版本的代码(该代码更长一些,并且需要C++ 11兼容的编译器来用于operator()函数的重载InnerProduct中使用的lambda表达式):
#include <Rcpp.h>
#include <RcppParallel.h>
// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::plugins(cpp11)]]

/*
 * based on: http://gallery.rcpp.org/articles/parallel-vector-sum/
 */
struct Sum : public RcppParallel::Worker {

  const RcppParallel::RVector<double> input;
  double value;

  Sum(const Rcpp::NumericVector input)
  : input(input), value(0) {}
  Sum(const Sum& sum, RcppParallel::Split)
  : input(sum.input), value(0) {}

  void operator()(std::size_t begin, std::size_t end) {
    value += std::accumulate(input.begin() + begin,
                             input.begin() + end,
                             0.0);
  }

  void join(const Sum& rhs) {
    value += rhs.value;
  }

};

/*
 * based on: http://gallery.rcpp.org/articles/parallel-inner-product/
 */
struct InnerProduct : public RcppParallel::Worker {

  const RcppParallel::RVector<double> x;
  const RcppParallel::RVector<double> y;
  double mean;
  double product;

  InnerProduct(const Rcpp::NumericVector x,
               const Rcpp::NumericVector y,
               const double mean)
  : x(x), y(y), mean(mean), product(0) {}

  InnerProduct(const InnerProduct& innerProduct,
               RcppParallel::Split)
  : x(innerProduct.x), y(innerProduct.y),
    mean(innerProduct.mean), product(0) {}

  void operator()(std::size_t begin, std::size_t end) {
    product += std::inner_product(x.begin() + begin,
                                  x.begin() + end,
                                  y.begin() + begin,
                                  0.0, std::plus<double>(),
                [&](double lhs, double rhs)->double {
                  return ( (lhs-mean)*(rhs-mean) );
                });
  }

  void join(const InnerProduct& rhs) {
    product += rhs.product;
  }

};

// [[Rcpp::export]]
double par_sd(const Rcpp::NumericVector& x_)
{
  int N = x_.size();
  Rcpp::NumericVector y_(x_);
  Sum sum(x_);
  RcppParallel::parallelReduce(0, x_.length(), sum);

  double mean = sum.value / N;
  InnerProduct innerProduct(x_, y_, mean);
  RcppParallel::parallelReduce(0, x_.length(), innerProduct);

  return std::sqrt( innerProduct.product / (N-1) );
}

最佳答案

您在实例化Armadillo对象方面犯了一个细微的错误-这导致复制并因此降低了性能。

改用const arma::colvec & invec接口(interface),一切都很好:

R> sourceCpp("/tmp/sd.cpp")

R> library(microbenchmark)

R> X <- rexp(500)

R> microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
Unit: microseconds
       expr    min      lq  median      uq    max neval
  armaSD(X)  3.745  4.0280  4.2055  4.5510 19.375   100
 armaSD2(X)  3.305  3.4925  3.6400  3.9525  5.154   100
      sd(X) 22.463 23.6985 25.1525 26.0055 52.457   100
   cppSD(X)  3.640  3.9495  4.2030  4.8620 13.609   100

R> X <- rexp(5000)

R> microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
Unit: microseconds
       expr    min      lq  median      uq    max neval
  armaSD(X) 18.627 18.9120 19.3245 20.2150 34.684   100
 armaSD2(X) 14.583 14.9020 15.1675 15.5775 22.527   100
      sd(X) 54.507 58.8315 59.8615 60.4250 84.857   100
   cppSD(X) 18.585 19.0290 19.3970 20.5160 22.174   100

R> X <- rexp(50000)

R> microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
Unit: microseconds
       expr     min      lq  median      uq     max neval
  armaSD(X) 186.307 187.180 188.575 191.825 405.775   100
 armaSD2(X) 142.447 142.793 143.207 144.233 155.770   100
      sd(X) 382.857 384.704 385.223 386.075 405.713   100
   cppSD(X) 181.601 181.895 182.279 183.350 194.588   100
R>

这是基于我的代码版本的,其中所有内容都是一个文件,并且按照我的建议定义了armaSD2 -从而赢得了成功。
#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
#include <vector>
#include <cmath>
#include <numeric>

// [[Rcpp::export]]
double cppSD(Rcpp::NumericVector rinVec) {
  std::vector<double> inVec(rinVec.begin(),rinVec.end());
  int n = inVec.size();
  double sum = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  double mean = sum / inVec.size();

  for(std::vector<double>::iterator iter = inVec.begin();
      iter != inVec.end();
      ++iter){
    double temp = (*iter - mean)*(*iter - mean);
    *iter = temp;
  }

  double sd = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  return std::sqrt( sd / (n-1) );
}

// [[Rcpp::export]]
double armaSD(arma::colvec inVec) {
  return arma::stddev(inVec);
}

//  [[Rcpp::export]]
double armaSD2(const arma::colvec & inVec) { return arma::stddev(inVec); }

/*** R
library(microbenchmark)
X <- rexp(500)
microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))

X <- rexp(5000)
microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))

X <- rexp(50000)
microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
*/

关于c++ - R stats::sd()与arma::stddev()与Rcpp实现的性能,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/24253228/

10-09 19:55
查看更多