在R中高效地计算一个点与一组点之间的所有距离

Question

首先，我刚开始学习R（昨天开始）。

我有两组点，data和centers，第一组大小为n，第二组大小为K（例如，n=3823，K=10），对于第一组中的每个i，我需要找到第二组中距离最小的j。

我的想法很简单：对于每个i，让dist[j]成为i和j之间的距离，我只需要使用which.min(dist)来找到我要找的内容。

每个点都是由64个double数组成的。

> dim(data)
[1] 3823   64
> dim(centers)
[1] 10 64

我已经尝试过了

for (i in 1:n) {
  for (j in 1:K) {
    d[j] <- sqrt(sum((centers[j,] - data[i,])^2))
  }
  S[i] <- which.min(d)
}

这段代码非常缓慢（当n = 200时，需要超过40秒!!）。我写的最快解决方案是：

distance <- function(point, group) {
  return(dist(t(array(c(point, t(group)), dim=c(ncol(group), 1+nrow(group)))))[1:nrow(group)])
}

for (i in 1:n) {
  d <- distance(data[i,], centers)
  which.min(d)
}

即使它进行了很多我不使用的计算（因为dist(m)计算了m的所有行之间的距离），它比另一个更快（有人能解释为什么吗？），但对于我需要的速度来说还不够快，因为它不仅仅会被使用一次。此外，distance代码非常丑陋。我试图用以下内容替换它：

distance <- function(point, group) {
  return (dist(rbind(point,group))[1:nrow(group)])
}

但这似乎慢了两倍。我还尝试使用dist来计算每对数据，但它也很慢。

我不知道该怎么办了。看起来我做错了什么。你有没有更高效的方法？

附：我需要手动实现k-means（这是我的任务的一部分）。我相信只需要欧几里得距离，但我还不确定，所以我希望有一些代码，可以轻松地替换距离计算。 stats::kmeans可以在不到一秒钟内完成所有计算。

- dbarbosa

---

这里的人们有点不喜欢做作业...所以尽量专注于一个具体的问题。 - aL3xa

Answer 1

与其遍历数据点，你可以将其简化为矩阵操作，这意味着你只需要遍历 K 即可。

# Generate some fake data.
n <- 3823
K <- 10
d <- 64
x <- matrix(rnorm(n * d), ncol = n)
centers <- matrix(rnorm(K * d), ncol = K)

system.time(
  dists <- apply(centers, 2, function(center) {
    colSums((x - center)^2)
})
)

运行于：

utilisateur     système      écoulé 
      0.100       0.008       0.108 

在我的笔记本电脑上。

- Jonathan Chang

---

+1 比我的计算 dists 矩阵的方法更好。这是一个很棒的技巧，通过自动复制向量添加或从矩阵中减去。 - Marek

我正在尝试使用您的解决方案，但是您的矩阵已经转置了。是否有一种像您对列所做的那样减去行的方法？ - dbarbosa

我尝试使用apply对线条进行减法，但速度不如您的解决方案。我现在正在转置矩阵并使用您的代码，速度非常快！非常感谢!!! 还要感谢您提供完整答案和小例子，以及使用system.time。Merci beaucoup :) - dbarbosa

Answer 2

rdist()是{fields}包中的一个R函数，能够快速计算矩阵格式下两组点之间的距离。

https://www.image.ucar.edu/~nychka/Fields/Help/rdist.html

使用方法：

library(fields)
#generating fake data
n <- 5
m <- 10
d <- 3

x <- matrix(rnorm(n * d), ncol = d)
y <- matrix(rnorm(m * d), ncol = d)

rdist(x, y)
          [,1]     [,2]      [,3]     [,4]     [,5]
 [1,] 1.512383 3.053084 3.1420322 4.942360 3.345619
 [2,] 3.531150 4.593120 1.9895867 4.212358 2.868283
 [3,] 1.925701 2.217248 2.4232672 4.529040 2.243467
 [4,] 2.751179 2.260113 2.2469334 3.674180 1.701388
 [5,] 3.303224 3.888610 0.5091929 4.563767 1.661411
 [6,] 3.188290 3.304657 3.6668867 3.599771 3.453358
 [7,] 2.891969 2.823296 1.6926825 4.845681 1.544732
 [8,] 2.987394 1.553104 2.8849988 4.683407 2.000689
 [9,] 3.199353 2.822421 1.5221291 4.414465 1.078257
[10,] 2.492993 2.994359 3.3573190 6.498129 3.337441

- Deuterium

Answer 3

我的解决方案：

# data is a matrix where each row is a point
# point is a vector of values
euc.dist <- function(data, point) {
  apply(data, 1, function (row) sqrt(sum((point - row) ^ 2)))
}

您可以尝试一下，例如：

x <- matrix(rnorm(25), ncol=5)
euc.dist(x, x[1,])

- Adriano Rivolli

Answer 4

你可能想要查看一下apply函数。

例如，这段代码

for (j in 1:K)
    {
    d[j] <- sqrt(sum((centers[j,] - data[i,])^2))
    }

可以轻松地被类似的东西替换

dt <- data[i,]
d <- apply(centers, 1, function(x){ sqrt(sum(x-dt)^2)})

你肯定可以进一步优化它，但我希望你能理解重点。

- nico

---

谢谢...这段代码比我写的第一段要快，但还远远不及使用“distance”的奇怪代码。 - dbarbosa

@dbarbosa：嗯，显然 stats::kmeans 包使用了编译后的代码，因此速度更快。只需键入 kmeans，即可查看其源代码。 :) - nico

Answer 5

dist 之所以运行速度快，是因为它向量化并调用内部 C 函数。
您的循环代码可以通过多种方式进行向量化。

例如，要计算 data 和 centers 之间的距离，您可以使用 outer：

diff_ij <- function(i,j) sqrt(rowSums((data[i,]-centers[j,])^2))
X <- outer(seq_len(n), seq_len(K), diff_ij)

这将给你一个距离的 n x K 矩阵。并且应该比循环要快得多。

然后，您可以使用 max.col 在每行中找到最大值（请参见帮助文档，在有多个最大值时有一些细微差别）。X 必须被否定，因为我们正在寻找最小值。

CL <- max.col(-X)

为了在R中提高效率，您应该尽可能地使用向量化。在许多情况下，循环可以被向量化的替代方法所取代。请查看帮助文档，了解rowSums（也描述了rowMeans、colSums、rowSums）、pmax、cumsum等内容。您可以搜索SO，例如https://stackoverflow.com/search?q=[r]+avoid+loop（复制并粘贴此链接，我不知道如何使其可点击）以获取一些示例。

- Marek

---

你好，我正在尝试使用你的代码，但它无法正常工作。我尝试使用与@Jonathan Chang编写的相同代码，并添加了：system.time(outer(seq_len(n), seq_len(K), function(i,j) sqrt(rowSums((x[,i]-centers[,j])^2))))，但是我遇到了以下错误：Error in dim(robj) <- c(dX, dY) : dims [product 38230] do not match the length of object [64]你看出问题在哪里了吗？ - dbarbosa

其实我一开始并没有理解 outer（我以为它会为每对参数调用一次函数）。现在我明白了，谢谢你，这很有用！还有，感谢你告诉我关于 max.col。 - dbarbosa