R-最快的方法在RCPP -OpenMP，rcppparallal或rcppthread中对所有列或所有矩阵的所有列进行

我正在使用此RCPP代码对值的向量进行快速选择，即从O（n）时间中的矢量获取kth最大元素（我将其保存为qselect.cpp）：

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;
// [[Rcpp::export]]
double qSelect(arma::vec& x, const int k) {
  // ARGUMENTS
  // x: vector to find k-th largest element in
  // k: desired k-th largest element
  // safety copy since nth_element modifies in place
  arma::vec y(x.memptr(), x.n_elem);
  // partially sort y in O(n) time
  std::nth_element(y.begin(), y.begin() + k - 1, y.end());
  // the k-th largest value
  const double kthValue = y(k-1);
  return kthValue;
}

我正在使用它作为计算所需百分位数的快速方法。例如

n = 50000
set.seed(1)
x = rnorm(n=n, mean=100, sd=20)
tau = 0.01 # desired percentile
k = tau*n+1 # here we will get the 6th largest element
library(Rcpp)
Rcpp::sourceCpp('qselect.cpp')
library(microbenchmark)
microbenchmark(qSelect(x,k)) # 53.32917, 548 µs
microbenchmark(sort(x, partial=k)[k]) # 53.32917, 694 µs = pure R solution

[这看起来可能已经很快，但是我需要在我的应用程序中进行数百万次的时间]

现在，我想修改此RCPP函数，以便在所有列或R矩阵的所有列或所有行上进行多线程QuickSelect，并将结果作为向量返回。由于我在RCPP中有点新手，所以我想提供一些建议，尽管在哪个框架上可能是最快的。最容易代码（必须轻松地工作跨平台，并且我需要对线程的NR进行良好的控制）。使用OpenMP，rcppparallel或rcppthread？甚至更好 - 如果某人也许可以表现出一种快速而优雅的方法？

是的有效零复制的方式意味着它是R内存。

因此，您可能需要交易额外的数据副本（例如， RMatrix type rcppparallel使用），并行执行。

但是，由于您的算法很简单且列，因此您还可以在上面的功能中使用一个OpenMP循环进行试验：将其传递给矩阵，请使用#pragma for将其循环在列上。

按照下面的建议，我尝试使用OpenMP进行多线程，这似乎可以使用笔记本电脑上的8个线程给出不错的加速。我将qselect.cpp文件修改为：

// [[Rcpp::depends(RcppArmadillo)]]
#define RCPP_ARMADILLO_RETURN_COLVEC_AS_VECTOR
#include <RcppArmadillo.h>
using namespace arma;
// [[Rcpp::export]]
double qSelect(arma::vec& x, const int k) {
  // ARGUMENTS
  // x: vector to find k-th largest element in
  // k: k-th statistic to look up
  // safety copy since nth_element modifies in place
  arma::vec y(x.memptr(), x.n_elem);
  // partially sorts y
  std::nth_element(y.begin(), y.begin() + k - 1, y.end());
  // the k-th largest value
  const double kthValue = y(k-1);
  return kthValue;
}

// [[Rcpp::export]]
arma::vec qSelectMbycol(arma::mat& M, const int k) {
  // ARGUMENTS
  // M: matrix for which we want to find the k-th largest elements of each column
  // k: k-th statistic to look up
  arma::mat Y(M.memptr(), M.n_rows, M.n_cols);
  // we apply over columns
  int c = M.n_cols;
  arma::vec out(c);
  int i;
  for (i = 0; i < c; i++) {
      arma::vec y = Y.col(i);
      std::nth_element(y.begin(), y.begin() + k - 1, y.end());
      out[i] = y(k-1); // the k-th largest value of each column
  }
  return out;
}
#include <omp.h>
// [[Rcpp::plugins(openmp)]]
// [[Rcpp::export]]
arma::vec qSelectMbycolOpenMP(arma::mat& M, const int k, int nthreads) {
  // ARGUMENTS
  // M: matrix for which we want to find the k-th largest elements of each column
  // k: k-th statistic to look up
  // nthreads: nr of threads to use
  arma::mat Y(M.memptr(), M.n_rows, M.n_cols);
  // we apply over columns
  int c = M.n_cols;
  arma::vec out(c);
  int i;
  omp_set_num_threads(nthreads);
#pragma omp parallel for shared(out) schedule(dynamic,1)
  for (i = 0; i < c; i++) {
    arma::vec y = Y.col(i);
    std::nth_element(y.begin(), y.begin() + k - 1, y.end());
    out(i) = y(k-1); // the k-th largest value of each column
  }
  return out;
}

基准：

n = 50000
set.seed(1)
x = rnorm(n=n, mean=100, sd=20)
M = matrix(rnorm(n=n*10, mean=100, sd=20), ncol=10)
tau = 0.01 # desired percentile
k = tau*n+1 # we will get the 6th smallest element
library(Rcpp)
Rcpp::sourceCpp('qselect.cpp')
library(microbenchmark
microbenchmark(apply(M, 2, function (col) sort(col, partial=k)[k]),
               apply(M, 2, function (col) qSelect(col,k)),
               qSelectMbycol(M,k),
               qSelectMbycolOpenMP(M,k,nthreads=8))[,1:4]
Unit: milliseconds
                                                 expr      min       lq      mean    median        uq        max neval cld
 apply(M, 2, function(col) sort(col, partial = k)[k]) 8.937091 9.301237 11.802960 11.828665 12.718612  43.316107   100   b
           apply(M, 2, function(col) qSelect(col, k)) 6.757771 6.970743 11.047100  7.956696  9.994035 133.944735   100   b
                                  qSelectMbycol(M, k) 5.370893 5.526772  5.753861  5.641812  5.826985   7.124698   100  a 
              qSelectMbycolOpenMP(M, k, nthreads = 8) 2.695924 2.810108  3.005665  2.899701  3.061996   6.796260   100  a

我对在RCPP中进行申请的速度的Ca 2倍增益感到惊讶

关于可能的代码优化的任何建议，但欢迎...

对于小型n（n＆lt; 1000），OpenMP版本的速度不快，也许是因为个人工作太小了。例如。对于n=500：

Unit: microseconds
                                                 expr     min       lq      mean   median       uq      max neval cld
 apply(M, 2, function(col) sort(col, partial = k)[k]) 310.477 324.8025 357.47145 337.8465 361.5810 1782.885   100   c
           apply(M, 2, function(col) qSelect(col, k)) 103.921 114.8255 141.59221 119.3155 131.9315 1990.298   100  b 
                                  qSelectMbycol(M, k)  24.377  32.2885  44.13873  35.2825  39.3440  900.210   100 a  
              qSelectMbycolOpenMP(M, k, nthreads = 8)  76.123  92.1600 130.42627  99.8575 112.4730 1303.059   100  b

相关内容

最新更新

热门标签：