检查矩阵中是否存在多个向量，需要速度答案

【问题标题】：Checking whether a number of vectors exist in a matrix, need speed检查矩阵中是否存在多个向量，需要速度
【发布时间】：2021-03-25 04:54:55
【问题描述】：

我有两个大型数值矩阵，需要检查其中一个中的哪些行存在于另一个中（存在等于相等）。这是我的代码：

myMatrix1 <- rbind(c(1,2,3),c(4,5,6),c(7,8,9))
myMatrix2 <- rbind(c(10,11,12),c(4,5,6),c(13,14,15))

logicalMatrix <- apply(myMatrix1,1,checkForEquality)
result <- apply(logicalMatrix,1,any)

checkForEquality <- function(x){
  apply(myMatrix2, 1, innerFcn, oneRow = x)
}
innerFcn <- function(x, oneRow){
  isTRUE(all.equal(x, oneRow))
}

结果是

[1] FALSE  TRUE FALSE

使用两个 2067*198 矩阵，这在我的机器上需要 350 秒。通过 CPU 并行化，我想我可以将它降低到 15 秒左右。不幸的是，任何超过 1 秒的时间都是不可接受的。我需要一些方向。如果重要的话，矩阵只包含 0、1 和 2。

【问题讨论】：

在我倾向于用来确定“相同性”的所有基本方法中，all.equal 是最慢的。考虑将all.equal(x,y) 替换为all(x==y)。特别是因为您需要返回一个简单的真/假，一旦all.equal 发现它们不相等，仍然会在一定程度上进行比较以报告差异。看来您真的想要一个短路操作，而不是“向我报告所有差异”。

标签： r performance matrix

【解决方案1】：

您可以使用欧几里得距离，可以表示为矢量化/线性代数运算：

dist_xy <- outer(rowSums(x^2), rowSums(y^2), '+') - tcrossprod(x, 2 * y))

基准测试：

nr = 1e5
nc = 200
x = t(sample(0:2, size = nc, replace = TRUE))
y = matrix(sample(0L:2L, size = nc * nr, replace = TRUE), nrow = nr)

all.equal(apply(y, 1, function(z) identical(z, x)),
          drop(outer(rowSums(x^2), rowSums(y^2), '+') == tcrossprod(x, 2 * y)))
# [1] TRUE

microbenchmark::microbenchmark(
  A = apply(y, 1, function(z) identical(z, x)),
  B = apply(y, 1, function(z) all(z == x)),
  C = drop(outer(rowSums(x^2), rowSums(y^2), '+') == tcrossprod(x, 2 * y)),
  times = 3
)
Unit: milliseconds
 expr      min       lq     mean   median       uq      max neval
    A 543.1362 559.9647 585.3627 576.7931 606.4760 636.1589     3
    B 609.7400 636.5922 667.3405 663.4445 696.1408 728.8370     3
    C 368.2808 416.4194 441.9118 464.5580 478.7273 492.8965     3

如果您有更大的矩阵，这应该会受益更多。

【讨论】：

我真的很惊讶这比 identical 或 all( == ) 方法更快 - 你有直觉为什么会这样吗？
我已经尝试过您的示例 (C)，实际上它比其他示例要快得多，因为您可以一次将其应用于整个矩阵。使用一个 2067*198 矩阵和另一个 6201*198 矩阵，它比相同（）快约 15 倍。您在示例中仅使用了一个向量 (x)。无论如何，它仍然不够快，我将尝试用 GPU 并行化做同样的事情。
是的，线性代数很快，尤其是当您将 R 与一些快速矩阵库（例如 OpenBLAS 或 MKL）链接时。

【解决方案2】：

计算距离矩阵：

X <- rbind(c(1,2,3),c(4,5,6),c(7,8,9))
Y <- rbind(c(10,11,12),c(4,5,6),c(13,14,15))

library(pracma)
ind <- which(distmat(X, Y) == 0L, arr.ind = TRUE)
#     row col
#[1,]   2   2

X[ind[, 1],]
#[1] 4 5 6

Y[ind[, 2],]
#[1] 4 5 6

如果你需要考虑浮点精度，使用这个：

ind <- which(abs(distmat(X, Y)) < tol, arr.ind = TRUE)

【讨论】：

【解决方案3】：

从all.equal 切换后，您可以获得大约 10 倍的加速。 all(x == y) 和 pracma::distmat 和 identical(x, y) 都快 10 倍左右，其中 identical(x, y) 最快（注意：这是用于长度为 200 的比较，如您的数据中一样。对于更长的向量，all(x == y) 可能是更快！）。

# sample data
## demo on single row vs matrix comparison
nr = 1e5
nc = 200
x = sample(0:2, size = nc, replace = TRUE)
y = matrix(sample(0L:2L, size = nc * nr, replace = TRUE), nrow = nr)


library(pracma)
microbenchmark::microbenchmark(
  apply(y, 1, function(z) identical(z, x)),
  apply(y, 1, function(z) all(z == x)),
  apply(y, 1, function(z) all.equal(z, x)),
  which(distmat(x, y) == 0L, arr.ind = TRUE),
  times = 2
)
# Unit: milliseconds
#                                        expr       min        lq      mean    median        uq       max
#    apply(y, 1, function(z) identical(z, x))  619.1912  619.1912  644.3034  644.3034  669.4156  669.4156
#        apply(y, 1, function(z) all(z == x))  759.4982  759.4982  789.4765  789.4765  819.4548  819.4548
#    apply(y, 1, function(z) all.equal(z, x)) 7618.6853 7618.6853 7657.7665 7657.7665 7696.8477 7696.8477
#  which(distmat(x, y) == 0L, arr.ind = TRUE)  824.8337  824.8337  899.2349  899.2349  973.6360  973.6360

您显然可以使用 Rcpp 更快，但即使在 for 循环中进行逐行比较在 R 中也会非常快（感谢 JIT 编译）。

我要尝试优化的下一个点是提前停止 - 如果您找到匹配项，则无需再检查与该行的任何比较，因此您可以继续下一行。使用apply 无法很好地做到这一点，但使用for 循环可以更好地控制。

【讨论】：

使用更大的向量，all(x==y) 在性能上超过了identical。
@r2evans 很有趣，很高兴知道——我想知道为什么会这样。我会在答案中添加这个警告。
请注意，比较算法只能应用于 rowSums 相等的行簇。对矩阵中的项进行平方，然后取 rowSums 将进一步减少公共行数。

【解决方案4】：

从a中减去b，然后求绝对值的rowSum，以防正负差抵消。零和行将是相同的：

diff <- a - b
idrows <- rowSum(abs(diff))

【讨论】：

您假设这两个矩阵中的同一行是相等的。这不是我阅读问题的方式。

【解决方案5】：

试试吧：

asplit(m2,1) %in% asplit(m1,1)

一些中等矩阵的例子：

dims<-c(1920,198)
set.seed(4)
m1<-matrix(sample(0:2,prod(dims),TRUE),ncol=dims[2])
m2<-matrix(sample(0:2,prod(dims),TRUE),ncol=dims[2])
#simulate some rows in m2 which are equal to some rows in m1
samind<-sample(nrow(m1),50)
samind2<-sample(nrow(m2),50)
m2[samind2,]<-m1[samind,]
system.time(res<-asplit(m2,1) %in% asplit(m1,1))
#   user  system elapsed 
#  0.407   0.000   0.406
#Check whether the result is correct
identical(which(res),sort(samind2))
#[1] TRUE

【讨论】：