Computing matriz de distancia dispersa por pares en R

Question

Tengo una matriz NxM y quiero calcular la matriz NxN de distancias euclidianas entre los puntos M. En mi problema, N es de aproximadamente 100,000. Como planeo usar esta matriz para un algoritmo vecino k-más cercano, solo necesito mantener las distancias más pequeñas de k, por lo que la matriz NxN resultante es muy escasa. Esto está en contraste con lo que sale de dist(), por ejemplo, que daría como resultado una matriz densa (y probablemente problemas de almacenamiento para mi tamaño N).

Los paquetes para kNN que he encontrado hasta ahora (knnflex, kknn, etc.) parecen usar matrices densas. Además, el paquete Matrix no ofrece una función de distancia pairwise.

Más cerca de mi objetivo, veo que el paquete spam tiene una función que permite tener en cuenta solo las distancias menores que cierto umbral, delta. En mi caso, sin embargo, un valor particular de delta puede producir demasiadas distancias (por lo que tengo que almacenar la matriz NxN densamente) o muy pocas distancias (de modo que no puedo usar kNN).

He visto comentarios previos sobre intentar realizar k-means clustering utilizando los paquetes bigmemory/biganalytics, pero no parece que pueda aprovechar estos métodos en este caso.

¿Alguien sabe una función/implementación que calculará una matriz de distancia de forma dispersa en R? Mi (temido) plan de copia de seguridad es tener dos for bucles y guardar resultados en un objeto Matrix.

Answer 1

Bueno, no puede tener que recurrir a la para-loops, ahora podemos :)

Hay, por supuesto la pregunta de cómo representar la matriz dispersa. Una forma simple es hacer que solo contenga los índices de los puntos más cercanos (y recalcular según sea necesario). Pero en la solución por debajo de, puse tanto la distancia ('d1' etc) y el índice ('i1' etc) en una sola matriz:

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 
} 

Tratando a cabo en 9 2d-puntos:

> 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 

Y probándolo en un problema mayor (10k puntos). Aún así, en 100k puntos y más dimensiones tomará mucho tiempo (como 15-30 minutos).

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

P.S. Solo noté que publicó una respuesta mientras escribía esto: la solución aquí es aproximadamente el doble de rápida porque no calcula la misma distancia dos veces (la distancia entre los puntos 1 y 13 es la misma que entre los puntos 13 y 1).

Answer 2

Por ahora estoy usando lo siguiente, inspirado en this answer. La salida es una matriz n x k donde el elemento (i,k) es el índice del punto de datos que es el k th más cercano a 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 

Si uno está preocupado acerca de cómo se manejan los lazos y todo eso, tal vez rank() deberían incorporarse.

El código anterior parece algo rápido, pero estoy seguro de que podría mejorarse (aunque no tengo tiempo para ir a la ruta C o fortran). Así que todavía estoy abierto a implementaciones rápidas y dispersas de lo anterior.

A continuación incluyo una versión parallelized que terminé usando:

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

Si desea mantener la lógica de su función min.k.dist y devolver distancias duplicadas, le recomendamos que la modifique un poco. Parece inútil devolver la primera línea con 0 de distancia, ¿verdad? ... y al incorporar algunos de los trucos en mi otra respuesta, puedes acelerar tu versión en un 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