两个数字向量上的 All-to-all setdiff，具有接受匹配的数字阈值

Question

我想做的或多或少是以下两个线程中讨论的问题的组合：

Perform non-pairwise all-to-all comparisons between two unordered character vectors --- The opposite of intersect --- all-to-all setdiff
Merge data frames based on numeric rownames within a chosen threshold and keeping unmatched rows as well

我有两个数值向量：

b_1 <- c(543.4591, 489.36325, 12.03, 896.158, 1002.5698, 301.569)
b_2 <- c(22.12, 53, 12.02, 543.4891, 5666.31, 100.1, 896.131, 489.37)

我想将 b_1 中的所有元素与 b_2 中的所有元素进行比较，反之亦然。

如果 b_1 中的 element_i 是 NOT 等于 中的任何 个数 range element_j ± 0.045 in b_2 then element_i 必须报告。

同样，如果 b_2 中的 element_j 是 NOT 等于 any 中的数字 range element_i ± 0.045 in b_1 则必须报element_j.

因此，基于上面提供的向量的示例答案将是：

### based on threshold = 0.045
in_b1_not_in_b2 <- c(1002.5698, 301.569)
in_b2_not_in_b1 <- c(22.12, 53, 5666.31, 100.1)

是否有 R 函数可以做到这一点？

Answer 1

矢量化野兽：

D <- abs(outer(b_1, b_2, "-")) > 0.045

in_b1_not_in_b2 <- b_1[rowSums(D) == length(b_2)]
#[1] 1002.570  301.569

in_b2_not_in_b1 <- b_2[colSums(D) == length(b_1)]
#[1]   22.12   53.00 5666.31  100.10

小时后...

Henrik shared a question complaining the memory explosion when using outer for long vectors: 。但是，outer 的内存瓶颈可以通过阻塞轻松消除。

f <- function (b1, b2, threshold, chunk.size = 5000) {

  n1 <- length(b1)
  n2 <- length(b2)
  chunk.size <- min(chunk.size, n1, n2)

  RS <- numeric(n1)  ## rowSums, to be accumulated
  CS <- numeric(n2)  ## colSums, to be accumulated

  j <- 0
  while (j < n2) {
    chunk.size_j <- min(chunk.size, n2 - j)
    ind_j <- (j + 1):(j + chunk.size_j)
    b2_j <- b2[ind_j]
    i <- 0
    while (i < n1) {
      chunk.size_i <- min(chunk.size, n1 - i)
      ind_i <- (i + 1):(i + chunk.size_i)
      M <- abs(outer(b1[ind_i], b2_j, "-")) > threshold
      RS[ind_i] <- RS[ind_i] + rowSums(M)
      CS[ind_j] <- CS[ind_j] + colSums(M)
      i <- i + chunk.size_i
      }
    j <- j + chunk.size_j
    }

  list(in_b1_not_in_b2 = b1[RS == n2], in_b2_not_in_b1 = b2[CS == n1])
  }

有了这个函数，outer 使用的内存永远不会超过存储两个 chunk.size x chunk.size 矩阵。现在让我们做点疯狂的事吧。

b1 <- runif(1e+5, 0, 10000)
b2 <- b1 + runif(1e+5, -1, 1)

如果我们做一个简单的outer，我们需要内存来存储两个1e+5 x 1e+5矩阵，最多149GB。但是，在我的 Sandy Bridge (2011) 笔记本电脑上，只有 4 GB RAM，计算是可行的。

system.time(oo <- f(b1, b2, 0.045, 5000))
#   user  system elapsed 
#365.800 167.348 533.912

鉴于我们一直在使用非常糟糕的算法，性能实际上已经足够好了

这里的所有答案都是穷举搜索，复杂度 length(b1) x length(b2)。如果我们处理排序数组，我们可以将其减少到 length(b1) + length(b2)。但是这种深度优化的算法只能用编译型语言来实现才能获得效率。

Answer 2

这是另一种方法

in_b1_not_in_b2 <- b_1[sapply(b_1, function(x) !any(abs(x - b_2) <= 0.045))]
in_b1_not_in_b2
#[1] 1002.570  301.569

in_b2_not_in_b1 <- b_2[sapply(b_2, function(x) !any(abs(x - b_1) <= 0.045))]
in_b2_not_in_b1
#[1]   22.12   53.00 5666.31  100.10

Answer 3

如果您乐于使用非 base 包，data.table::inrange 是一个方便的功能。

x1[!inrange(x1, x2 - 0.045, x2 + 0.045)]
# [1] 1002.570  301.569

x2[!inrange(x2, x1 - 0.045, x1 + 0.045)]
# [1]   22.12   53.00 5666.31  100.10

inrange 在更大的数据集上也很有效。在例如1e5 向量，inrange 比其他两个替代方案快 > 700 倍：

n <- 1e5
b1 <- runif(n, 0, 10000)
b2 <- b1 + runif(n, -1, 1)

microbenchmark(
  f1 = f(b1, b2, 0.045, 5000),
  f2 = list(in_b1_not_in_b2 = b1[sapply(b1, function(x) !any(abs(x - b2) <= 0.045))],
       in_b2_not_in_b1 = b2[sapply(b2, function(x) !any(abs(x - b1) <= 0.045))]),
  f3 = list(in_b1_not_in_b2 = b1[!inrange(b1, b2 - 0.045, b2 + 0.045)],
       in_b2_not_in_b1 = b2[!inrange(b2, b1 - 0.045, b1 + 0.045)]),
  unit = "relative", times = 10)
# Unit: relative
#  expr      min       lq     mean   median        uq       max neval
#    f1 1976.931 1481.324 1269.393 1103.567 1173.3017 1060.2435    10
#    f2 1347.114 1027.682  858.908  766.773  754.7606  700.0702    10
#    f3    1.000    1.000    1.000    1.000    1.0000    1.0000    10

是的，他们给出了相同的结果：

n <- 100
b1 <- runif(n, 0, 10000)
b2 <- b1 + runif(n, -1, 1)

all.equal(f(b1, b2, 0.045, 5000),
          list(in_b1_not_in_b2 = b1[sapply(b1, function(x) !any(abs(x - b2) <= 0.045))],
               in_b2_not_in_b1 = b2[sapply(b2, function(x) !any(abs(x - b1) <= 0.045))]))
# TRUE

all.equal(f(b1, b2, 0.045, 5000),
          list(in_b1_not_in_b2 = b1[!inrange(b1, b2 - 0.045, b2 + 0.045)],
               in_b2_not_in_b1 = b2[!inrange(b2, b1 - 0.045, b1 + 0.045)]))
# TRUE

searching for inrange on SO 时的几个相关的、可能有用的答案。

两个数字向量上的 All-to-all setdiff，具有接受匹配的数字阈值

All-to-all setdiff on two numeric vectors with a numeric threshold for accepting matches

r

vector

compare

set-difference