在 Rcpp 中对矩阵的所有列或所有行进行快速选择的多线程最快方法 - OpenMP、RcppParallel 或 RcppThread
Fastest way to multithread doing quickselect on all columns or all rows of a matrix in Rcpp - OpenMP, RcppParallel or RcppThread
我正在使用此 Rcpp 代码对值向量执行 quickselect，即在 O(n) 时间内从向量中获取第 k 个最大元素（我将其保存为 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 矩阵的所有列或所有行执行多线程快速选择，并且 return 结果作为向量。因为我是 Rcpp 的新手，所以我想得到一些建议，尽管哪个框架可能最快并且最容易编码（它必须轻松跨平台工作并且我需要很好地控制 nr要使用的线程）。使用 OpenMP, RcppParallel or RcppThread？或者甚至更好 - 如果有人可以展示一种快速而优雅的方法来做到这一点？
是的，这将是多线程变体的候选者——但正如 RcppParallel 文档将告诉您的那样，并行代码的一个要求是非 R 内存，这里我们在高效 中使用 RcppArmadillo零拷贝 方式意味着它是 R 内存。
因此您可能需要权衡并行执行的额外数据副本（例如，RMatrix 类型 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
令我惊讶的是，在不使用多线程（qSelectMbycol 函数）的情况下，在 Rcpp 中执行应用程序的速度提高了约 2 倍，并且使用 OpenMP 多线程（qSelectMbycolOpenMP）速度进一步提高了 2 倍。
欢迎就可能的代码优化提出任何建议...
对于小的 n (n<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
我正在使用此 Rcpp 代码对值向量执行 quickselect，即在 O(n) 时间内从向量中获取第 k 个最大元素（我将其保存为 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 矩阵的所有列或所有行执行多线程快速选择，并且 return 结果作为向量。因为我是 Rcpp 的新手，所以我想得到一些建议，尽管哪个框架可能最快并且最容易编码（它必须轻松跨平台工作并且我需要很好地控制 nr要使用的线程）。使用 OpenMP, RcppParallel or RcppThread？或者甚至更好 - 如果有人可以展示一种快速而优雅的方法来做到这一点？
是的，这将是多线程变体的候选者——但正如 RcppParallel 文档将告诉您的那样，并行代码的一个要求是非 R 内存，这里我们在高效 中使用 RcppArmadillo零拷贝 方式意味着它是 R 内存。
因此您可能需要权衡并行执行的额外数据副本（例如，RMatrix 类型 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
令我惊讶的是，在不使用多线程（qSelectMbycol 函数）的情况下，在 Rcpp 中执行应用程序的速度提高了约 2 倍，并且使用 OpenMP 多线程（qSelectMbycolOpenMP）速度进一步提高了 2 倍。
欢迎就可能的代码优化提出任何建议...
对于小的 n (n<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