加速重复的函数调用答案

【问题标题】：Speeding up a repeated function call加速重复的函数调用
【发布时间】：2013-01-24 07:05:05
【问题描述】：

我正在尝试计算以下值：

1/N * sum[i=0 to N-1]( log(abs(r_i - 2 * r_i * x_i)) )

其中 x_i 是递归计算的：

x_{i+1} = r_i * x_i * (1 - x_i)

所有r_is 都被给出（尽管它们随着i 而变化），x_0 被给出。（据我所知，没有任何棘手的数学方法可以将此计算简化为非迭代公式以加快速度）。

我的问题是它非常慢，我想知道一些外部观点是否可以帮助我加快速度。

# x0: a scalar. rs: a numeric vector, length N
# N: typically ~5000
f <- function (x0, rs, N) {
    lambda <- 0                                                                 
    x <- x0                                                                     
    for (i in 1:N) {                                                            
        r <- rs[i]                                                              
        rx <- r * x                                                             
        lambda <- lambda + log(abs(r - 2 * rx))                                 
        # calculate the next x value
        x <- rx - rx * x                                                        
    }                                                                           
    return(lambda / N)                                                          
}

现在这个函数本身就相当快了，但是我想调用它大约 4,000,000 次（对于 2000 x 2000 矩阵中的每个单元格一次），每个都有不同的 @987654327 @向量。

但如果我只调用它 2500 次（N=1000），它需要大约 25 秒，使用以下配置文件：

      self.time self.pct total.time total.pct
"f"       19.98    81.22      24.60    100.00
"*"        2.00     8.13       2.00      8.13
"-"        1.32     5.37       1.32      5.37
"+"        0.70     2.85       0.70      2.85
"abs"      0.56     2.28       0.56      2.28
":"        0.04     0.16       0.04      0.16

有谁知道我可以如何加快速度？看起来乘法需要一段时间，但我已经预先缓存了所有重复的乘法。

我还尝试利用sum( log(stuff(i)) ) 与log(prod(stuff(i)) 相同来减少对log 和abs 的调用，但结果证明这是不可行的，因为stuff 是长度为N 的向量（以千计）并且典型值至少为 1，因此 prod(stuff) 最终成为 Inf 到 R。

【问题讨论】：

标签： r vectorization

【解决方案1】：

在我看来，瓶颈是函数中的for 循环。

我用 Rcpp 重写如下：

# x0: a scalar. rs: a numeric vector, length N
# N: typically ~5000
x0 <- runif(1)
N <- 5000
rs <- rnorm(5000)
f <- function (x0, rs, N) {
    lambda <- 0                                                                 
    x <- x0                                                                     
    for (i in 1:N) {                                                            
        r <- rs[i]                                                              
        rx <- r * x                                                             
        lambda <- lambda + log(abs(r - 2 * rx))                                 
        # calculate the next x value
        x <- rx - rx * x                                                        
    }                                                                           
    return(lambda / N)                                                          
}

library(inline)
library(Rcpp)
f1 <- cxxfunction(sig=c(Rx0="numeric", Rrs="numeric"), plugin="Rcpp", body='
  double x0 = as<double>(Rx0);
  NumericVector rs(Rrs);
  int N = rs.size();
  double lambda = 0, x = x0, r, rx;
  for(int i = 0;i < N;i++) {
    r = rs[i];
    rx = r * x;
    lambda = lambda + log( fabs(r - 2 * rx) );
    x = rx - rx * x;
  }
  lambda /= N;
  return wrap(lambda);
  ')
f(x0, rs, N)
f1(x0, rs)

library(rbenchmark)

benchmark(f(x0, rs, N), f1(x0, rs))

在我上次的测试中，f1 比 f 快 140 倍。

【讨论】：