非負の成分を得る最も速い方法

Question

ダブルベクトルの非負成分を得るより速い方法は何ですか？それは

pmax(x, 0) 

私の試みはRcpp使用され、次のとおりです。

//' @title Parallel maximum 
//' @description A faster \code{pmax()}. 
//' 
//' @name pmaxC 
//' @param x A numeric vector. 
//' @param a A single numeric value. 
//' @return The parallel maximum of the input values. 
//' @note This function will always be faster than \code{pmax(x, a)} when \code{a} is a single value, but can be slower than \code{pmax.int(x, a)} when \code{x} is short. Use this function when comparing a numeric vector with a single value. 
//' @export pmaxC 

#include <Rcpp.h> 
using namespace Rcpp; 

// [[Rcpp::export]] 
NumericVector pmaxC(NumericVector x, double a) { 
    int n = x.length(); 
    NumericVector out(n); 

    for (int i = 0; i < n; ++i) { 
    double xi = x[i]; 
    if (xi < a) { 
     out[i] = a; 
    } else { 
     out[i] = xi; 
    } 
    } 

    return out; 
} 

これはささやかな改善である：

set.seed(5) 
x <- rnorm(1e6) 

microbenchmark(pmax(x, 0), pmaxC(x, 0)) 
Unit: milliseconds 
     expr  min  lq  mean median  uq  max neval cld 
    pmax(x, 0) 8.500419 8.621341 11.09672 10.132045 10.791020 58.44972 100 a 
pmaxC(x, 0) 5.624480 5.709262 8.83968 7.598093 7.907853 53.91339 100 a 

どちらも許容できないほど遅いですが、一般的なシナリオそれは与えられた、Iパッケージがより速いアプローチを開発したかどうか疑問に思っていました。

Answer 1

実行している操作はかなりシンプルなので、上記のアルゴリズムに関して改善の余地があるとは思えません。しかし、の場合、実際には余分なパフォーマンスを引き出す必要がありますが、これは並列化に適した候補のようです。ここでRcppParallelを使用して可能な実装である：

// [[Rcpp::depends(RcppParallel)]] 
#include <RcppParallel.h> 
#include <Rcpp.h> 

struct Pmax : public RcppParallel::Worker { 
    struct Apply { 
     double mx; 
     Apply(double mx_) 
      : mx(mx_) 
     {} 

     double operator()(const double x) const 
     { 
      return x > mx ? x : mx; 
     } 
    }; 

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

    Apply f; 

    Pmax(const Rcpp::NumericVector input_, 
     Rcpp::NumericVector output_, 
     double mx_) 
    : input(input_), output(output_), f(mx_) 
    {} 

    void operator()(std::size_t begin, std::size_t end) 
    { 
     std::transform(
      input.begin() + begin, 
      input.begin() + end, 
      output.begin() + begin, 
      f 
     ); 
    } 
}; 

// [[Rcpp::export]] 
Rcpp::NumericVector par_pmax(Rcpp::NumericVector x, double y) 
{ 
    Rcpp::NumericVector res = Rcpp::no_init_vector(x.size()); 
    Pmax p(x, res, y); 
    RcppParallel::parallelFor(0, x.size(), p); 

    return res; 
} 

---

テストこれはあなたの例のデータを、私は合理的な改善を得る：

set.seed(5) 
x <- rnorm(1e6) 

all.equal(pmax(x, 0), par_pmax(x, 0)) 
#[1] TRUE 

microbenchmark::microbenchmark(
    pmax(x, 0), 
    pmaxC(x, 0), 
    par_pmax(x, 0), 
    times = 500L 
) 
# Unit: milliseconds 
#   expr  min  lq  mean median  uq  max neval 
#  pmax(x, 0) 11.843528 12.193126 14.972588 13.030448 16.799250 102.09895 500 
#  pmaxC(x, 0) 7.804883 8.036879 10.462070 8.772635 12.407587 69.08290 500 
# par_pmax(x, 0) 2.244691 2.443971 4.552169 2.624008 6.359027 65.99233 500