【发布时间】：2014-06-16 22:33:32
【问题描述】：
为了进行 C++ / 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) );
}
我还决定测试 Armadillo 库中的 stddev 函数，考虑到它可以在向量上调用：
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
double armaSD(arma::colvec inVec)
{
return arma::stddev(inVec);
}
然后，我将这两个函数与基本 R 函数 sd() 进行基准测试，以获得一些不同大小的向量：
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
对于较小的样本，我的 C++ 函数 cppSD() 比 stats::sd() 快，但对于较大的向量则较慢，因为 stats::sd() 是矢量化的，我对此并不感到非常惊讶。但是，我没想到会从arma::stddev() 函数中看到这样的性能下降，因为它似乎也在以矢量化方式运行。我使用arma::stdev() 的方式是否存在问题，或者仅仅是stats::sd() 的编写方式（我假设为C）可以更有效地处理更大的向量？任何意见将不胜感激。
更新：
虽然我的问题最初是关于 arma::stddev 的正确使用，而不是试图找到最有效的方法来计算样本标准差，但看到 Rcpp::sd 糖函数执行非常有趣那么好。为了让事情变得更有趣，我将下面的 arma::stddev 和 Rcpp::sd 函数与我改编自 JJ Allaire 的 Rcpp Gallery 的两个帖子的 RcppParallel 版本进行了基准测试 - here 和 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
这是在我的笔记本电脑上运行 64 位 linux 和 i5-4200U CPU @ 1.60GHz 处理器；但我猜par_sd 和sugar_sd 之间的区别在Windows 机器上不会那么明显。
RcppParallel 版本的代码（相当长，并且需要 C++11 兼容的编译器来编译用于 InnerProduct 函子的重载 operator() 的 lambda 表达式）：
#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) );
}
【问题讨论】：
-
除了下面 Dirk Eddelbuettel 指出的问题之外，您的 cppSD() 函数还有其他问题。这不是计算标准偏差的可靠方法，因为它很容易受到 sum 变量中的浮点溢出的影响。 std::accumulate() 不是为数字设计的。相反，您的代码应该使用运行统计信息作为主要计算方法，或者作为后备。 cppSD() 函数也是低效的：它不必要地将数据复制到一个新的向量，然后重新写入向量。相反，它应该以只读方式使用原始数据
-
这值得重复：使用诸如 std::accumulate() 或 std::inner_product() 之类的函数来计算平均值和标准差是一个非常糟糕的主意：它们对浮点溢出不鲁棒或下溢。 Armadillo 的标准差函数要安全得多，因为它考虑了这些潜在问题。使用上述并行化方法，您只会更快地得到不好的结果。
标签： c++ r performance rcpp armadillo