R stats::sd()与arma::stddev()与Rcpp实现的性能

只是为了在我的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) );
}

我还决定测试Armadillo库中的stddev函数，考虑到它可以在矢量上调用:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
double armaSD(arma::colvec inVec)
{
  return arma::stddev(inVec);
}

然后我将这两个函数与基本R函数sd()进行基准测试，用于几个大小不同的向量:

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()是矢量化的，我对此并不感到非常惊讶。然而，我没有预料到arma::stddev()函数会出现这样的性能下降，因为它似乎也以矢量化的方式运行。我使用arma::stdev()的方式是否有问题，或者仅仅是stats::sd()被写成(我假设是在C中)这样一种方式，它可以更有效地处理更大的向量?如有任何意见，将不胜感激。

:

虽然我的问题最初是关于arma::stddev的正确使用，而不是试图找到最有效的方法来计算样本标准偏差，但看到Rcpp::sd糖函数表现如此之好是非常有趣的。为了让事情变得更有趣，我将下面的arma::stddev和Rcpp::sd函数与我从JJ Allaire的Rcpp画廊的两篇文章中改编的RcppParallel版本进行了基准测试:

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_sd和sugar_sd之间的差异在Windows机器上不会太大。

RcppParallel版本的代码(它相当长，并且需要一个c++ 11兼容的编译器来处理InnerProduct函子的重载operator()中使用的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) );
}