Производительность R stats:: sd() vs. arma:: stddev() по сравнению с реализацией Rcpp

Просто для того, чтобы работать над моим программированием на С++/Rcpp, я сделал попытку реализовать стандартную функцию отклонения (выборки):

Я также решил проверить функцию stddev из библиотеки Armadillo, учитывая, что ее можно вызвать по вектору:

Затем я сравнил эти две функции с базовой функцией R sd() для нескольких векторов разного размера:

Меня не удивило, что моя С++-функция cppSD() была быстрее, чем stats::sd() для более мелких образцов, но медленнее для векторов большего размера, так как stats::sd() векторизован. Тем не менее, я не ожидал такого снижения производительности от функции arma::stddev(), поскольку он, похоже, также работает в векторном режиме. Есть ли проблема с тем, как я использую arma::stdev(), или это просто, что stats::sd() был написан (в C я предполагаю) таким образом, что он может обрабатывать более крупные векторы гораздо эффективнее? Любой вход будет оценен.

Хотя мой вопрос был первоначально о правильном использовании arma::stddev, а не о попытке найти наиболее эффективный способ расчета стандартного отклонения выборки, очень интересно видеть, что функция сахара Rcpp::sd выполняет так хорошо. Чтобы сделать вещи более интересными, я сравнивал функции arma::stddev и Rcpp::sd ниже с версией RcppParallel, которую я адаптировал из двух сообщений JJ Allaire Rcpp Gallery - здесь и здесь:

Это было на моем ноутбуке с 64-битным Linux с процессором i5-4200U с процессором с частотой 1,60 ГГц; но я предполагаю, что разница между par_sd и sugar_sd будет менее существенной на машине Windows.

И код для версии RcppParallel (который значительно длиннее и требует компилятора, совместимого с С++ 11 для выражения лямбда, используемого в перегруженном operator() функтора InnerProduct):

Ответы

Ответ 1

Вы сделали тонкую ошибку в том, как вы создаете объект Armadillo, что приводит к копиям и, следовательно, к ухудшению производительности.

Вместо этого используйте интерфейс const arma::colvec & 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

Я думаю, что функция sd, построенная в салоне Rcpp, намного эффективнее. См. Код ниже:

   #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