在R中计算稀疏的成对距离矩阵

Question

我有一个大小为NxM的矩阵，想要计算M个点之间的欧几里得距离，生成一个大小为NxN的矩阵。在我的问题中，N大约为100,000。由于我打算将此矩阵用于k最近邻算法，因此我只需要保留k个最小距离，所以生成的NxN矩阵非常稀疏。这与例如dist()生成的结果形成对比，后者会导致密集矩阵（并且可能对我的大小N造成存储问题）。

到目前为止，我发现的kNN包（例如knnflex、kknn等）似乎都使用密集矩阵。此外，Matrix包没有提供成对距离函数。

更接近我的目标的是，我看到spam包有一个nearest.dist()函数，可以只考虑小于某个阈值delta的距离。然而，在我的情况下，特定值的delta可能会产生太多的距离（使我必须密集地存储NxN矩阵），或者距离太少（使我无法使用kNN）。

我看到以前有关尝试使用bigmemory/biganalytics包执行k-means聚类的讨论，但在这种情况下似乎无法利用这些方法。

有人知道如何在R中以稀疏方式计算距离矩阵的函数/实现吗？我的备选方案是使用两个for循环并将结果保存在Matrix对象中。

- Christopher DuBois

---

只是确认一下...你知道 dist http://stat.ethz.ch/R-manual/R-patched/library/stats/html/dist.html 这个吗？ - Benjamin

抱歉，我没有清楚地解释为什么dist()对我的情况不够适用。它会生成一个密集矩阵，而且存储N x N矩阵有点麻烦。 - Christopher DuBois

你应该考虑接受这里的一个回答，如果你认为它实际上回答了问题（如果你认为自己的回答最好，那就接受自己的回答），或者编辑你的问题以澄清为什么它们不能回答。 - Tommy

“有点烦人”是轻描淡写了——如果N为100,000，那就是一个480Gb的矩阵。 - MichaelChirico

Answer 1

好的，我们不能让您诉诸使用for循环，不是吗 :)

当然，问题在于如何表示稀疏矩阵。一种简单的方法是仅将其包含最接近点的索引（并根据需要重新计算）。但在下面的解决方案中，我将距离（“d1”等）和索引（“i1”等）放入单个矩阵中：

sparseDist <- function(m, k) {
    m <- t(m)
    n <- ncol(m)
    d <- vapply( seq_len(n-1L), function(i) { 
        d<-colSums((m[, seq(i+1L, n), drop=FALSE]-m[,i])^2)
        o<-sort.list(d, na.last=NA, method='quick')[seq_len(k)]
        c(sqrt(d[o]), o+i) 
        }, numeric(2*k)
    )
    dimnames(d) <- list(c(paste('d', seq_len(k), sep=''),
        paste('i', seq_len(k), sep='')), colnames(m)[-n])
    d
}

在9个二维点上尝试它：

> m <- matrix(c(0,0, 1.1,0, 2,0, 0,1.2, 1.1,1.2, 2,1.2, 0,2, 1.1,2, 2,2),
              9, byrow=TRUE, dimnames=list(letters[1:9], letters[24:25]))
> print(dist(m), digits=2)
    a   b   c   d   e   f   g   h
b 1.1                            
c 2.0 0.9                        
d 1.2 1.6 2.3                    
e 1.6 1.2 1.5 1.1                
f 2.3 1.5 1.2 2.0 0.9            
g 2.0 2.3 2.8 0.8 1.4 2.2        
h 2.3 2.0 2.2 1.4 0.8 1.2 1.1    
i 2.8 2.2 2.0 2.2 1.2 0.8 2.0 0.9
> print(sparseDist(m, 3), digits=2)
     a   b   c   d   e   f   g   h
d1 1.1 0.9 1.2 0.8 0.8 0.8 1.1 0.9
d2 1.2 1.2 1.5 1.1 0.9 1.2 2.0  NA
d3 1.6 1.5 2.0 1.4 1.2 2.2  NA  NA
i1 2.0 3.0 6.0 7.0 8.0 9.0 8.0 9.0
i2 4.0 5.0 5.0 5.0 6.0 8.0 9.0  NA
i3 5.0 6.0 9.0 8.0 9.0 7.0  NA  NA

尝试在一个更大的问题上使用它（10k个点）。但是，在100k个点和更多维度上，它将需要很长时间（大约15-30分钟）。

n<-1e4; m<-3; m=matrix(runif(n*m), n)
system.time( d <- sparseDist(m, 3) ) # 9 seconds on my machine...

顺便提一下，在我写这篇文章的时候，我注意到你已经发布了答案：这里的解决方案大约快了两倍，因为它不会计算相同的距离两次（点1和点13之间的距离与点13和点1之间的距离相同）。

- Tommy

---

谢谢您的回答。我同意它大约快了两倍。然而，对于我的应用程序（kNN），我认为仅有距离矩阵的上三角实际上略微不方便。我想我可能会坚持使用我提交的代码的并行版本。再次感谢！ - Christopher DuBois

Answer 2

目前我正在使用以下方法，受到此答案的启发。 输出是一个n x k矩阵，其中元素(i,k)是距离第i个数据点第k近的数据点的索引。

n <- 10
d <- 3
x <- matrix(rnorm(n * d), ncol = n)

min.k.dists <- function(x,k=5) {
  apply(x,2,function(r) {
    b <- colSums((x - r)^2)
    o <- order(b)
    o[1:k]
  })
}

min.k.dists(x)  # first row should be 1:ncol(x); these points have distance 0
dist(t(x))      # can check answer against this

如果有人担心如何处理关系和其他事情，也许应该加入rank（）函数。

上面的代码似乎相当快，但我确定它可以改进（虽然我没有时间去使用C或Fortran）。因此，我仍然希望能够快速且稀疏地实现上述内容。

下面是我最终使用的并行化版本：

min.k.dists <- function(x,k=5,cores=1) {
  require(multicore)
  xx <- as.list(as.data.frame(x))
  names(xx) <- c()
  m <- mclapply(xx,function(r) {
    b <- colSums((x - r)^2)
    o <- order(b)
    o[1:k]
  },mc.cores=cores)
  t(do.call(rbind,m))
}

- Christopher DuBois

---

你需要执行 dist(t(x)) 来获得可比较的答案。 - Tommy

Answer 3

如果您想保留 min.k.dist 函数的逻辑并返回重复距离，您可能需要考虑对其进行一些修改。返回第一行的 0 距离似乎是毫无意义的，对吧？...通过结合我其他答案中的一些技巧，您可以将您的版本加速约 30%：

min.k.dists2 <- function(x, k=4L) {
  k <- max(2L, k + 1L)
  apply(x, 2, function(r) {
    sort.list(colSums((x - r)^2), na.last=NA, method='quick')[2:k]
  })
}

> n<-1e4; m<-3; m=matrix(runif(n*m), n)
> system.time(d <- min.k.dists(t(m), 4)) #To get 3 nearest neighbours and itself
   user  system elapsed 
  17.26    0.00   17.30 
> system.time(d <- min.k.dists2(t(m), 3)) #To get 3 nearest neighbours
   user  system elapsed 
   12.7     0.0    12.7 

- Tommy