三次函数计算的速度优化

Question

我需要(sin(x)/x)^3在R中多次评估函数.最快的方法是什么:

1. (sin(x)/x)^3
2. (sin(x)/x)^3L
3. { y=sin(x)/x; y*y*y } ？

Answer 1

tl; dr第三个选项(y*y*y)是最快的; 转换为C++有一点帮助,但不是我们所期望的数量级(只有大约20%的改进),因为R在向量化时已经非常有效地完成了这项任务.

使用microbenchmark包找出...

library(microbenchmark)
x <- 10
m <- microbenchmark((sin(x)/x)^3,(sin(x)/x)^3L,
                     {y=sin(x)/x; y*y*y}, times=1e4)

## Unit: nanoseconds
##                  expr min   lq     mean median   uq    max neval
##          (sin(x)/x)^3   1 1524 1795.508   1576 1654 220730 10000
##         (sin(x)/x)^3L   1 1503 1766.368   1558 1633 216711 10000
## { y=sin(x)/x; y*y*y }   2 1623 1925.608   1692 1785 243385 10000

现在尝试矢量化版本(对于长度为10 ^ 5的向量),包括Rcpp-ized版本:

set.seed(101)
x <- rnorm(1e5)
library(Rcpp)
sourceCpp("cubebench.cpp")
m2 <- microbenchmark((sin(x)/x)^3,(sin(x)/x)^3L,
                 {y=sin(x)/x; y*y*y},
                 sin_cube(x),
                 sin_cubepow(x), times=100)

## Unit: milliseconds
##                    expr   min     lq   mean median     uq     max
##            (sin(x)/x)^3 9.512 10.284 10.685 10.492 10.785  13.212
##           (sin(x)/x)^3L 9.956 10.480 11.902 10.735 11.125 105.164
##  { y=sin(x)/x;  y*y*y } 2.455  2.855  3.348  3.063  3.541   5.356
##             sin_cube(x) 1.906  2.278  2.611  2.355  2.785   4.732
##          sin_cubepow(x) 8.331  9.180  9.804  9.515  9.960  13.931

令人惊讶的是,对于更长的向量,第三种选择更快.C++变体与相应的R版本没有什么不同.

照片:

comb <- rbind(data.frame(w="short",as.data.frame(m)),
              data.frame(w="long",as.data.frame(m2)))
library(ggplot2); theme_set(theme_bw())
library(ggstance)
ggplot(comb,aes(time,expr))+geom_violinh(fill="gray")+
    scale_x_log10()+
    labs(x="time (ns)",y="")+
    facet_grid(.~w,scale="free")

• 通过迄今为止最重要的性能提示是向量化的计算尽可能(您节省至少10这样的一个因素)
• 如果每次计算1000到10000纳秒之间的差异对您很重要,您可能需要使用"更接近金属"的计算平台(C/C++/Rcpp或Julia或......)

这是cubebench.cpp:

// Hacked from http://gallery.rcpp.org/articles/run_sum-benchmark/
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector sin_cube(NumericVector x) {

    int sz = x.size();
    double y;
    NumericVector res(sz);

    // loop through the vector
    for(int i = 0; i < sz; i++){
        y = sin(x[i])/x[i];
        res[i] = y*y*y;
    }
    return res;
}

// [[Rcpp::export]]
NumericVector sin_cubepow(NumericVector x) {

    int sz = x.size();
    double y;
    NumericVector res(sz);

    // loop through the vector
    for(int i = 0; i < sz; i++){
        y = sin(x[i])/x[i];
        res[i] = pow(y,3.0);
    }
    return res;
}