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

Posted

技术标签:

【中文标题】R stats::sd() 与 arma::stddev() 与 Rcpp 实现的性能【英文标题】:Performance of R stats::sd() vs. arma::stddev() vs. Rcpp implementation 【发布时间】:2014-06-16 22:33:32 【问题描述】:

为了进行 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::stddevRcpp::sd 函数与我改编自 JJ Allaire 的 Rcpp Gallery 的两个帖子的 RcppParallel 版本进行了基准测试 - here 和 here:

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 兼容的编译器来编译用于 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) );

【问题讨论】:

除了下面 Dirk Eddelbuettel 指出的问题之外,您的 cppSD() 函数还有其他问题。这不是计算标准偏差的可靠方法,因为它很容易受到 sum 变量中的浮点溢出的影响。 std::accumulate() 不是为数字设计的。相反,您的代码应该使用运行统计信息作为主要计算方法,或者作为后备。 cppSD() 函数也是低效的:它不必要地将数据复制到一个新的向量,然后重新写入向量。相反,它应该以只读方式使用原始数据 这值得重复:使用诸如 std::accumulate() 或 std::inner_product() 之类的函数来计算平均值和标准差是一个非常糟糕的主意:它们对浮点溢出不鲁棒或下溢。 Armadillo 的标准差函数要安全得多,因为它考虑了这些潜在问题。使用上述并行化方法,您只会更快地得到不好的结果。 【参考方案1】:

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

改用const arma::colvec &amp; invec的接口,一切都好:

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))
*/

【讨论】:

【参考方案2】:

我认为 Rcpp 糖中内置的sd 函数效率更高。请看下面的代码:

   #include <RcppArmadillo.h>
   //[[Rcpp::depends(RcppArmadillo)]]
   #include <vector>
   #include <cmath>
   #include <numeric>
   using namespace Rcpp;
   //[[Rcpp::export]]                                                                                                                                                         
  double sd_cpp(NumericVector& xin)
 std::vector<double> xres(xin.begin(),xin.end());
 int n=xres.size();
 double sum=std::accumulate(xres.begin(),xres.end(),0.0);
 double mean=sum/n;
 for(std::vector<double>::iterator iter=xres.begin();iter!=xres.end();++iter)
   double tmp=(*iter-mean)*(*iter-mean);
   *iter=tmp;
 
   double sd=std::accumulate(xres.begin(),xres.end(),0.0);
   return std::sqrt(sd/(n-1));
 
  //[[Rcpp::export]]
 double sd_arma(arma::colvec& xin)
 return arma::stddev(xin);

 //[[Rcpp::export]]
 double sd_sugar(NumericVector& xin)
 return sd(xin);


> sourcecpp("sd.cpp")

> microbenchmark(sd(X),sd_cpp(X),sd_arma(X),sd_sugar(X))
   Unit: microseconds
      expr    min      lq     mean  median      uq     max neval
      sd(X) 47.655 49.4120 51.88204 50.5395 51.1950 113.643   100
  sd_cpp(X) 28.145 28.4410 29.01541 28.6695 29.4570  37.118   100
  sd_arma(X) 23.706 23.9615 24.65931 24.1955 24.9520  50.375   100
 sd_sugar(X) 19.197 19.478 20.38872 20.0785 21.2015  28.664   100

【讨论】:

使用 std::accumulate() 计算平均值和标准差是一个非常糟糕的主意:它对浮点上溢或下溢不鲁棒。 Armadillo 的标准差函数要安全得多,因为它考虑了这些潜在问题。【参考方案3】:

关于Rcpp函数的2个问题:

1) 您想要没有均值的标准差的可能性有多大?如果要求没有均值的 SD 并不常见,为什么不返回两者,而不是要求 R 用户进行 2 次函数调用,这实际上计算了两次均值。

2) 有什么理由在函数内部克隆向量吗?

using namespace Rcpp;
// [[Rcpp::plugins(cpp14)]]

// [[Rcpp::export]]
NumericVector cppSD(NumericVector rinVec)

  double n = (double)rinVec.size();
  double sum = 0;
  for (double& v : rinVec)
     sum += v;
  double mean = sum / n;
  double varianceNumerator = 0;
  for(double& v : rinVec)
     varianceNumerator += (v - mean) * (v - mean);
  return NumericVector::create(_["std.dev"] = sqrt(varianceNumerator/ (n - 1.0)),
                               _["mean"] = mean);

【讨论】:

以上是关于R stats::sd() 与 arma::stddev() 与 Rcpp 实现的性能的主要内容,如果未能解决你的问题,请参考以下文章

R语言入门 - R的安装与帮助

java 选择排序与冒泡排序

一R语言简介R与统计

origin2021与r语言哪个版本匹配

r 与R的KMeans

r GA与R