匹配两个具有容差的非常大的向量（快速！但节省工作空间）

Posted 2023-02-21

【中文标题】匹配两个具有容差的非常大的向量（快速！但节省工作空间）

【英文标题】：Matching two very very large vectors with tolerance (fast! but working space sparing)

【发布时间】：2018-04-08 00:39:54

【问题描述】：

考虑我有两个向量。一个是包含所有感兴趣值的参考向量/列表和一个可以包含任何可能值的样本向量。现在我想在参考列表中找到我的样本的匹配项，具有一定的容差，该容差不是固定的，并且依赖于向量内的比较值：

matches: abs(((referencelist - sample[i])/sample[i])*10^6)) < 0.5

将两个向量四舍五入是不行的！

例如考虑：

referencelist <- read.table(header=TRUE, text="value  name
154.00312  A
154.07685  B
154.21452  C
154.49545  D
156.77310  E
156.83991  F
159.02992  G
159.65553  H
159.93843  I")

sample <- c(154.00315, 159.02991, 154.07688, 156.77312)

所以我得到了结果：

    name value      reference
1    A   154.00315  154.00312
2    G   159.02991  159.02992
3    B   154.07688  154.07685
4    E   156.77312  156.77310

我能做的是使用例如像

这样的外部函数

myDist <- outer(referencelist, sample, FUN=function(x, y) abs(((x - y)/y)*10^6))
matches <- which(myDist < 0.5, arr.ind=TRUE)
data.frame(name = referencelist$name[matches[, 1]], value=sample[matches[, 2]])

或者我可以使用for() 循环。

但我的特殊问题是，参考向量大约有 1*10^12 个条目，而我的样本向量大约有 1*10^7。因此，通过使用 outer()，我可以轻松破坏所有工作空间限制，并通过使用 for() 或链式 for() 循环，这将需要数天/数周才能完成。

有没有人知道如何在 R 中快速做到这一点，仍然很精确，但在计算机上工作消耗最大。 64 GB 内存？

感谢您的帮助！

祝福

【问题讨论】：

参考向量是否已排序？如果是这样，您可以使用二分法来找到最接近的数字。如果没有 - 排序！另外，它是如何存储的？它太大而无法放入您的 RAM。是单个文件还是跨越多个文件？

在您给出的示例中，所有值似乎都在一个足够窄的范围内（154、160）。我们可以通过计算(x - y) / 154 而不是(x - y) / y 来简化问题吗？

@john Coleman：是的，它按升序排序。你能举例说明什么是二分法吗？以前没听过，不好意思。参考向量适合我的 RAM，然后我有大约 64 GB 的可用空间。

@aurele 我不知道 100% 是什么意思，但如果我理解正确，那么不，我们不能，因为两个向量的范围从大约 100 到 3000。

我的意思是“二分搜索”。看基础R函数findInterval()

【参考方案1】：

使用data.table（并从@eddi 的binary search 复制粘贴（也称为二分法，参见@John Coleman 的评论））：

library(data.table)

dt <- as.data.table(referencelist)
setattr(dt, "sorted", "value")

tol <- 0.5
dt2 <- dt[J(sample), .(.I, ref = value, name), roll = "nearest", by = .EACHI]
dt2[, diff := abs(ref - value) / value * 1e6]
dt2[diff <= tol]

#       value I      ref name       diff
# 1: 154.0032 1 154.0031    A 0.19480121
# 2: 159.0299 7 159.0299    G 0.06288125
# 3: 154.0769 2 154.0769    B 0.19470799
# 4: 156.7731 5 156.7731    E 0.12757289

我没有对内存使用和执行时间进行基准测试，但data.table 以在这两个方面都非常出色而闻名。如果它不适合您，请说出来，也许我会尝试对事物进行基准测试。

注意：我对data.table 的使用非常幼稚。

还有一个使用findInterval 的解决方案就在下面：https://***.com/a/29552922/6197649，但我预计它的性能会更差（同样：需要基准测试）。

【讨论】：

非常感谢，这听起来很不错。我会尽快对其进行测试，但我很遗憾地说，由于向量的大小，这可能需要几个小时/几天。很快回来，再次感谢！

只是一个简短的问题，sample 需要是单个向量还是可以是数据框或矩阵，可以选择包含其中值的列样本？

@JmO 线性搜索是O(n)。对于n = 10^12，这是令人望而却步的，尤其是当您想这样做10^7 次时。另一方面，二分查找是O(log_2(n))。 10^12 的以 2 为底的对数约为 40。请注意，每次搜索 40 步而不是 10^12 是 250 亿倍的加速。不需要几天甚至几个小时。只需几分钟就足够了。

非常感谢你们！我今天将执行它。对于上面的答案：由于我没有真正使用过 data.table 这是否只给我最近的匹配或公差范围内的每个匹配？：考虑我在参考列表中有一个 F = 154.0033。现在对于示例中的 156.0032，在公差范围内的参考列表中有 A 和 F，我希望将两者都放在单独的行中的结果列表中。这种方法有可能吗，还是只给我最近的？

像 roll="nearest" +- 2 这样的解决方案就足够了。有什么方法可以实现吗？会非常非常好。非常感谢这里的每一个帮助。

【参考方案2】：

你的匹配条件

abs(((referencelist - sample[i])/sample[i])*10^6)) < 0.5

可以改写为

sample[i] * (1 - eps) < referencelist < sample[i] * (1 + eps)

eps = 0.5E-6.

使用它，我们可以使用 non-equi-join 在referencelist 中为每个sample 查找所有 匹配项（不仅是最近的！）：

library(data.table)
options(digits = 10)
eps <- 0.5E-6 # tol * 1E6
setDT(referencelist)[.(value = sample, 
                       lower = sample * (1 - eps), 
                       upper = sample * (1 + eps)), 
                     on = .(ref > lower, ref < upper), .(name, value, reference = x.ref)]

再现了预期的结果：

   name     value reference
1:    A 154.00315 154.00312
2:    G 159.02991 159.02992
3:    B 154.07688 154.07685
4:    E 156.77312 156.77310

作为对OP's comment 的回应，假设我们将referencelist2 修改为F = 154.00320，那么这也会被捕获：

setDT(referencelist2)[.(value = sample, 
                       lower = sample * (1 - eps), 
                       upper = sample * (1 + eps)), 
                     on = .(ref > lower, ref < upper), .(name, value, reference = x.ref)]

   name     value reference
1:    A 154.00315 154.00312
2:    F 154.00315 154.00320
3:    G 159.02991 159.02992
4:    B 154.07688 154.07685
5:    E 156.77312 156.77310