R提高嵌套for（）循环的效率，以在大型数据集中进行简单的距离计算

Question

我有两组点（带有 x、y、z 坐标）数据 dtmT（113k 观测值）和 ptmT（200k 观测值）。对于 dtmT 中的每个点，我都希望计算到 ptmT 中某个点的最短距离。我对 R 非常陌生，没有其他编程背景，所以我嵌套了 for 循环，因此对于 dtmT 的每个点，它计算该点与 ptmT 中每个点之间的距离并将其存储在矩阵中（ptmTDistM） .我使用后循环将矩阵中每一列的最小值作为向量获取，然后使用 cbind 将其附加回 dtmT，以便最终产品是 dtmT，其中 x,y,z,Dist 表示可能的最短距离dtmT 指向 ptmT 内的任何点。这对于 5 次观察和最多 500 次非常有效，但是当我尝试使用 5,000 次并且完整的数据集是 dtmT 中的 113K 观察和 ptmT 中的 200k 时，它会挂起。我最初使用数据框对此进行了编程，但已经阅读了一些让我尝试使用矩阵的问题和答案。我也明白使用向量和 lapply 组是最好的，我不确定如何将嵌套的 for 循环转换到 lapply 组，特别是因为索引对于我现在如何获得它非常重要。我也看过 Dist() 但不确定如何应用它来获得我需要的东西。

提供了每个数据集的前 5 个观察结果以及我到目前为止所做的工作。

非常感谢您的帮助！

#first 5 observations of ptmT dataset
ptmT <- c(621019.2, 621024.2, 621023.7, 621018.3, 621019.2, 2701229.1, 
2701231.2, 2701231.9, 2701230.2, 2701229.1, 2071.5, 2080.0, 2080.0, 2071.5, 
2071.5)
dim(ptmT) <- c(5,3)
colnames(ptmT) <- c("XP", "YP", "ZP")

#first 5 observations of dtmT dataset
dtmT <- c( 621757.360, 621757.360, 621757.419, 621757.536, 
621757.540,2701071.810, 2701071.810, 2701071.814, 2701071.843, 2701071.844, 
2089.210, 2088.110, 2070.435, 2053.536, 2052.951)
dim(dtmT) <- c(5,3)
colnames(dtmT) <- c("X", "Y", "Z")

dtmTDist <- 0
ptmTDist <- 0
ptmTDistM <- matrix(data = NA, nrow = length(ptmT[,1]), ncol = 
length(dtmT[,1]))

require (svMisc)
for (row in 1:nrow(dtmT))   {
    progress(row)               
    X <- dtmT[row, "X"]
    Y <- dtmT[row, "Y"]
    Z <- dtmT[row, "Z"]

        for (i in 1:nrow(ptmT)) {
            X2 <- ptmT[i, "XP"]
            Y2 <- ptmT[i, "YP"]
            Z2 <- ptmT[i, "ZP"]

            D <- sqrt((X - X2)^2 + (Y - Y2)^2 + (Z - Z2)^2)
            ptmTDistM[i,row] <- D
            }
    }
Dist <- apply(ptmTDistM, 2, min)
dtmT2 <- cbind(dtmT,Dist)

Answer 1

您可以使用最近邻搜索包，例如 https://github.com/jefferis/RANN，它将为每个查询点返回最近点及其与参考点的距离（使用有效的空间索引）

P <- 200000
ptmT <- data.frame(x=runif(P),y=runif(P),z=runif(P))
D <- 113000
dtmT <- data.frame(x=runif(D),y=runif(D),z=runif(D))
library(RANN)
res <- nn2(ptmT,dtmT,1)

Answer 2

通过利用 R 的矢量算法等功能，您可能能够获得一些性能改进。但是，任何需要检查集合 A 中的每个点与集合 B 中的每个点的方法都将变得非常苛刻，因为这两个集合都变大了，因为要进行的比较的数量与 O(m*n) 成比例，其中 m 和 n 是大小两组。

有时有助于解决此问题的一个技巧是按地理位置对您的集合进行分块，并使用该分块来确定您实际测试的对。

例如，在 2D 中：

• 从 A 中随机挑选 100 个点。对于每个点，通过与 B 中的每个点进行比较，找出与 B 中最近邻居的距离。（总计：100*n 次比较。）
• 设 h = 上面的最大值。
• 将您的空间分成块，大小为 2h x 2h。对于 A 中的任何点，您可以几乎确定它在 B 中的最近邻居将位于其自己的块中，或者位于 8 个相邻块之一中。
• 对于 B 中的每个点，确定它位于哪个块中，并设置一个索引或向量列表，以便您可以轻松引用“B 中位于块 [x,y] 中的所有点”。
• 对于 A 中的每个点 P，找出它所在的块，并注意它与该块的最近边界的距离（称为 d），然后针对 B 中位于同一块中的所有点对其进行测试. （这些是您可以利用矢量算术的地方。）
• 如果你在 B 中发现一个点接近或等于 d，那么这肯定是最近的邻居，你可以停下来。
• 否则，如果您找到的最近点比 d 更远，或者您的搜索区域中根本没有来自 B 的点，则将搜索扩展到相邻块，并设置 d
• 重复直到找到最近的点，然后转到下一个 P 直到完成。

