2014-06-16 17 views
15

Solo per lo scopo di lavorare sulla mia programmazione C++/Rcpp, ho preso una decisione sull'implementazione di una funzione di deviazione standard (campione) :Prestazioni di statistiche R :: sd() vs. arma :: stddev() rispetto all'implementazione di 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)); 
} 

ho anche deciso di testare la funzione stddev dalla libreria Armadillo, considerando che può essere chiamato un vettore:

#include <RcppArmadillo.h> 

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

using namespace Rcpp; 

// [[Rcpp::export]] 

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

Poi Benchmarked queste due funzioni contro la funzione di base R sd() per alcuni vettori di dimensioni variabili:

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 

non ero terribilmente sorpreso per il fatto che il mio C++ funzione cppSD() era più veloce di stats::sd() per i campioni più piccoli, ma più lento per i vettori di dimensioni più grandi in quanto è stats::sd() Vettorializzare. Tuttavia, non mi aspettavo di vedere un tale degrado delle prestazioni dalla funzione arma::stddev() poiché sembra funzionare anche in modo vettorializzato. C'è un problema con il modo in cui sto usando arma::stdev(), o è semplicemente che stats::sd() è stato scritto (in C sto assumendo) in modo tale da poter gestire i vettori più grandi in modo molto più efficiente? Qualsiasi input sarebbe apprezzato.


Aggiornamento:

Anche se la mia domanda era in origine circa l'uso corretto dei arma::stddev, e non tanto di cercare trovare il modo più efficiente possibile calcolare la deviazione standard del campione, è molto interessante per vedere che la funzione zucchero Rcpp::sd si comporta così bene. Per rendere le cose un po 'più interessante, ho benchmark le arma::stddev e Rcpp::sd funzioni basso contro una versione RcppParallel che ho adattato da due dei messaggi Rcpp Galleria di JJ Allaire - here e 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 

Questo è stato il mio portatile in esecuzione 64 -bit linux con processore i5-4200U a 1.60GHz; ma suppongo che la differenza tra par_sd e sugar_sd sia meno sostanziale su una macchina Windows.

E il codice per la versione RcppParallel (che è notevolmente più lungo, e richiede un compilatore compatibile C++ 11 per l'espressione lambda utilizzato in sovraccarico operator() del InnerProduct funtore):

#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)); 
} 
+0

In aggiunta al problema indicato da Dirk Eddelbuettel di seguito, la funzione cppSD() ha ulteriori problemi. Non è un metodo robusto per calcolare la deviazione standard, poiché è facilmente influenzato da overflow in virgola mobile nella variabile sum. std :: accumulate() non è progettato per la numerics. Invece, il tuo codice dovrebbe utilizzare statistiche in esecuzione come metodo di calcolo primario o come fallback. La funzione cppSD() è anche inefficiente: copia inutilmente i dati su un vettore nuovo e quindi riscrive il vettore. Invece, dovrebbe usare i dati originali in modo di sola lettura – mtall

+1

Questo vale ripetere: l'uso di funzioni come std :: accumulate() o std :: inner_product() per calcolare la media e la deviazione standard è una pessima idea: non sono robusti per overflow o underflow in virgola mobile. La funzione di deviazione standard di Armadillo è molto più sicura, in quanto tiene conto di questi potenziali problemi. Usando il metodo di parallelizzazione sopra riportato otterrai risultati cattivi più rapidamente. – mtall

risposta

13

hai fatto un piccolo errore nel modo in cui istanziate l'oggetto Armadillo - che porta a copie e quindi a prestazioni degradate.

utilizzare un'interfaccia di const arma::colvec & invec invece, e tutto è buono:

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> 

che si basa sulla mia versione di codice in cui tutto è un unico file e armaSD2 è definito come ho suggerito - che porta alla vincita prestazione.

#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

penso che la funzione sd costruito in zucchero Rcpp è molto più efficiente.Vedere il codice di seguito:

#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 
+1

L'utilizzo di std :: accumulate() per calcolare la media e la deviazione standard è una pessima idea: non è robusto per gli overflow o i underflow in virgola mobile. La funzione di deviazione standard di Armadillo è molto più sicura, in quanto tiene conto di questi potenziali problemi. – mtall