21

我有一个NxM矩阵,我想计算点NxN之间的欧几里得距离矩阵M。在我的问题中,N大约是 100,000。由于我计划将此矩阵用于 k 近邻算法,因此我只需要保持k最小距离,因此得到的NxN矩阵非常稀疏。这与来自 的结果形成对比dist(),例如,这将导致密集矩阵(并且对于我的 size 可能存在存储问题N)。

到目前为止,我发现的 kNN 包(knnflexkknn等)似乎都使用了密集矩阵。此外,该Matrix软件包不提供成对距离功能。

更接近我的目标,我看到该spam软件包具有一项nearest.dist()功能,该功能允许人们仅考虑小于某个阈值的距离,delta. 然而,在我的例子中,一个特定的值delta可能会产生太多的距离(所以我必须NxN密集地存储矩阵)或太少的距离(所以我不能使用 kNN)。

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

有人知道在 R 中以稀疏方式计算距离矩阵的函数/实现吗?我(可怕的)备份计划是有两个for循环并将结果保存在一个Matrix对象中。

4

3 回答 3

7

好吧,我们不能让你诉诸 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 个 2d 点上试一试:

> 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...

PS刚刚注意到你在我写这篇文章时发布了一个答案:这里的解决方案大约快两倍,因为它不会两次计算相同的距离(点1和13之间的距离与点13和1之间的距离相同) .

于 2011-04-06T16:19:07.147 回答
3

现在我正在使用以下内容,受此答案的启发。输出是一个n x k矩阵,其中 element(i,k)k最接近的数据点的索引i

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()应该合并。

上面的代码看起来有点快,但我确信它可以改进(虽然我没有时间去Corfortran路线)。所以我仍然对上述的快速和稀疏的实现持开放态度。

下面我包含了一个我最终使用的并行版本:

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))
}
于 2011-04-06T16:03:00.573 回答
1

如果您想保留 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 
于 2011-04-07T15:54:59.147 回答