加速循环中指数随机变量的重复生成答案

【问题标题】：Speeding up repeated generation of Exponential random variables in loop加速循环中指数随机变量的重复生成
【发布时间】：2017-08-25 05:37:22
【问题描述】：

我正在实现一个算法，作为其中的一部分，我需要生成指数随机变量。不幸的是，我无法真正避免循环，因为每个生成的随机变量都取决于前一个，所以我认为矢量化是不可能的。我围绕生成进行了一些计算，但瓶颈（目前）是生成。在这一点上，我假设 N 会很大（N >= 1,000,000）。

下面是一些示例代码：

N <- 1e7
#Preallocate
x <- rep(0, times=N)

#Set a starting seed
x[1] <- runif(1)

for(i in 2:N) {
    #Do some calculations
    x[i] <- x[i-1] + rexp(1, x[i-1])  #Bottleneck
    #Do some more calculations
}

如何加快速度？我已经尝试在 Rcpp 中实现，但在这种情况下它似乎并没有做太多。是否有另一种聪明的方法可以绕过每次迭代中的 rexp() 调用？

【问题讨论】：

标签： r performance random rcpp

【解决方案1】：

蛮力 Rcpp 解决方案：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector genExp(int N) {

  NumericVector res(N);
  double prev;

  res[0] = unif_rand();
  for (int i = 1; i < N; i++) {
    prev = res[i-1];
    res[i] = prev + exp_rand() / prev;
  }

  return res;
}

使用N = 1e6 进行微基准测试：

Unit: milliseconds
             expr        min         lq       mean     median         uq        max neval
    draw_up_front  167.17031  168.57345  170.62292  170.18072  171.73782  175.46868    20
 draw_one_at_time 1415.01898 1465.57139 1510.81220 1502.15753 1550.07829 1623.70831    20
             rcpp   28.25466   29.33682   33.52528   29.89636   30.74908   94.38009    20

N = 1e7:

Unit: milliseconds
          expr      min       lq      mean    median        uq       max neval
 draw_up_front 1698.730 1708.739 1737.8633 1716.1345 1752.3276 1923.3940    20
          rcpp  297.142  319.794  338.6755  327.6626  364.6308  398.1554    20

【讨论】：

谢谢。我已经使用这个答案来帮助加快我的代码速度，但是真正的进步是由于上面 Gregor 的答案，所以接受了那个。

【解决方案2】：

我们可以使用 if X ~ Exp(λ) then kX ~ Exp(λ/k) (source: Wikipedia) 来加速代码。这样我们就可以在前面使用rate = 1 进行所有随机绘制，然后在循环内分割以适当地缩放它们。

draws = rexp(N, rate = 1)
x <- rep(0, times = N)
x[1] <- runif(1)
for(i in 2:N) {
    #Do some calculations
    x[i] <- x[i-1] + draws[i] / x[i-1] 
    #Do some more calculations
}

具有N = 1e6 值的微基准测试表明这大约快 14 倍：

N <- 1e6
draws = rexp(N, rate = 1)
x <- rep(0, times = N)
x[1] <- runif(1)

microbenchmark::microbenchmark(
    draw_up_front = {
        draws = rexp(N, rate = 1)
        for (i in 2:N)
            x[i] <- x[i - 1] + draws[i] / x[i - 1]
    },
    draw_one_at_time = {
        for (i in 2:N)
            x[i] <- x[i - 1] + rexp(1, x[i - 1])
    },
    times = 10
)
# Unit: milliseconds
#              expr       min        lq      mean    median        uq       max neval cld
#     draw_up_front  153.9547  156.6552  159.9622  160.1901  161.9803  167.2831    10  a 
#  draw_one_at_time 2207.1997 2212.0460 2280.1265 2236.5197 2332.9913 2478.5104    10   b

【讨论】：

嗯...有趣。我会试试这个，看看效果如何。感谢您的回复。
如果你没有把所有的函数都放在基准测试中，那就是作弊。