加速简单的 R 代码（矢量化？-解网

问：

我有两个正整数向量，指定范围的开始和结束“位置”

starts <- sample(10^6,replace = T)
ends <- starts+sample(100:1000,length(starts),replace=T)

因此，它们指定了 1000000 个长度为 100 到 1000 个单位的范围。现在我想知道一个位置（正整数）被一个范围“覆盖”了多少次。为此，我这样做：

coverage <- integer(max(ends))
for(i in seq(length(starts))) {
      coverage[starts[i]:ends[i]] <- coverage[starts[i]:ends[i]] + 1 
}

但是由于for循环，它相对较慢。对于数十亿个范围，可能需要很长时间。我找不到矢量化此代码的方法。我可以拆分工作并使用多个 CPU，但速度提升微乎其微。apply、lapply 和其他元函数不会提高速度（正如预期的那样）。例如

coverage <- tabulate(unlist(Map(':', starts,ends)))

由于“地图”部分，速度也很慢。我担心这也需要更多的内存。

有什么想法吗？

R 性能

您的主要问题不是循环或，而是函数。矢量化它确实似乎是 SO 上一个非常频繁的请求......为可重现的例子和很好的尝试点赞。虽然你真的不需要创建一个最小的可重现的例子。在使用使用随机种子的函数创建数据集时，添加 a 始终是一种很好的做法。此外，显示所需的输出使其更易于阅读（尽管在您的情况下这没什么大不了的）。Map:10^6set.seed

0赞 jeanlain 5/5/2016

如果我在通话中替换，它仍然比简单的要慢得多。因此，与矢量化代码相比，Map 速度较慢。:+Mapstarts+ends

0赞 David Arenburg 5/5/2016

这不是关于具体的，而是关于你需要评估一个函数的次数+该函数的效率。如果您更换，在我的机器上大约需要一秒钟（相比之下，需要 7 秒）。计算函数 1e6 次还不错。在最坏的情况下，您可以在 Rcpp 中编写一个矢量化版本并完成。Map:+::

答：

3赞 Gergely Danyi 5/5/2016 #1

您可以保留从任何特定索引开始和结束的范围的计数，然后对这些索引的差值应用累积总和。

聚合从每个索引开始的范围数
聚合在每个索引之前的一个位置结束的范围数（如果包括在内）ends
计算净变化：count of starts - count of ends
遍历索引并累积汇总净变化。这将给出早于此索引开始且尚未在此索引处结束的数字范围。

“覆盖”数字等于每个指数的累积总和。

我尝试了这种方法，使用稀疏向量来减少内存使用。虽然使用正常向量可能会更快，但不确定。对于给定示例，它比循环方法快 5.7 倍。sparseVector

library(Matrix)

set.seed(123)

starts <- sample(10^6,replace = T)
ends <- starts+sample(100:1000,length(starts),replace=T)

v.cov <- NULL
fun1 <- function() {
  coverage <- integer(max(ends))
  for(i in seq(length(starts))) {
    coverage[starts[i]:ends[i]] <- coverage[starts[i]:ends[i]] + 1 
  }
  v.cov <<- coverage
}
# Testing "for loop" approach
system.time(fun1())
# user  system elapsed 
# 21.84    0.00   21.83 

v.sum <- NULL
fun2 <- function() {      
  # 1. Aggregate the number of ranges that start at each index
  t.starts <- table(starts)
  i.starts <- strtoi(names(t.starts))
  x.starts <- as.vector(t.starts)
  sv.starts <- sparseVector(x=x.starts, i=i.starts, length=max(ends)+1)  # to match length of sv.ends below
  # 2. Aggregate the number of ranges that end at one position before each index
  t.ends <- table(ends)
  i.ends <- strtoi(names(t.ends))+1  # because "ends" are inclusive 
  x.ends <- as.vector(t.ends)
  sv.ends <- sparseVector(x=x.ends, i=i.ends, length=max(ends)+1)

  sv.diff <- sv.starts - sv.ends
  v.sum <<- cumsum(sv.diff)[1:max(ends)]  # drop last element
}
# Testing "cumulative sum" approach
system.time(fun2())
# user  system elapsed 
# 3.828   0.000   3.823

identical(v.cov, v.sum)
# TRUE

此外，可能有一种比使用构造函数更好的方法来提取 x 和 i，这可能会进一步提高速度。sparseVectortablestrtoi(names(x))

编辑

避免改用 1 列strtoisparseMatrix

v.sum.mat <- NULL
fun3 <- function() {
  v.ones <- rep(1, length(starts))
  m.starts <- sparseMatrix(i=starts, j=v.ones, x=v.ones, dims=c(max(ends)+1,1))
  m.ends <- sparseMatrix(i=ends+1, j=v.ones, x=v.ones, dims=c(max(ends)+1,1))
  m.diff <- m.starts - m.ends
  v.sum.mat <<- cumsum(m.diff[,1])[1:max(ends)]
}
# Testing "cumulative sum" approach using matrix
system.time(fun3())
#   user  system elapsed 
#  0.456   0.028   0.486 

identical(v.cov, v.sum.mat)
# TRUE

编辑 2 - 超快，超短

根据 @alexis_laz 的评论，谢谢！

fun4 <- function() {
  cumsum(tabulate(starts, max(ends) + 1L) - tabulate(ends + 1L, max(ends) + 1L))[1:max(ends)]
}
system.time(v.sum.tab <- fun4())
# user  system elapsed 
# 0.040   0.000   0.041 

identical(as.integer(v.cov), v.sum.tab)
# TRUE

加速简单的 R 代码（矢量化？

Speed up simple R code (vectorize?)

评论

评论