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

Question

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

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)))

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

有什么想法吗？

Answer 1

您可以对在任何特定索引处开始和结束的范围进行计数，然后对它们的差值应用累积总和。

1. 聚合从每个索引开始的范围数
2. 聚合在每个索引前一个位置结束的范围数（如果包含ends）
3. 计算净变化：count of starts - count of ends
4. 循环索引并累计净变化。这将给出早于该索引开始但尚未在该索引处结束的数字范围。

“覆盖”数等于每个索引处的累积总和。

我尝试使用稀疏向量来减少内存使用这种方法。尽管使用法线向量可能会更快，但不确定。 使用sparseVector，它比给定示例的循环方法快 5.7 倍。

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

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

编辑

避免 strtoi 使用 1 列 sparseMatrix 代替

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

EDIT 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