（速度挑战）根据欧几里得距离计算两个矩阵的行之间的距离矩阵有什么更快的方法？

Question

首先，这不是calculating Euclidean distance between two matrices的问题。

假设我有两个矩阵x 和y，例如，

set.seed(1)
x <- matrix(rnorm(15), ncol=5)
y <- matrix(rnorm(20), ncol=5)

在哪里

> x
           [,1]       [,2]      [,3]       [,4]       [,5]
[1,] -0.6264538  1.5952808 0.4874291 -0.3053884 -0.6212406
[2,]  0.1836433  0.3295078 0.7383247  1.5117812 -2.2146999
[3,] -0.8356286 -0.8204684 0.5757814  0.3898432  1.1249309

> y
            [,1]       [,2]        [,3]       [,4]        [,5]
[1,] -0.04493361 0.59390132 -1.98935170 -1.4707524 -0.10278773
[2,] -0.01619026 0.91897737  0.61982575 -0.4781501  0.38767161
[3,]  0.94383621 0.78213630 -0.05612874  0.4179416 -0.05380504
[4,]  0.82122120 0.07456498 -0.15579551  1.3586796 -1.37705956

然后我想得到维度为 3×4 的距离矩阵 distmat，其中元素 distmat[i,j] 是来自 norm(x[1,]-y[2,],"2") 或 dist(rbind(x[1,],y[2,])) 的值。

• 我的代码

distmat <- as.matrix(unname(unstack(within(idx<-expand.grid(seq(nrow(x)),seq(nrow(y))), d <-sqrt(rowSums((x[Var1,]-y[Var2,])**2))), d~Var2)))

给了

> distmat
         [,1]     [,2]     [,3]     [,4]
[1,] 3.016991 1.376622 2.065831 2.857002
[2,] 4.573625 3.336707 2.698124 1.412811
[3,] 3.764925 2.235186 2.743056 3.358577

但我认为我的代码在使用大量行的x 和y 时不够优雅或高效。

• 目标

我期待更快和更优雅的代码 base R 来实现这个目标。提前感谢！

• 基准模板（更新中）

为方便起见，您可以使用以下基准测试您的代码是否更快：

set.seed(1)
x <- matrix(rnorm(15000), ncol=5)
y <- matrix(rnorm(20000), ncol=5)
# my customized approach
method_ThomasIsCoding_v1 <- function() {
  as.matrix(unname(unstack(within(idx<-expand.grid(seq(nrow(x)),seq(nrow(y))), d <-sqrt(rowSums((x[Var1,]-y[Var2,])**2))), d~Var2)))
}
method_ThomasIsCoding_v2 <- function() {
  `dim<-`(with(idx<-expand.grid(seq(nrow(x)),seq(nrow(y))), sqrt(rowSums((x[Var1,]-y[Var2,])**2))),c(nrow(x),nrow(y)))
}
method_ThomasIsCoding_v3 <- function() {
  `dim<-`(with(idx1<-list(Var1 = rep(1:nrow(x), nrow(y)), Var2 = rep(1:nrow(y), each = nrow(x))), sqrt(rowSums((x[Var1,]-y[Var2,])**2))),c(nrow(x),nrow(y)))
}
# approach by AllanCameron
method_AllanCameron <- function()
{
  `dim<-`(sqrt(rowSums((x[rep(1:nrow(x), nrow(y)),] - y[rep(1:nrow(y), each = nrow(x)),])^2)), c(nrow(x), nrow(y)))
}
# approach by F.Prive
method_F.Prive <- function() {
  sqrt(outer(rowSums(x^2), rowSums(y^2), '+') - tcrossprod(x, 2 * y))
}
# an existing approach by A. Webb from https://stackoverflow.com/a/35107198/12158757
method_A.Webb <- function() {
  euclidean_distance <- function(p,q) sqrt(sum((p - q)**2))
  outer(
    data.frame(t(x)),
    data.frame(t(y)),
    Vectorize(euclidean_distance)
  )
}



bm <- microbenchmark::microbenchmark(
  method_ThomasIsCoding_v1(),
  method_ThomasIsCoding_v2(),
  method_ThomasIsCoding_v3(),
  method_AllanCameron(),
  method_F.Prive(),
  # method_A.Webb(),
  unit = "relative",
  check = "equivalent",
  times = 10
)
bm

