Calcul de la matrice des distances par paires en R
j'ai une matrice NxM et je veux calculer la matrice NxN des distances euclidiennes entre les points M . Dans mon problème, N est d'environ 100 000. Comme j'ai l'intention d'utiliser cette matrice pour un algorithme de K-plus proche voisin, je n'ai besoin de garder les k plus petites distances, de sorte que la matrice NxN résultant est très clairsemée. Ceci est en contraste avec ce qui sort de dist() , par exemple, qui aboutirait à une matrice dense (et probablement problèmes de stockage pour ma taille N ).
les paquets pour kNN que j'ai trouvés jusqu'à présent ( knnflex , kknn , etc) semblent tous utiliser des matrices denses. De plus, le paquet Matrix n'offre pas de fonction de distance en paires.
plus près de mon objectif, je vois que le paquet spam a une fonction nearest.dist() qui permet de ne considérer que les distances inférieures à un certain seuil, delta . Dans mon cas, cependant, un une valeur particulière de delta peut produire trop de distances (de sorte que je dois stocker la matrice NxN densément) ou trop peu de distances (de sorte que je ne peux pas utiliser kNN).
j'ai déjà vu des discussions sur le fait d'essayer d'effectuer k-means clustering en utilisant les paquets bigmemory/biganalytics , mais il ne semble pas que je puisse tirer parti de ces méthodes dans ce cas.
est-ce que quelqu'un connaît une fonction/implémentation qui va calculer un matrice de distance de façon éparse en R? Mon plan de sauvegarde (redouté) est d'avoir deux for boucles et enregistrer des résultats dans un Matrix objet.
3 réponses
Eh bien, nous ne pouvons pas vous avoir recours à for-loops, maintenant pouvons-nous:)
il y a bien sûr la question de savoir comment représenter la matrice clairsemée. Une façon simple est de ne contenir que les indices des points les plus proches (et de recalculer au besoin). Mais dans la solution ci-dessous, j'ai mis à la fois la distance ('d1' etc) et l'indice ('i1' etc) dans une matrice unique:
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
}
l'Essayer sur 9 2d-points:
> 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
Et l'essayer sur un plus grand problème (10k points). Néanmoins, sur 100k points et plus de dimensions, cela prendra beaucoup de temps (comme 15-30 minutes).
n<-1e4; m<-3; m=matrix(runif(n*m), n)
system.time( d <- sparseDist(m, 3) ) # 9 seconds on my machine...
P. vient de noter que vous avez posté une réponse pendant que j'écrivais ceci: la solution ici est environ deux fois plus rapide parce qu'elle ne calcule pas la même distance deux fois (la distance entre les points 1 et 13 est la même qu'entre les points 13 et 1).
pour l'instant j'utilise ce qui suit, inspiré par cette réponse . Le résultat est une matrice n x k où l'élément (i,k) est l'indice du point de données qui est le k e le plus proche de 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 l'on s'inquiète de la façon dont les liens sont traités et autres, peut-être que rank() devrait être incorporé.
le code ci-dessus semble un peu rapide, mais je suis sûr qu'il pourrait être amélioré (bien que je n'ai pas le temps d'aller sur la route C ou fortran ). Je suis donc toujours ouvert aux implémentations rapides et éparses de ce qui précède.
ci-dessous, j'inclus une version parallélisée que j'ai fini par utiliser:
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))
}
Si vous voulez garder la logique de votre min.K. fonction dist et retour des distances dupliquées, vous voudrez peut-être envisager de le modifier un peu. Il semble inutile de retourner la première ligne avec 0 distance, non? ...et en incorporant certains des trucs dans mon autre réponse, vous pouvez accélérer votre version d'environ 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