R语言双层for循环区间计数脚本运行速度优化方案

问题背景 #

我编写了一段基于双层for循环的R脚本，希望优化代码缩短运行时间。
以下是简化后的可复现数据，以及我在自有数据上实际使用的代码：

• nuc：存储位点位置的向量，示例中用nuc2=100:200演示，实际数据长度达9304567
• tel：存储坐标对aa（区间起点）、bb（区间终点）和对应ID的数据框，当前实际数据有53组坐标，后续会增长至数百组
  计算目标为统计每个位置落在各组aa、bb坐标区间内的次数，例如位置111共落在G、I、J三组坐标对应的区间中。

示例数据与原始实现代码 #

# 构造示例数据
tel=data.frame(aa=c(153,113,163,117,193,162,110,109,186,103),
               bb=c(189,176,185,130,200,189,156,123,198,189),
               ID=c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J"))

# 查看示例tel数据
> tel
    aa  bb ID
1  153 189  A
2  113 176  B
3  163 185  C
4  117 130  D
5  193 200  E
6  162 189  F
7  110 156  G
8  109 123  H
9  186 198  I
10 103 189  J

nuc2=100:200

# 原始双层for循环实现
count_occ=0
count_occ_int=NULL
count_occ_fin=NULL

for (j in 1:length(nuc2)){
  for (i in 1:nrow(tel)) {
    
    if (nuc2[j]< tel$bb[i] & nuc2[j]>tel$aa[i]) 
    {count_occ=count_occ+1}
  }
  count_occ_int=count_occ
  count_occ_fin=c(count_occ_fin,count_occ_int)
  count_occ=0
}

nuc_occ=data.frame(nuc=nuc2, occ=count_occ_fin)

# 查看前20行计算结果
> head(nuc_occ,20)
   nuc occ
1  100   0
2  101   0
3  102   0
4  103   0
5  104   1
6  105   1
7  106   1
8  107   1
9  108   1
10 109   1
11 110   2
12 111   3
13 112   3
14 113   3
15 114   4
16 115   4
17 116   4
18 117   4
19 118   5
20 119   5

实际数据场景下，现有代码运行耗时超过60小时。此前考虑过使用apply函数改写，但不确定如何实现原双层for循环的计算逻辑。

优化方案 #

按改造成本从低到高、性能提升幅度从小到大排序，可按需选择：

方案1：向量化运算替换R层循环（零依赖，改造成本极低） #

R的向量化运算为底层C实现，比R层面的for循环快3~4个数量级，无需嵌套循环即可完成计算：

# 核心逻辑：构造逻辑判断矩阵，按行求和得到每个位点的区间命中次数
count_occ_fin <- rowSums(outer(nuc2, tel$aa, `>`) & outer(nuc2, tel$bb, `<`))
nuc_occ <- data.frame(nuc = nuc2, occ = count_occ_fin)

• 注意：如果位点长度达900万、区间数增长到数百，逻辑矩阵会占用约1GB内存，内存足够的前提下，该方案可将运行时间从60小时压缩到10秒以内。

方案2：data.table非等连接（内存友好，适合大数据量） #

如果内存不足以存储大逻辑矩阵，可使用data.table的非等连接实现区间计数，无需构造大矩阵，内存占用极低，速度与向量化方案接近：

library(data.table)
setDT(tel)
nuc_dt <- data.table(nuc = nuc2)
# 非等连接统计命中次数
nuc_occ <- tel[nuc_dt, on = .(aa < nuc, bb > nuc), .N, by = .EACHI]
setnames(nuc_occ, c("nuc", "drop_col", "occ"))
nuc_occ[, drop_col := NULL]

该方案在位点长度千万级、区间数上千的场景下，内存占用仅数百MB，运行时间在秒级。

方案3：差分法（性能最优，内存占用极低） #

差分法是区间计数场景的专用最优算法，仅需两次线性扫描即可得到结果，不受位点长度、区间数量增长的影响：

# 对齐位点索引范围
min_nuc <- min(nuc2)
max_nuc <- max(nuc2)
count_vec <- integer(length = max_nuc - min_nuc + 3) # 预留冗余位避免索引越界

# 遍历区间做差分标记（原逻辑为开区间，有效范围是aa+1 到 bb-1）
for (i in seq_len(nrow(tel))) {
  count_vec[tel$aa[i] - min_nuc + 2] <- count_vec[tel$aa[i] - min_nuc + 2] + 1
  count_vec[tel$bb[i] - min_nuc + 1] <- count_vec[tel$bb[i] - min_nuc + 1] - 1
}
# 累计求和得到最终计数
count_occ_fin <- cumsum(count_vec)[seq_along(nuc2)]
nuc_occ <- data.frame(nuc = nuc2, occ = count_occ_fin)

• 该方案在位点长度千万级、区间数上万的场景下，运行时间不超过1秒，内存占用仅几MB，是当前场景的最优选择。

不建议用apply系列函数改写原始循环，apply本质仍是R层面的循环，不会带来本质性能提升，只有将R层循环替换为底层实现的向量化逻辑或专用算法，才能实现数百到数千倍的速度提升。

内容的提问来源于stack exchange，提问作者Aurelia Kurtis