Fastest pairwise distance metric in python

Question

I have an 1D array of numbers, and want to calculate all pairwise euclidean distances. I have a method (thanks to SO) of doing this with broadcasting, but it's inefficient because it calculates each distance twice. And it doesn't scale well.

Here's an example that gives me what I want with an array of 1000 numbers.

import numpy as np import random r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)]) dists = np.abs(r - r[:, None])

What's the fastest implementation in scipy/numpy/scikit-learn that I can use to do this, given that it has to scale to situations where the 1D array has >10k values.

Note: the matrix is symmetric, so I'm guessing that it's possible to get at least a 2x speedup by addressing that, I just don't know how.

There's a function for that: scipy.spatial.distance.pdist. I dunno whether this is the fastest option, since it needs to have checks for multidimensional data, non-Euclidean norms, and other things, but it's built in. — user2357112
– user2357112, Commented Nov 29, 2013 at 3:51
How fast do you need this to be? It's never going to scale better than O(n^2), since you have to populate n^2 entries of output. Your existing solution is O(n^2), and there doesn't seem to be much room for major optimizations. — user2357112
– user2357112, Commented Nov 29, 2013 at 3:56
This seems to scale to >10k values well enough already when I try it. Remember that you need to populate 100 million entries of output. That's almost half a gigabyte of pairwise distances. — user2357112
– user2357112, Commented Nov 29, 2013 at 4:10
@askewchan I don't think it does... If you follow the source code, in the end this is the function getting called. Not only is there no fancy optimization, but for 1D vectors it is squaring and taking the square root to compute the absolute value. Probably worse than the OP's code for his particular use case. — Jaime
– Jaime, Commented Nov 29, 2013 at 5:45
@CTZhu If I'm not mistaken, scipy is always compiled with BLAS, it's not optional as with numpy. — askewchan
– askewchan, Commented Nov 29, 2013 at 5:46

roblanf · Accepted Answer · 2013-12-02 00:59:10Z

Neither of the other answers quite answered the question - 1 was in Cython, one was slower. But both provided very useful hints. Following up on them suggests that scipy.spatial.distance.pdist is the way to go.

Here's some code:

import numpy as np import random import sklearn.metrics.pairwise import scipy.spatial.distance r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)]) c = r[:, None] def option1(r): dists = np.abs(r - r[:, None]) def option2(r): dists = scipy.spatial.distance.pdist(r, 'cityblock') def option3(r): dists = sklearn.metrics.pairwise.manhattan_distances(r)

Timing with IPython:

In [36]: timeit option1(r) 100 loops, best of 3: 5.31 ms per loop In [37]: timeit option2(c) 1000 loops, best of 3: 1.84 ms per loop In [38]: timeit option3(c) 100 loops, best of 3: 11.5 ms per loop

I didn't try the Cython implementation (I can't use it for this project), but comparing my results to the other answer that did, it looks like scipy.spatial.distance.pdist is roughly a third slower than the Cython implementation (taking into account the different machines by benchmarking on the np.abs solution).

I assume this is as fast as: scikit-learn.org/stable/modules/generated/… ? the version in sklearn?
In the scipy.spatial.distance.pdist case, should it be c instead of r?

Saullo G. P. Castro · Accepted Answer · 2013-11-29 11:08:46Z

Here is a Cython implementation that gives more than 3X speed improvement for this example on my computer. This timing should be reviewed for bigger arrays tough, because the BLAS routines can probably scale much better than this rather naive code.

I know you asked for something inside scipy/numpy/scikit-learn, but maybe this will open new possibilities for you:

File my_cython.pyx:

import numpy as np cimport numpy as np import cython cdef extern from "math.h": double abs(double t) @cython.wraparound(False) @cython.boundscheck(False) def pairwise_distance(np.ndarray[np.double_t, ndim=1] r): cdef int i, j, c, size cdef np.ndarray[np.double_t, ndim=1] ans size = sum(range(1, r.shape[0]+1)) ans = np.empty(size, dtype=r.dtype) c = -1 for i in range(r.shape[0]): for j in range(i, r.shape[0]): c += 1 ans[c] = abs(r[i] - r[j]) return ans

The answer is a 1-D array containing all non-repeated evaluations.

To import into Python:

import numpy as np import random import pyximport; pyximport.install() from my_cython import pairwise_distance r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)], dtype=float) def solOP(r): return np.abs(r - r[:, None])

Timing with IPython:

In [2]: timeit solOP(r) 100 loops, best of 3: 7.38 ms per loop In [3]: timeit pairwise_distance(r) 1000 loops, best of 3: 1.77 ms per loop

cyborg · Accepted Answer · 2013-11-30 23:38:41Z

Using half the memory, but 6 times slower than np.abs(r - r[:, None]):

triu = np.triu_indices(r.shape[0],1) dists2 = abs(r[triu[1]]-r[triu[0]])

Collectives™ on Stack Overflow

Fastest pairwise distance metric in python

3 Answers 3

3 Comments

1 Comment

Comments

Linked

Hot Network Questions