这意味着对于 A 中的每个点，您只针对 B 中的少量附近点进行测试，而不是测试地图上的所有内容。尽管搜索方法更复杂，但对于较大的 m & n，您应该会看到更好的搜索时间。

如果您的数据点分布非常不均匀，您可能需要使用网格形状；理想情况下，“块”被设计成每个只包含 B 的几个成员。

另外，一个次要的经济：请注意，最小化距离平方也会最小化距离。因此，您可以执行 sqrt(min(dist^2)) 而不是找到 min(dist)，这将为您节省大量平方根运算，这是值得的。

Answer 3

由于我们无法避免计算两点之间的距离，（除非之前已计算出完全相同的点对）您肯定需要进行 113,000*200,000 次计算。

加快速度的唯一方法是尝试使计算尽可能并行。

您绝对应该尝试 cmets 中建议的并行包。

这是我在 R 中使用 apply 函数的解决方案，它尝试尽可能多地进行矢量化和计算。

#Function to calculate Euclidean distance. We can simply use matrix algebra here.
computeDistance <- function(P,Q){
  D <- sqrt(sum((P-Q)^2))
  return(D)
}

#We use one apply row-wise on dtmT and for compute distance with each row in ptmT.
#Since this is a perfectlly parallel process, apply will be substantially faster than a for loop
distMat <- apply(dtmT, MARGIN = 1, function(p){apply(ptmT,MARGIN = 1,FUN = function(q){computeDistance(p,q)})})

#Calculate minimum of each column to get the minimum distance
minDist <- apply(distMat,2,min)

#Attach to dtmT
dtmTFinal <- cbind(dtmT,"Minimum_Distance" = minDist)

我在 5000*5000 的情况下尝试了这个，在普通笔记本电脑上花了大约一分钟。

希望这会有所帮助。

Answer 4

这里的一个主要问题是内存，因为您的 113k x 200k 矩阵将占用大约 170 GB 内存。但是，您永远不需要完整的矩阵。相反，您只需要每行的最小值。此外，您可以以矢量化方式计算此最小值，只留下一个循环：

Dist <- vector(length = nrow(dtmT), mode = "numeric")
for (row in 1:nrow(dtmT))   {
  X <- dtmT[row, "X"]
  Y <- dtmT[row, "Y"]
  Z <- dtmT[row, "Z"]

  Dist[row] <- sqrt(min((X - ptmT[ ,"XP"])^2 + (Y - ptmT[ ,"YP"])^2 + (Z - ptmT[ , "ZP"])^2))
}
cbind(dtmT,Dist)

现在这个循环是“令人尴尬的并行”，您可以使用 foreach 进行并行化：

library(foreach)
library(doParallel)
registerDoParallel(cores = 4)
Dist <- foreach (row = 1:nrow(dtmT), .combine = c) %dopar% {
  X <- dtmT[row, "X"]
  Y <- dtmT[row, "Y"]
  Z <- dtmT[row, "Z"]

  sqrt(min((X - ptmT[ ,"XP"])^2 + (Y - ptmT[ ,"YP"])^2 + (Z - ptmT[ , "ZP"])^2))  
}
cbind(dtmT,Dist)

将应用使用for 循环的替代方法。将其与更紧凑的符号相结合，我们得到：

apply(dtmT, 1, function(x) sqrt(min(colSums((x-t(ptmT))^2))))

同样，apply 可以轻松并行化。将其应用于尺寸小 10 倍的问题在双核机器上给出：

library(parallel)
cl <- makeForkCluster(2)
dtmT <- matrix(runif(3 * 11300), ncol = 3)
ptmT <- matrix(runif(3 * 200000), ncol = 3)
system.time(Dist <- parApply(cl, dtmT, 1, function(x) sqrt(min(colSums((x-t(ptmT))^2)))))
#>        User      System verstrichen 
#>       0.021       0.004      34.474
head(cbind(dtmT, Dist))
#>                                            res
#> [1,] 0.9111543 0.5971182 0.8725145 0.010714792
#> [2,] 0.4893960 0.3321890 0.7440035 0.008545801
#> [3,] 0.3637524 0.6051168 0.7955850 0.003792442
#> [4,] 0.6684364 0.1819622 0.2487011 0.017937629
#> [5,] 0.6761877 0.1731773 0.3214378 0.011912805
#> [6,] 0.8060648 0.7789117 0.1673685 0.012680877