R语言解决矩阵距离的问题

发表于 | 分类于 ，

刚刚群里有人提问，“请教个问题，有两个矩阵A和B，A矩阵n行2列，2列分别代表经度和纬度，B矩阵相同（行数不同，A有1万行，B有350万行），需要找出B矩阵中每一行数据，在A矩阵中距离最近的行号，该采用什么方法。”

乍一看，如果矩阵A和矩阵B是同型矩阵的话，那么可以直接对矩阵A和矩阵B做减法，然后再计算距离。但是，在A和B非同型矩阵的情况下。貌似只能用循环的方式解决了。当然，用lapply()等泛函应该也可以。我想到一种解法，不过十分繁琐，而且用到了最不推荐的循环方式。鉴于目前没有想到什么更好的办法，姑且把目前的想法记下来。

由于提问者给的矩阵维度过大，算完估计黄花菜都凉了。所以我把数据维度缩小几个数量级，以下是代码。

1
2
3
4
5
6
7
8
# 生成数据
set.seed(12)
A <- matrix(rnorm(1e3), ncol = 2)
B <- matrix(runif(1e2, 2, 6), ncol = 2)
dfA <- as.data.frame(A)
dfB <- as.data.frame(B)
head(dfA)
1
2
3
4
5
6
7
#>           V1          V2
#> 1 -1.4805676 -1.54930140
#> 2  1.5771695  0.87145355
#> 3 -0.9567445  0.04878753
#> 4 -0.9200052  0.17027381
#> 5 -1.9976421  0.07466159
#> 6 -0.2722960  0.28756053
1
head(dfB)
1
2
3
4
5
6
7
#>         V1       V2
#> 1 2.023439 4.884397
#> 2 5.519670 3.213451
#> 3 4.788179 5.259515
#> 4 4.394472 4.067191
#> 5 5.772587 5.793892
#> 6 2.257782 2.378813
1
2
3
4
5
6
7
8
9
10
# 计算矩阵A每一行到矩阵B每一行的距离
distance <- c()
nm <- c()
for (i in 1:nrow(dfA)) for (j in 1:nrow(dfB)) {
  distance <- c(sqrt((dfA[[1]][i] - dfB[[1]][j]) ^ 2 + 
                        (dfA[[2]][i] - dfB[[2]][j]) ^ 2), distance)
  nm <- c(nm, paste("A", i, "B", j, sep = "-"))
}
names(distance) <- nm
head(distance)
1
2
#>  A-1-B-1  A-1-B-2  A-1-B-3  A-1-B-4  A-1-B-5  A-1-B-6 
#> 5.092386 6.108986 6.713148 7.358928 3.383767 6.072581

上面的代码运行效率其实是很低的，因为每经过一轮迭代，R就要为distance和nm两个向量重新分配内存。这里为什么不用distance <- numeric(nrow(dfA) * nrow(dfB))这种循环呢？可以看到,distance的长度是矩阵A和矩阵B行数的积，但是在for循环中，索引向量只有i和j，i和j均小于distance的长度。这就导致了循环无法继续。

1
2
3
4
5
6
7
8
9
10
11
12
13
library(dplyr)
library(tidyr)
dist_df <- data.frame(nm = names(distance), dst = unname(distance))
dist_df %>%
  separate(nm, into = c("matA", "rowA", "matB", "rowB"),
           sep = "-") %>%
  group_by(rowA) %>%
  mutate(min_dist = min(dst),
         diff_dist = dst - min_dist) %>%
  # print(n = 150)
  filter(diff_dist == 0) %>%
  select(1:4) %>%
  print(n = 20)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#> # A tibble: 500 x 4
#> # Groups:   rowA [500]
#>     matA  rowA  matB  rowB
#>    <chr> <chr> <chr> <chr>
#>  1     A     1     B     5
#>  2     A     2     B     5
#>  3     A     3     B     5
#>  4     A     4     B     5
#>  5     A     5     B     5
#>  6     A     6     B     5
#>  7     A     7     B     5
#>  8     A     8     B     5
#>  9     A     9     B     5
#> 10     A    10     B     5
#> 11     A    11     B     5
#> 12     A    12     B     5
#> 13     A    13     B     5
#> 14     A    14     B     5
#> 15     A    15     B     5
#> 16     A    16     B     5
#> 17     A    17     B     5
#> 18     A    18     B     5
#> 19     A    19     B     5
#> 20     A    20     B     5
#> # ... with 480 more rows

解决上面的问题，一共写了接近30行的代码。可以说已经是很繁琐了。但现在想想也就这一个办法(⊙o⊙)…。如果以后再想到其他的更简单的做法，再重写一篇博客。