这样

Unit: relative
                       expr      min       lq     mean   median       uq      max neval
 method_ThomasIsCoding_v1() 9.471806 8.838704 7.308433 7.567879 6.989114 5.429136    10
 method_ThomasIsCoding_v2() 4.623405 4.469646 3.817199 4.024436 3.703473 2.854471    10
 method_ThomasIsCoding_v3() 4.881620 4.832024 4.070866 4.134011 3.924366 3.367746    10
      method_AllanCameron() 5.654533 5.279920 4.436071 4.772527 4.184927 3.157814    10
           method_F.Prive() 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000    10

Answer 1

method_XXX <- function() {
  sqrt(outer(rowSums(x^2), rowSums(y^2), '+') - tcrossprod(x, 2 * y))
}

Unit: relative
                       expr       min        lq     mean    median        uq      max
 method_ThomasIsCoding_v1() 12.151624 10.486417 9.213107 10.162740 10.235274 5.278517
 method_ThomasIsCoding_v2()  6.923647  6.055417 5.549395  6.161603  6.140484 3.438976
 method_ThomasIsCoding_v3()  7.133525  6.218283 5.709549  6.438797  6.382204 3.383227
      method_AllanCameron()  7.093680  6.071482 5.776172  6.447973  6.497385 3.608604
               method_XXX()  1.000000  1.000000 1.000000  1.000000  1.000000 1.000000

Answer 2

我在管理一个可提供 3 倍加速的单线时保持它的简单和基本 R。

`dim<-`(sqrt(rowSums((x[rep(1:nrow(x), nrow(y)),] - y[rep(1:nrow(y), each = nrow(x)),])^2)), c(nrow(x), nrow(y)))
#          [,1]     [,2]     [,3]     [,4]
# [1,] 3.016991 1.376622 2.065831 2.857002
# [2,] 4.573625 3.336707 2.698124 1.412811
# [3,] 3.764925 2.235186 2.743056 3.358577

不幸的是，我的 PC 内存在使用大矩阵的测试时阻塞，所以我不得不将尺寸减少一个数量级来运行测试。

显示完整代码：

set.seed(1)
x <- matrix(rnorm(1500), ncol=5)
y <- matrix(rnorm(2000), ncol=5)
# my customized approach
method_ThomasIsCoding <- function() {
  as.matrix(unname(unstack(within(idx<-expand.grid(seq(nrow(x)),seq(nrow(y))), d <-sqrt(rowSums((x[Var1,]-y[Var2,])**2))), d~Var2)))
}
# an existing approach by A. Webb from https://stackoverflow.com/a/35107198/12158757
method_A.Webb <- function() {
  euclidean_distance <- function(p,q) sqrt(sum((p - q)**2))
  outer(
    data.frame(t(x)),
    data.frame(t(y)),
    Vectorize(euclidean_distance)
  )
}
# your approach
method_AllanCameron <- function()
{
  `dim<-`(sqrt(rowSums((x[rep(1:nrow(x), nrow(y)),] - y[rep(1:nrow(y), each = nrow(x)),])^2)), c(nrow(x), nrow(y)))
}

microbenchmark::microbenchmark(
  method_ThomasIsCoding(),
  method_A.Webb(),
  method_AllanCameron(),
  times = 10
)

结果

# Unit: milliseconds
#                     expr       min        lq      mean    median        uq       max neval
#  method_ThomasIsCoding()  63.08587  64.70988  69.59648  67.73379  75.90281  76.92903    10
#          method_A.Webb() 330.44824 349.90977 376.36962 368.52164 392.11780 446.57269    10
#    method_AllanCameron()  16.29938  18.20057  21.02634  20.45267  22.41767  31.28646    10

Answer 3

proxy 包有一个功能。

library(proxy)
dist(x, y)

     [,1]     [,2]     [,3]     [,4]    
[1,] 3.016991 1.376622 2.065831 2.857002
[2,] 4.573625 3.336707 2.698124 1.412811
[3,] 3.764925 2.235186 2.743056 3.358577