Computing sparse pairwise distance matrix in R

Question

I have a NxM matrix and I want to compute the NxN matrix of Euclidean distances between the M points. In my problem, N is about 100,000. As I plan to use this matrix for a k-nearest neighbor algorithm, I only need to keep the k smallest distances, so the resulting NxN matrix is very sparse. This is in contrast to what comes out of dist(), for example, which would result in a dense matrix (and probably storage problems for my size N).

The packages for kNN that I've found so far (knnflex, kknn, etc) all appear to use dense matrices. Also, the Matrix package does not offer a pairwise distance function.

Closer to my goal, I see that the spam package has a nearest.dist() function that allows one to only consider distances less than some threshold, delta. In my case, however, a particular value of delta may produce too many distances (so that I have to store the NxN matrix densely) or too few distances (so that I can't use kNN).

I have seen previous discussion on trying to perform k-means clustering using the bigmemory/biganalytics packages, but it doesn't seem like I can leverage these methods in this case.

Does anybody know a function/implementation that will compute a distance matrix in a sparse fashion in R? My (dreaded) backup plan is to have two for loops and save results in a Matrix object.

Answer 1

Well, we can't have you resorting to for-loops, now can we :)

There is of course the question of how to represent the sparse matrix. A simple way is to have it only contain the indices of the points that are closest (and recalculate as needed). But in the solution below, I put both distance ('d1' etc) and index ('i1' etc) in a single matrix:

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
}

Trying it out on 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

And trying it on a larger problem (10k points). Still, on 100k points and more dimensions it will take a long time (like 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.S. Just noted that you posted an answer as I was writing this: the solution here is roughly twice as fast because it doesn't calculate the same distance twice (the distance between points 1 and 13 is the same as between points 13 and 1).