Вычисление разреженной матрицы попарных расстояний в R

Question

У меня есть матрица NxM, и я хочу вычислить матрицу NxN евклидовых расстояний между M точками. В моей проблеме N составляет около 100 000. Поскольку я планирую использовать эту матрицу для алгоритма k-ближайшего соседа, мне нужно только сохранить наименьшее расстояние k, поэтому полученная матрица NxN будет очень разреженной. Это в отличие от того, что получается, например, dist(), что может привести к плотной матрице (и, вероятно, к проблемам с хранилищем для моего размера N).

Пакеты для kNN, которые я нашел до сих пор (knnflex, kknn и т. Д.), Похоже, используют плотные матрицы. Кроме того, пакет Matrix не предлагает функцию попарного расстояния.

Ближе к моей цели, я вижу, что пакет spam имеет функцию nearest.dist(), которая позволяет учитывать только расстояния, меньшие некоторого порога delta. В моем случае, однако, конкретное значение delta может привести к слишком большому количеству расстояний (так что я должен плотно хранить матрицу NxN) или к слишком небольшим расстояниям (чтобы я не мог использовать kNN).

Я видел предыдущее обсуждение попытки выполнить кластеризацию k-средних с использованием пакетов bigmemory/biganalytics, но не похоже, что я могу использовать эти методы в этом случае.

Кто-нибудь знает функцию / реализацию, которая будет редко вычислять матрицу расстояний в R? Мой (страшный) план резервного копирования должен иметь два цикла for и сохранять результаты в объекте Matrix.

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

И пробовал это на более крупной проблеме (10 тыс. Баллов). Тем не менее, на 100 тыс. Точек и более измерениях это займет много времени (например, 15-30 минут).

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

P.S. Только что отметил, что вы отправили ответ, как я писал: решение здесь примерно в два раза быстрее, потому что оно не рассчитывает одинаковое расстояние дважды (расстояние между точками 1 и 13 такое же, как между точками 13 и 1).

Answer 2

Пока я использую следующее, вдохновленное этим ответом . Выходные данные представляют собой матрицу n x k, где элемент (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().

Приведенный выше код кажется довольно быстрым, но я уверен, что его можно улучшить (хотя у меня нет времени идти по маршруту 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))
}

Answer 3

Если вы хотите сохранить логику своей функции min.k.dist и вернуть дубликаты расстояний, вы можете немного ее изменить. Кажется бессмысленным возвращать первую строку с нулевым расстоянием, верно? ... и включив некоторые уловки в мой другой ответ, вы можете ускорить вашу версию примерно на 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