计算Hawk过程梯度的有效方法
Efficient way to calculate Hawk's process gradient
我有兴趣计算以下数量
B(i) = \sum_{j < i}(x_i-x_j)exp^{-\beta(x_i - x_j)}
这是计算 Hawk 过程可能性参数之一的梯度的一部分（更多信息可在此处找到：http://www.ism.ac.jp/editsec/aism/pdf/031_1_0145.pdf）。
Beta 只是一个常数，用于抖动问题 x_i 是我的第 i 个数据点。
我正在尝试使用以下代码块在 RCPP 中计算上述数量：
for( int i = 1; i< x.size();i++) {
double temp=0;
for(int j=0; j<=i-1;j++){
temp+=(x[i]-x[j])*exp(-beta*(x[i]-x[j]));
}
但它非常低效且缓慢。关于如何加速这个公式的任何建议？
这是一个 O(N^2) 操作，没有考虑 exp 的成本。任何调整都可能产生最小的改进。
一些简单的建议：
- 缓存 外循环中
x[i] 的值，因为您在内循环中重复对其进行子集化。
- 从使用
exp(-beta * ..) 切换到 1/exp(beta*(x ... ))
- 使用
++i 而不是 i++ 来避免 a slight performance hiccup 因为你避免了后者所做的 i 的副本。
原代码：
#include<Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_org(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
for (int i = 1; i < n; i++) {
double temp = 0;
for (int j = 0; j <= i - 1; j++) {
temp += (x[i] - x[j]) * exp(-beta * (x[i] - x[j]));
}
B(i - 1) = temp;
}
return B;
}
修改后的代码：
#include<Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
double x_i;
for (int i = 1; i < n; ++i) {
double temp = 0;
x_i = x[i];
for (int j = 0; j <= i - 1; ++j) {
temp += (x_i - x[j]) * 1 / exp(beta * (x_i - x[j]));
}
B(i - 1) = temp;
}
return B;
}
测试
set.seed(111)
x = rnorm(1e4)
all.equal(
hawk_process_org(x),
hawk_process_cache(x)
)
#> [1] TRUE
bench_func = microbenchmark::microbenchmark(
hawk_process_org(x),
hawk_process_cache(x)
)
bench_func
#> Unit:milliseconds
#> expr min lq mean median uq max neval
#> hawk_process_org(x) 436.5349 465.9674 505.9606 481.4703 500.6652 894.7477 100
#> hawk_process_cache(x) 446.0499 454.9098 485.3830 468.6580 494.9457 799.0940 100
因此，根据建议，您会获得略微更好的结果。
标准操作在 C++ 中非常快（+、- 等）。
然而，exp 计算起来更复杂，所以更慢。
因此，如果我们想要提高一些性能，则更有可能预先计算 exp 计算。
这里，B(i) = \sum_{j < i}(x_i-x_j)exp^{-\beta(x_i - x_j)} 等同于 B(i) = \sum_{j < i}(x_i-x_j) / exp^{\beta x_i} * exp^{\beta x_j} 这样你就可以预先计算每个索引的 exp （并且也把依赖于 i 的那个放在循环）。通过重构它，您可以进行其他预计算。所以，我把之前的两个解决方案放在这里，然后是我的增量解决方案：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_org(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
for (int i = 1; i < n; i++) {
double temp = 0;
for (int j = 0; j <= i - 1; j++) {
temp += (x[i] - x[j]) * exp(-beta * (x[i] - x[j]));
}
B(i - 1) = temp;
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
double x_i;
for (int i = 1; i < n; ++i) {
double temp = 0;
x_i = x[i];
for (int j = 0; j <= i - 1; ++j) {
temp += (x_i - x[j]) * 1 / exp(beta * (x_i - x[j]));
}
B(i - 1) = temp;
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_2(Rcpp::NumericVector x,
double beta = 3) {
int i, j, n = x.size();
Rcpp::NumericVector B(n);
Rcpp::NumericVector x_exp = exp(beta * x);
double temp;
for (i = 1; i < n; i++) {
temp = 0;
for (j = 0; j < i; j++) {
temp += (x[i] - x[j]) * x_exp[j] / x_exp[i];
}
B[i] = temp;
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_3(Rcpp::NumericVector x,
double beta = 3) {
int i, j, n = x.size();
Rcpp::NumericVector B(n);
Rcpp::NumericVector x_exp = exp(beta * x);
double temp;
for (i = 1; i < n; i++) {
temp = 0;
for (j = 0; j < i; j++) {
temp += (x[i] - x[j]) * x_exp[j];
}
B[i] = temp / x_exp[i];
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_4(Rcpp::NumericVector x,
double beta = 3) {
Rcpp::NumericVector exp_pre = exp(beta * x);
Rcpp::NumericVector exp_pre_cumsum = cumsum(exp_pre);
Rcpp::NumericVector x_exp_pre_cumsum = cumsum(x * exp_pre);
return (x * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_5(Rcpp::NumericVector x,
double beta = 3) {
int n = x.size();
NumericVector B(n);
double exp_pre, exp_pre_cumsum = 0, x_exp_pre_cumsum = 0;
for (int i = 0; i < n; i++) {
exp_pre = exp(beta * x[i]);
exp_pre_cumsum += exp_pre;
x_exp_pre_cumsum += x[i] * exp_pre;
B[i] = (x[i] * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
return B;
}
/*** R
set.seed(111)
x = rnorm(1e3)
all.equal(
hawk_process_org(x),
hawk_process_cache(x)
)
all.equal(
hawk_process_org(x),
hawk_process_cache_2(x)[-1]
)
all.equal(
hawk_process_org(x),
hawk_process_cache_3(x)[-1]
)
all.equal(
hawk_process_org(x),
hawk_process_cache_4(x)[-1]
)
all.equal(
hawk_process_org(x),
hawk_process_cache_5(x)[-1]
)
microbenchmark::microbenchmark(
hawk_process_org(x),
hawk_process_cache(x),
hawk_process_cache_2(x),
hawk_process_cache_3(x),
hawk_process_cache_4(x),
hawk_process_cache_5(x)
)
*/
x = rnorm(1e3) 的基准：
Unit: microseconds
expr min lq mean median uq max neval cld
hawk_process_org(x) 19801.686 20610.0365 21017.89339 20816.1385 21157.4900 25548.042 100 d
hawk_process_cache(x) 20506.903 21062.1370 21534.47944 21297.8710 21775.2995 26030.106 100 e
hawk_process_cache_2(x) 1895.809 2038.0105 2087.20696 2065.8220 2103.0695 3212.874 100 c
hawk_process_cache_3(x) 430.084 458.3915 494.09627 474.2840 503.0885 1580.282 100 b
hawk_process_cache_4(x) 50.657 55.2930 71.60536 57.6105 63.5700 1190.260 100 a
hawk_process_cache_5(x) 43.373 47.0155 60.43775 49.6640 55.6235 842.288 100 a
这比试图从小的优化中获得纳秒要有效得多，小的优化可能会使您的代码更难阅读。
不过，让我们试试@coatless 在我最后的解决方案中提出的优化：
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_6(Rcpp::NumericVector x,
double beta = 3) {
int n = x.size();
NumericVector B = Rcpp::no_init(n);
double x_i, exp_pre, exp_pre_cumsum = 0, x_exp_pre_cumsum = 0;
for (int i = 0; i < n; ++i) {
x_i = x[i];
exp_pre = exp(beta * x_i);
exp_pre_cumsum += exp_pre;
x_exp_pre_cumsum += x_i * exp_pre;
B[i] = (x_i * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
return B;
}
x = rnorm(1e6) 的基准：
Unit: milliseconds
expr min lq mean median uq max neval cld
hawk_process_cache_5(x) 42.52886 43.53653 45.28427 44.46688 46.74129 57.38046 100 a
hawk_process_cache_6(x) 42.14778 43.19054 45.93252 44.28445 46.51052 153.30447 100 a
还是不太有说服力..
有趣的问题。在我的测试中，将这两个答案结合起来确实可以进一步提高性能（进一步降低基准）：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector hawk_process_cache_combined(NumericVector x,
double beta = 3) {
int n = x.size();
NumericVector B = Rcpp::no_init(n-1);
double exp_pre(exp(beta * x[0]));
double exp_pre_cumsum(exp_pre);
double x_exp_pre_cumsum(x[0] * exp_pre);
double x_i;
for (int i = 1; i < n; ++i) {
x_i = x[i];
exp_pre = exp(beta * x_i);
exp_pre_cumsum += exp_pre;
x_exp_pre_cumsum += x_i * exp_pre;
B[i-1] = (x_i * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
return B;
}
all.equal(
hawk_process_org(x),
hawk_process_cache_combined(x)
)
#> [1] TRUE
虽然原来的公式是 "embarrassingly parallel"，但此表达式不再是这种情况。但是，像 cumsum 这样的前缀扫描算法也可以并行化。像 ArrayFire provide interfaces to such algorithms using the GPU. Using RcppArrayFire 这样的库可以基于 F.Privé 的 hawk_process_cached_4:
编写
// [[Rcpp::depends(RcppArrayFire)]]
#include <RcppArrayFire.h>
// [[Rcpp::export]]
af::array hawk_process_af(RcppArrayFire::typed_array<f32> x,
double beta = 3) {
af::array exp_pre = exp(beta * x);
af::array exp_pre_cumsum = af::accum(exp_pre);
af::array x_exp_pre_cumsum = af::accum(x * exp_pre);
af::array result = (x * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
return result(af::seq(1, af::end));
}
这里的结果并不完全相等，因为我的driver/card只支持单精度浮点数：
all.equal(
hawk_process_org(x),
hawk_process_af(x)
)
#> [1] "Mean relative difference: 3.437819e-07"
如果使用双精度，可以在上面写 f64 并获得相同的结果。现在进行基准测试：
set.seed(42)
x <- rnorm(1e3)
microbenchmark::microbenchmark(
hawk_process_af(x),
hawk_process_cache_combined(x),
hawk_process_cache_5(x)[-1]
)
#> Unit: microseconds
#> expr min lq mean median uq max neval
#> hawk_process_af(x) 245.281 277.4625 338.92232 298.5410 346.576 1030.045 100
#> hawk_process_cache_combined(x) 35.343 39.0120 43.69496 40.7770 45.264 84.242 100
#> hawk_process_cache_5(x)[-1] 52.408 57.8580 65.55799 60.5265 67.965 125.864 100
x <- rnorm(1e6)
microbenchmark::microbenchmark(
hawk_process_af(x),
hawk_process_cache_combined(x),
hawk_process_cache_5(x)[-1]
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> hawk_process_af(x) 27.54936 28.42794 30.93452 29.20025 32.40667 49.41888 100
#> hawk_process_cache_combined(x) 34.00380 36.84497 40.74862 39.03649 41.85902 111.51628 100
#> hawk_process_cache_5(x)[-1] 47.02501 53.24702 57.94747 55.35018 58.42097 130.89737 100
所以对于小向量，组合方法更快，而一旦卸载到 GPU 得到回报，时间更长。所有这一切都不是使用一些高功率 GPU，而是简单的板载图形：
RcppArrayFire::arrayfire_info()
#> ArrayFire v3.5.1 (OpenCL, 64-bit Linux, build 0a675e8)
#> [0] BEIGNET: Intel(R) HD Graphics Skylake ULT GT2, 4096 MB
我有兴趣计算以下数量
B(i) = \sum_{j < i}(x_i-x_j)exp^{-\beta(x_i - x_j)}
这是计算 Hawk 过程可能性参数之一的梯度的一部分（更多信息可在此处找到：http://www.ism.ac.jp/editsec/aism/pdf/031_1_0145.pdf）。
Beta 只是一个常数，用于抖动问题 x_i 是我的第 i 个数据点。
我正在尝试使用以下代码块在 RCPP 中计算上述数量：
for( int i = 1; i< x.size();i++) {
double temp=0;
for(int j=0; j<=i-1;j++){
temp+=(x[i]-x[j])*exp(-beta*(x[i]-x[j]));
}
但它非常低效且缓慢。关于如何加速这个公式的任何建议？
这是一个 O(N^2) 操作，没有考虑 exp 的成本。任何调整都可能产生最小的改进。
一些简单的建议：
- 缓存 外循环中
x[i]的值，因为您在内循环中重复对其进行子集化。 - 从使用
exp(-beta * ..)切换到1/exp(beta*(x ... )) - 使用
++i而不是i++来避免 a slight performance hiccup 因为你避免了后者所做的i的副本。
原代码：
#include<Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_org(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
for (int i = 1; i < n; i++) {
double temp = 0;
for (int j = 0; j <= i - 1; j++) {
temp += (x[i] - x[j]) * exp(-beta * (x[i] - x[j]));
}
B(i - 1) = temp;
}
return B;
}
修改后的代码：
#include<Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
double x_i;
for (int i = 1; i < n; ++i) {
double temp = 0;
x_i = x[i];
for (int j = 0; j <= i - 1; ++j) {
temp += (x_i - x[j]) * 1 / exp(beta * (x_i - x[j]));
}
B(i - 1) = temp;
}
return B;
}
测试
set.seed(111)
x = rnorm(1e4)
all.equal(
hawk_process_org(x),
hawk_process_cache(x)
)
#> [1] TRUE
bench_func = microbenchmark::microbenchmark(
hawk_process_org(x),
hawk_process_cache(x)
)
bench_func
#> Unit:milliseconds
#> expr min lq mean median uq max neval
#> hawk_process_org(x) 436.5349 465.9674 505.9606 481.4703 500.6652 894.7477 100
#> hawk_process_cache(x) 446.0499 454.9098 485.3830 468.6580 494.9457 799.0940 100
因此，根据建议，您会获得略微更好的结果。
标准操作在 C++ 中非常快（+、- 等）。
然而，exp 计算起来更复杂，所以更慢。
因此，如果我们想要提高一些性能，则更有可能预先计算 exp 计算。
这里，B(i) = \sum_{j < i}(x_i-x_j)exp^{-\beta(x_i - x_j)} 等同于 B(i) = \sum_{j < i}(x_i-x_j) / exp^{\beta x_i} * exp^{\beta x_j} 这样你就可以预先计算每个索引的 exp （并且也把依赖于 i 的那个放在循环）。通过重构它，您可以进行其他预计算。所以，我把之前的两个解决方案放在这里，然后是我的增量解决方案：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_org(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
for (int i = 1; i < n; i++) {
double temp = 0;
for (int j = 0; j <= i - 1; j++) {
temp += (x[i] - x[j]) * exp(-beta * (x[i] - x[j]));
}
B(i - 1) = temp;
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache(Rcpp::NumericVector x, double beta = 3) {
int n = x.size();
Rcpp::NumericVector B = Rcpp::no_init( n - 1);
double x_i;
for (int i = 1; i < n; ++i) {
double temp = 0;
x_i = x[i];
for (int j = 0; j <= i - 1; ++j) {
temp += (x_i - x[j]) * 1 / exp(beta * (x_i - x[j]));
}
B(i - 1) = temp;
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_2(Rcpp::NumericVector x,
double beta = 3) {
int i, j, n = x.size();
Rcpp::NumericVector B(n);
Rcpp::NumericVector x_exp = exp(beta * x);
double temp;
for (i = 1; i < n; i++) {
temp = 0;
for (j = 0; j < i; j++) {
temp += (x[i] - x[j]) * x_exp[j] / x_exp[i];
}
B[i] = temp;
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_3(Rcpp::NumericVector x,
double beta = 3) {
int i, j, n = x.size();
Rcpp::NumericVector B(n);
Rcpp::NumericVector x_exp = exp(beta * x);
double temp;
for (i = 1; i < n; i++) {
temp = 0;
for (j = 0; j < i; j++) {
temp += (x[i] - x[j]) * x_exp[j];
}
B[i] = temp / x_exp[i];
}
return B;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_4(Rcpp::NumericVector x,
double beta = 3) {
Rcpp::NumericVector exp_pre = exp(beta * x);
Rcpp::NumericVector exp_pre_cumsum = cumsum(exp_pre);
Rcpp::NumericVector x_exp_pre_cumsum = cumsum(x * exp_pre);
return (x * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_5(Rcpp::NumericVector x,
double beta = 3) {
int n = x.size();
NumericVector B(n);
double exp_pre, exp_pre_cumsum = 0, x_exp_pre_cumsum = 0;
for (int i = 0; i < n; i++) {
exp_pre = exp(beta * x[i]);
exp_pre_cumsum += exp_pre;
x_exp_pre_cumsum += x[i] * exp_pre;
B[i] = (x[i] * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
return B;
}
/*** R
set.seed(111)
x = rnorm(1e3)
all.equal(
hawk_process_org(x),
hawk_process_cache(x)
)
all.equal(
hawk_process_org(x),
hawk_process_cache_2(x)[-1]
)
all.equal(
hawk_process_org(x),
hawk_process_cache_3(x)[-1]
)
all.equal(
hawk_process_org(x),
hawk_process_cache_4(x)[-1]
)
all.equal(
hawk_process_org(x),
hawk_process_cache_5(x)[-1]
)
microbenchmark::microbenchmark(
hawk_process_org(x),
hawk_process_cache(x),
hawk_process_cache_2(x),
hawk_process_cache_3(x),
hawk_process_cache_4(x),
hawk_process_cache_5(x)
)
*/
x = rnorm(1e3) 的基准：
Unit: microseconds
expr min lq mean median uq max neval cld
hawk_process_org(x) 19801.686 20610.0365 21017.89339 20816.1385 21157.4900 25548.042 100 d
hawk_process_cache(x) 20506.903 21062.1370 21534.47944 21297.8710 21775.2995 26030.106 100 e
hawk_process_cache_2(x) 1895.809 2038.0105 2087.20696 2065.8220 2103.0695 3212.874 100 c
hawk_process_cache_3(x) 430.084 458.3915 494.09627 474.2840 503.0885 1580.282 100 b
hawk_process_cache_4(x) 50.657 55.2930 71.60536 57.6105 63.5700 1190.260 100 a
hawk_process_cache_5(x) 43.373 47.0155 60.43775 49.6640 55.6235 842.288 100 a
这比试图从小的优化中获得纳秒要有效得多，小的优化可能会使您的代码更难阅读。
不过，让我们试试@coatless 在我最后的解决方案中提出的优化：
// [[Rcpp::export]]
Rcpp::NumericVector hawk_process_cache_6(Rcpp::NumericVector x,
double beta = 3) {
int n = x.size();
NumericVector B = Rcpp::no_init(n);
double x_i, exp_pre, exp_pre_cumsum = 0, x_exp_pre_cumsum = 0;
for (int i = 0; i < n; ++i) {
x_i = x[i];
exp_pre = exp(beta * x_i);
exp_pre_cumsum += exp_pre;
x_exp_pre_cumsum += x_i * exp_pre;
B[i] = (x_i * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
return B;
}
x = rnorm(1e6) 的基准：
Unit: milliseconds
expr min lq mean median uq max neval cld
hawk_process_cache_5(x) 42.52886 43.53653 45.28427 44.46688 46.74129 57.38046 100 a
hawk_process_cache_6(x) 42.14778 43.19054 45.93252 44.28445 46.51052 153.30447 100 a
还是不太有说服力..
有趣的问题。在我的测试中，将这两个答案结合起来确实可以进一步提高性能（进一步降低基准）：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector hawk_process_cache_combined(NumericVector x,
double beta = 3) {
int n = x.size();
NumericVector B = Rcpp::no_init(n-1);
double exp_pre(exp(beta * x[0]));
double exp_pre_cumsum(exp_pre);
double x_exp_pre_cumsum(x[0] * exp_pre);
double x_i;
for (int i = 1; i < n; ++i) {
x_i = x[i];
exp_pre = exp(beta * x_i);
exp_pre_cumsum += exp_pre;
x_exp_pre_cumsum += x_i * exp_pre;
B[i-1] = (x_i * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
}
return B;
}
all.equal(
hawk_process_org(x),
hawk_process_cache_combined(x)
)
#> [1] TRUE
虽然原来的公式是 "embarrassingly parallel"，但此表达式不再是这种情况。但是，像 cumsum 这样的前缀扫描算法也可以并行化。像 ArrayFire provide interfaces to such algorithms using the GPU. Using RcppArrayFire 这样的库可以基于 F.Privé 的 hawk_process_cached_4:
// [[Rcpp::depends(RcppArrayFire)]]
#include <RcppArrayFire.h>
// [[Rcpp::export]]
af::array hawk_process_af(RcppArrayFire::typed_array<f32> x,
double beta = 3) {
af::array exp_pre = exp(beta * x);
af::array exp_pre_cumsum = af::accum(exp_pre);
af::array x_exp_pre_cumsum = af::accum(x * exp_pre);
af::array result = (x * exp_pre_cumsum - x_exp_pre_cumsum) / exp_pre;
return result(af::seq(1, af::end));
}
这里的结果并不完全相等，因为我的driver/card只支持单精度浮点数：
all.equal(
hawk_process_org(x),
hawk_process_af(x)
)
#> [1] "Mean relative difference: 3.437819e-07"
如果使用双精度，可以在上面写 f64 并获得相同的结果。现在进行基准测试：
set.seed(42)
x <- rnorm(1e3)
microbenchmark::microbenchmark(
hawk_process_af(x),
hawk_process_cache_combined(x),
hawk_process_cache_5(x)[-1]
)
#> Unit: microseconds
#> expr min lq mean median uq max neval
#> hawk_process_af(x) 245.281 277.4625 338.92232 298.5410 346.576 1030.045 100
#> hawk_process_cache_combined(x) 35.343 39.0120 43.69496 40.7770 45.264 84.242 100
#> hawk_process_cache_5(x)[-1] 52.408 57.8580 65.55799 60.5265 67.965 125.864 100
x <- rnorm(1e6)
microbenchmark::microbenchmark(
hawk_process_af(x),
hawk_process_cache_combined(x),
hawk_process_cache_5(x)[-1]
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> hawk_process_af(x) 27.54936 28.42794 30.93452 29.20025 32.40667 49.41888 100
#> hawk_process_cache_combined(x) 34.00380 36.84497 40.74862 39.03649 41.85902 111.51628 100
#> hawk_process_cache_5(x)[-1] 47.02501 53.24702 57.94747 55.35018 58.42097 130.89737 100
所以对于小向量，组合方法更快，而一旦卸载到 GPU 得到回报，时间更长。所有这一切都不是使用一些高功率 GPU，而是简单的板载图形：
RcppArrayFire::arrayfire_info()
#> ArrayFire v3.5.1 (OpenCL, 64-bit Linux, build 0a675e8)
#> [0] BEIGNET: Intel(R) HD Graphics Skylake ULT GT2, 4096 MB