矢量化依赖于先前元素的产品计算？

Question

我正在尝试加速/矢量化时间序列中的一些计算.我可以在for循环中矢量化计算,这可能取决于早期迭代的结果吗？例如:

z <- c(1,1,0,0,0,0)
zi <- 2:6
for  (i in zi) {z[i] <- ifelse (z[i-1]== 1, 1, 0) }

使用前面步骤中更新的z [i]值:

> z
[1] 1 1 1 1 1 1

在我的努力中矢量化这个

z <- c(1,1,0,0,0,0)
z[zi] <- ifelse( z[zi-1] == 1, 1, 0)

逐个元素的操作不使用操作中更新的结果:

> z
[1] 1 1 1 0 0 0

因此,这种向量化操作以"并行"而非迭代方式运行.有没有办法我可以写/向量化这个来获得for循环的结果？

Answer 1

ifelse是矢量化的,如果你在for循环中一次在一个元素上使用它会有一点点惩罚.在您的示例中,您可以通过使用if而不是来获得相当好的加速ifelse.

fun1 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- ifelse(z[i-1]==1, 1, 0)
  }
  z
}

fun2 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- if(z[i-1]==1) 1 else 0
  }
  z
}

z <- c(1,1,0,0,0,0)
identical(fun1(z),fun2(z))
# [1] TRUE
system.time(replicate(10000, fun1(z)))
#   user  system elapsed 
#   1.13    0.00    1.32
system.time(replicate(10000, fun2(z)))
#   user  system elapsed 
#   0.27    0.00    0.26 

fun2通过编译,您可以获得额外的速度提升.

library(compiler)
cfun2 <- cmpfun(fun2)
system.time(replicate(10000, cfun2(z)))
#   user  system elapsed 
#   0.11    0.00    0.11

所以没有矢量化就有10倍的加速.正如其他人所说(有些人已经说明),有一些方法可以对你的例子进行矢量化,但这可能无法转化为你的实际问题.希望这足以适用.

filter如果您可以根据自回归或移动平均过程弄清楚如何表达问题,该功能对您也很有用.

Answer 2

这是一个很好的简单例子,Rcpp可以发光.

那么让我们首先重新编写函数1和2以及它们编译的对应函数:

library(inline)
library(rbenchmark)
library(compiler)

fun1 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- ifelse(z[i-1]==1, 1, 0)
  }
  z
}
fun1c <- cmpfun(fun1)


fun2 <- function(z) {
  for(i in 2:NROW(z)) {
    z[i] <- if(z[i-1]==1) 1 else 0
  }
  z
}
fun2c <- cmpfun(fun2)

我们非常容易地编写Rcpp变体:

funRcpp <- cxxfunction(signature(zs="numeric"), plugin="Rcpp", body="
  Rcpp::NumericVector z = Rcpp::NumericVector(zs);
  int n = z.size();
  for (int i=1; i<n; i++) {
    z[i] = (z[i-1]==1.0 ? 1.0 : 0.0);
  }
  return(z);
")

这使用内联包来动态编译,加载和链接五线程.

现在我们可以定义我们的测试日期,我们比原来的更长一些(因为只运行原始的太少次导致不可测量的时间):

R> z <- rep(c(1,1,0,0,0,0), 100)
R> identical(fun1(z),fun2(z),fun1c(z),fun2c(z),funRcpp(z))
[1] TRUE
R> 

所有答案都被视为相同.

最后,我们可以测试:

R> res <- benchmark(fun1(z), fun2(z),
+                  fun1c(z), fun2c(z),
+                  funRcpp(z),
+                  columns=c("test", "replications", "elapsed", 
+                            "relative", "user.self", "sys.self"),
+                  order="relative",
+                  replications=1000)
R> print(res)
        test replications elapsed relative user.self sys.self
5 funRcpp(z)         1000   0.005      1.0      0.01        0
4   fun2c(z)         1000   0.466     93.2      0.46        0
2    fun2(z)         1000   1.918    383.6      1.92        0
3   fun1c(z)         1000  10.865   2173.0     10.86        0
1    fun1(z)         1000  12.480   2496.0     12.47        0

对于最好的R版本,编译版本赢得了近400倍,而对于其字节编译版本,编译版本几乎达到100.对于函数1,字节编译的重要性要小得多,并且这两种变体都使C++的性能远远超过两千.

编写C++版本花了大约一分钟.速度提升表明它花了一分钟.

为了比较,这里是更经常调用的原始短向量的结果:

R> z <- c(1,1,0,0,0,0)
R> res2 <- benchmark(fun1(z), fun2(z),
+                  fun1c(z), fun2c(z),
+                  funRcpp(z),
+                  columns=c("test", "replications", 
+                            "elapsed", "relative", "user.self", "sys.self"),
+                  order="relative",
+                  replications=10000)
R> print(res2)
        test replications elapsed  relative user.self sys.self
5 funRcpp(z)        10000   0.046  1.000000      0.04        0
4   fun2c(z)        10000   0.132  2.869565      0.13        0
2    fun2(z)        10000   0.271  5.891304      0.27        0
3   fun1c(z)        10000   1.045 22.717391      1.05        0
1    fun1(z)        10000   1.202 26.130435      1.20        0

定性排名没有变化:Rcpp版本占主导地位,功能2 排名第二.字节编译版本的速度大约是普通R版本的两倍,但仍然比C++版本快三倍.并且相对差异较小:相对而言,函数调用开销较少,实际循环更重要:C++在较长向量中的实际循环操作中获得更大的优势.这是一个重要的结果,因为它表明更多真实大小的数据,编译版本可以获得更大的好处.

编辑纠正代码示例中的两个小疏忽.并再次编辑,感谢Josh捕获相对于fun2c的设置错误.

Answer 3

我认为这是作弊而不是概括,但是:根据你上面的规则,向量中出现1将使所有后续元素1(通过递归:z[i]如果z[i-1]等于1 则将1设置为1;因此z[i]将设置为1 z[i-2]等于1;等等).根据您真正想要做的事情,如果您仔细考虑,可能会有这样的递归解决方案......

z <- c(1,1,0,0,0,0)
first1 <- min(which(z==1))
z[seq_along(z)>first1] <- 1

编辑:这是错的,但是我不肯承认我的错误.基于一点点的游戏(并且思路较少),我认为这种递归的实际解决方案更加对称,甚至更简单:

rep(z[1],length(z))

测试用例:

z <- c(1,1,0,0,0,0)
z <- c(0,1,1,0,0,0)
z <- c(0,0,1,0,0,0)