0

我正在尝试优化以下代码。

dim <- c(10000,100)

m <- matrix(sample(0:10, prod(dim), replace = TRUE), nrow = dim[1], ncol = dim[2])

system.time({

  output <- matrix(0, nrow = dim[1], ncol = dim[2])

  for (i in 1:dim[1]){        
    output[i,1] <- m[i,1]       
    for (j in 2:dim[2]){          
      output[i,j] <- output[i, j-1] * 0.5 + m[i,j]          
    }        
  }
})

从概念上讲，它与简单的累积和非常相似：

system.time({

  output <- matrix(0, nrow = dim[1], ncol = dim[2])

  for (i in 1:dim[1]){       
      output[i,] <- cumsum(m[i,])        
  }
})

问题是，代码的第一部分慢了大约 100 倍。有没有办法构建一个定制版本的 cumsum() 来解决问题？

4

2 回答 2

2

您的情况与生成系数为 0.5 的 AR(1) 模型完全相同。您可以使用该filter函数生成数据。filter还支持高阶递归、卷积或它们的混合（想想 ARMA 模型）。你可以看看convolve其他卷积。此外，您可以编译代码以加快循环速度。在我的代码中，编译循环和未编译循环代码分别比过滤器慢 111 倍和 162 倍。

library(compiler)
library(rbenchmark)

CustomCumsum<-function(x,alpha){
out<-x[1]
for(i in 2:length(x))
    out[i] <- out[i-1]*alpha+x[i]
out
}

compiledCustomCumsum<-cmpfun(CustomCumsum)

FilterCustomCumsum<-function(x,alpha) as.numeric(filter(x,alpha, method = "recursive"))

x<-rnorm(1000)
# Test whether they are the same
identical(compiledCustomCumsum(x,0.5) , FilterCustomCumsum(x,0.5) )

benchmark(
CustomCumsum=CustomCumsum(x,0.5),compiledCustomCumsum=compiledCustomCumsum(x,0.5),          FilterCustomCumsum=FilterCustomCumsum(x,0.5)
)

输出：

                  test replications elapsed relative user.self sys.self user.child sys.child
2 compiledCustomCumsum          100    8.89  111.125      8.78     0.01         NA        NA
1         CustomCumsum          100   13.02  162.750     11.84     0.50         NA        NA
3   FilterCustomCumsum          100    0.08    1.000      0.08     0.00         NA        NA
于 2013-02-11T16:11:58.857 回答
2

用 C 编写我的自定义累积函数确实比其他所有方法都快：

sign <- signature(x="numeric", n="integer", d="numeric")

code <- "
  for (int i=1; i < *n; i++) {
    x[i] = x[i-1]*d[0] + x[i];
  }"

c_fn <- cfunction(sign,
                  code,
                  convention=".C"
)

CCustomCumsum <- function(vector, decay){
  c_fn(x=vector, n=length(vector), d=decay)$x
}

运行 Julian 的基准测试：

x<-rnorm(1000)
# Test whether they are the same
identical(compiledCustomCumsum(x,0.5) , FilterCustomCumsum(x,0.5) )

benchmark(
  CustomCumsum=CustomCumsum(x,0.5),
  compiledCustomCumsum=compiledCustomCumsum(x,0.5),
  FilterCustomCumsum=FilterCustomCumsum(x,0.5),
  CCustomCumsum = CCustomCumsum(x, 0.5)  
)

，我得到：

                  test replications elapsed relative user.self sys.self user.child sys.child
4        CCustomCumsum          100   0.002      1.0     0.002    0.000          0         0
2 compiledCustomCumsum          100   0.631    315.5     0.536    0.095          0         0
1         CustomCumsum          100   0.931    465.5     0.882    0.046          0         0
3   FilterCustomCumsum          100   0.036     18.0     0.033    0.003          0         0
于 2013-02-11T19:50:16.887 回答