使用环绕优化跨越三个连续索引的公式的性能答案

【问题标题】：Optimize performance of a formula spanning three consecutive indices, with wraparound使用环绕优化跨越三个连续索引的公式的性能
【发布时间】：2017-09-23 04:30:30
【问题描述】：

我想优化这个公式的实现。

这是公式：

x 是一个值数组。 i 从 1 变为 N，其中 N > 2400000。对于i=0，i-1 是最后一个元素，对于i=lastElement，i+1 是第一个元素。这是我写的代码：

   x <- 1:2400000
   re <- array(data=NA, dim = NROW(x))
   lastIndex = NROW(x)
   for(i in 1:lastIndex){
      if (i==1) {
        re[i] = x[i]*x[i] - x[lastIndex]*x[i+1]
      } else if(i==lastIndex) {
        re[i] = x[i]*x[i] - x[i-1]*x[1]
      } else {
        re[i] = x[i]*x[i] - x[i-1]*x[i+1]  
      }
    }

可以由apply 在 R 中完成吗？

【问题讨论】：

“优化它”是为了什么？紧凑性（见我的回答）性能？（使用data.table）
优化 w.r.t.性能。
Chetan，那么你必须编辑标题才能说出来。 “优化”本质上是模棱两可的。

标签： r performance formula apply indices

【解决方案1】：

您的公式的 lapply 实现如下所示：

x <- c(1:2400000) 
last <- length(x)

re <- lapply(x, function(i) {
    if(i == 1) {
        x[i]*x[i] - x[last]*x[i+1]
    } else if (i == last) {
        x[i]*x[i] - x[i-1]*x[1]
    } else {
        x[i]*x[i] - x[i-1]*x[i+1]  
    }
}) 

re <- unlist(re)

lapply 将返回一个列表，因此转换为向量是使用 unlist() 完成的

【讨论】：

使用sapply 而不是lapply，它不返回列表，而是返回向量/矩阵。甚至vapply 提前知道输出的大小和类型

【解决方案2】：

1) 您可以通过用最后一行和第一行的副本填充数组 x 的开头和结尾来避免计算中的所有特殊情况；像这样：

N <- NROW(x)
x <- rbind(x[N], x, x[1]) # pad start and end to give wraparound 

re <- lapply(2:N, function(i) { x[i]*x[i] - x[i-1]*x[i+1] } )
#re <- unlist(re) as andbov wrote

# and remember not to use all of x, just x[2:N], elsewhere

2) 直接矢量化，正如@Dason 的回答：

# Do the padding trick on x , then
x[2:N]^2 - x[1:N-1]*x[3:N+1]

3) 如果性能很重要，我怀疑在 i 上使用 data.table 或 for-loop 会更快，因为它引用了三个连续的行。

4) 要获得更多性能，use byte-compiling

5) 如果您需要更快的速度，使用 Rcpp 扩展（C++ 底层）How to use Rcpp to speed up a for loop?

请参阅我引用的那些问题，了解使用 lineprof 和微基准测试找出瓶颈所在的良好示例。

【讨论】：

在lapply 里面不应该是2:N 而不是x[2:N]？此外，这不是性能效率，需要很长时间才能运行。
我喜欢填充部分。聪明的举动:)
@Chetan：将一些随机种子数据添加到您的问题详细信息中，这样我们就可以实际运行苹果对苹果的比较。 “运行需要很长时间” 并不具体，我们其他任何人也无法验证它。对于N>240万，取一个实际值。我假设您没有超出内存限制；如果你是，所有的赌注都没有。
当然，2:N，而不是 x[2:N]，无论如何，代码的意图很明确。
致反对者：这方面做了很多工作，所以请告诉我您认为需要改进的地方。

【解决方案3】：

我们可以为此使用直接向量化

# Make fake data
x <- 1:10
n <- length(x)
# create vectors for the plus/minus indices
xminus1 <- c(x[n], x[-n])
xplus1 <- c(x[-1], x[1])

# Use direct vectorization to get re
re <- x^2 - xminus1*xplus1

【讨论】：

太棒了！谢谢达森:)
这是创建一个非常大的向量/数组的三个副本。您可以使用填充技巧避免副本，然后 x[2:N]^2 - x[1:N-1]*x[3:N+1]
@Chetan：这需要 3 倍的内存。如果 x 很大，那么当你用完内存时，它会降低性能。

【解决方案4】：

如果每个x[i] 都等于i 那么你可以做一些数学运算：
xi^2 - (xi-1)*(xi+1) = 1
所以结果的所有元素都是1（只有第一个和最后一个不是1）。
结果是：

c(1-2*N, rep(1, N-2), N*N-(N-1))

在一般情况下（x 中的任意值）您可以这样做（如 Dason 的回答）：

x*x - c(x[N], x[-N])*c(x[-1], x[1])

这是来自zoo 的rollapply() 的解决方案：

library("zoo")
rollapply(c(x[length(x)],x, x[1]), width=3, function(x) x[2]^2 - x[1]*x[3]) # or:
rollapply(c(tail(x,1), x, x[1]), width=3, function(x) x[2]^2 - x[1]*x[3])

这是基准：

library("microbenchmark")
library("zoo")

N <- 10000
x <- 1:N

microbenchmark(
  math=c(1-2*N, rep(1, N-2), N*N-(N-1)), # for the data from the question
  vect.i=x*x - c(x[N], x[-N])*c(x[-1], x[1]), # general data
  roll.i=rollapply(c(x[length(x)],x, x[1]), width=3, function(x) x[2]^2 - x[1]*x[3]), # or:
  roll.tail=rollapply(c(tail(x,1), x, x[1]), width=3, function(x) x[2]^2 - x[1]*x[3])
)
# Unit: microseconds
#      expr       min         lq        mean     median         uq        max neval cld
#      math    33.613    34.4950    76.18809    36.9130    38.0355   2002.152   100  a 
#    vect.i   188.928   192.5315   732.50725   197.1955   198.5245  51649.652   100  a 
#    roll.i 56748.920 62217.2550 67666.66315 68195.5085 71214.9785 109195.049   100   b
# roll.tail 57661.835 63855.7060 68815.91001 67315.5425 71339.6045 119428.718   100   b

【讨论】：