Oops. Looks like I forgot to attach the test program that generated
that output so you can tell what dist2g actually does.
Funny thing is -- despite being written in C, hypot doesn't actually
win any of the test cases for which it's applicable.
--bb
On 4/17/07, Bill Baxter <[EMAIL PROTECTED]> wrote:
Here's a bunch of dist matrix implementations and their timings.
The upshot is that for most purposes this seems to be the best or at
least not too far off (basically the cookbook solution Kier posted)
def dist2hd(x,y):
"""Generate a 'coordinate' of the solution at a time"""
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
for i in xrange(x.shape[1]):
diff2 = x[:,i,None] - y[:,i]
diff2 **= 2
d += diff2
npy.sqrt(d,d)
return d
The only place where it's far from the best is for a small number of
points (~10) with high dimensionality (~100), which does come up in
machine learning contexts. For those cases this does much better
(factor of :
def dist2b3(x,y):
d = npy.dot(x,y.T)
d *= -2.0
d += (x*x).sum(1)[:,None]
d += (y*y).sum(1)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
# in place sqrt
npy.sqrt(d,d)
return d
So given that, the obvious solution (if you don't want to delve into
non-numpy solutions) is to use a hybrid that just switches between
the two. Not sure what the proper switch is since it seems kind of
complicated, and probably depends some on cache specifics. But just
switching based on the dimension of the points seems to be pretty
effective:
def dist2hy(x,y):
if x.shape[1]<5:
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
for i in xrange(x.shape[1]):
diff2 = x[:,i,None] - y[:,i]
diff2 **= 2
d += diff2
npy.sqrt(d,d)
return d
else:
d = npy.dot(x,y.T)
d *= -2.0
d += (x*x).sum(1)[:,None]
d += (y*y).sum(1)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
# in place sqrt
npy.sqrt(d,d)
return d
All of this assumes 'C' contiguous data. All bets are off if you have
non-contiguous or 'F' ordered data. And maybe if x and y have very
different numbers of points.
--bb
On 4/17/07, Keir Mierle <[EMAIL PROTECTED]> wrote:
> On 4/13/07, Timothy Hochberg <[EMAIL PROTECTED]> wrote:
> > On 4/13/07, Bill Baxter <[EMAIL PROTECTED]> wrote:
> > > I think someone posted some timings about this before but I don't recall.
>
> http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/498246
>
> [snip]
> > I'm going to go out on a limb and contend, without running any timings, that
> > for large M and N, a solution using a for loop will beat either of those.
> > For example (untested):
> >
> > results = empty([M, N], float)
> > # You could be fancy and swap axes depending on which array is larger, but
> > # I'll leave that for someone else
> > for i, v in enumerate(x):
> > results[i] = sqrt(sum((v-y)**2, axis=-1))
> > Or something like that. The reason that I suspect this will be faster is
> > that it has better locality, completely finishing a computation on a
> > relatively small working set before moving onto the next one. The one liners
> > have to pull the potentially large MxN array into the processor repeatedly.
>
> In my experience, it is indeed the case that the for loop version is
> faster. The fastest of the three versions offered in the above url is
> the last:
>
> from numpy import mat, zeros, newaxis
> def calcDistanceMatrixFastEuclidean2(nDimPoints):
> nDimPoints = array(nDimPoints)
> n,m = nDimPoints.shape
> delta = zeros((n,n),'d')
> for d in xrange(m):
> data = nDimPoints[:,d]
> delta += (data - data[:,newaxis])**2
> return sqrt(delta)
>
> This is easily extended to two different nDimPoints matricies.
>
> Cheers,
> Keir
> _______________________________________________
> Numpy-discussion mailing list
> [email protected]
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
"""
This tests the performance of several distance matrix implementations
"""
import numpy as npy
import timeit
def dist2a(x,y):
d = (x*x).sum(1)[:,None] + (y*y).sum(1) - 2*npy.dot(x,y.T)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
d = npy.sqrt(d)
return d
def dist2b(x,y):
d = (x*x).sum(1)[:,None] + (y*y).sum(1) - 2*npy.dot(x,y.T)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
# in place sqrt
npy.sqrt(d,d)
return d
def dist2b2(x,y):
d = -2.0 * npy.dot(x,y.T)
d += (x*x).sum(1)[:,None]
d += (y*y).sum(1)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
# in place sqrt
npy.sqrt(d,d)
return d
def dist2b3(x,y):
d = npy.dot(x,y.T)
d *= -2.0
d += (x*x).sum(1)[:,None]
d += (y*y).sum(1)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
# in place sqrt
npy.sqrt(d,d)
return d
def dist2c(x,y):
d = npy.sqrt( npy.sum( (x[:,None]-y)**2, axis=-1) )
return d
def dist2c1(x,y):
d = npy.sqrt( ((x[:,None]-y)**2).sum(-1) )
return d
def dist2c2(x,y):
diff2 = x[:,None]-y
diff2 **= 2
return npy.sqrt( diff2.sum(-1) )
def dist2d(x,y):
N = x.shape[0]
d = npy.sum( (x[:,None]-y)**2, axis=-1)
# sqrt in place
npy.sqrt(d,d)
return d
def dist2e(x,y):
"""Generate a column of the output at a time"""
d = npy.empty((x.shape[0],y.shape[0]),x.dtype)
for i,yi in enumerate(y):
d[:,i] = ((x - yi)**2).sum(axis=1)
return npy.sqrt(d)
def dist2f(x,y):
"""Generate a row of the output at a time"""
d = npy.empty((x.shape[0],y.shape[0]),x.dtype)
for i,xi in enumerate(x):
d[i] = ((xi - y)**2).sum(axis=1)
return npy.sqrt(d)
def dist2g(x,y):
"""Use the built-in hypot function"""
d = x[:,None]-y
d = npy.hypot(d[...,0],d[...,1])
# this creates some -eps values, (and only works for 2d data)
d[d<0] = 0
return d
def dist2h(x,y):
"""Generate a 'coordinate' of the solution at a time"""
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
for i in xrange(x.shape[1]):
diff = x[:,i,None] - y[:,i]
d += diff**2
return npy.sqrt(d)
def dist2hb(x,y):
"""Generate a 'coordinate' of the solution at a time"""
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
for i in xrange(x.shape[1]):
diff = x[:,i,None] - y[:,i]
d += diff**2
npy.sqrt(d,d)
return d
def dist2hc(x,y):
"""Generate a 'coordinate' of the solution at a time"""
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
y = y.copy('f')
for i in xrange(x.shape[1]):
diff = x[:,i,None] - y[:,i]
d += diff**2
npy.sqrt(d,d)
return d
def dist2hd(x,y):
"""Generate a 'coordinate' of the solution at a time"""
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
for i in xrange(x.shape[1]):
diff2 = x[:,i,None] - y[:,i]
diff2 **= 2
d += diff2
npy.sqrt(d,d)
return d
def dist2hy(x,y):
"""Hybrid of two good solutions"""
if x.shape[1]<5:
d = npy.zeros((x.shape[0],y.shape[0]),dtype=x.dtype)
for i in xrange(x.shape[1]):
diff2 = x[:,i,None] - y[:,i]
diff2 **= 2
d += diff2
npy.sqrt(d,d)
return d
else:
d = npy.dot(x,y.T)
d *= -2.0
d += (x*x).sum(1)[:,None]
d += (y*y).sum(1)
# Rounding errors occasionally cause negative entries in d
d[d<0] = 0
# in place sqrt
npy.sqrt(d,d)
return d
def main():
fns = ['dist2a',
#'dist2b',
#'dist2b2',
'dist2b3',
#'dist2c',
#'dist2c1',
'dist2c2',
'dist2d',
'dist2e',
'dist2f',
'dist2g',
#'dist2h',
#'dist2hb',
#'dist2hc',
'dist2hd',
'dist2hy']
from itertools import izip
for N,trials in izip((10,100,1000,3000),(10000,500,50,10)):
dlist = [1,2,3]
if N<1000: dlist += [5,10]
if N<100: dlist += [100]
for d in dlist:
print "DIM=%d, N = %d ------------- " % (d,N)
timers = []
for f in fns:
if f=='dist2g' and d!=2:
timers.append( None )
continue
timers.append(
timeit.Timer('d = distfn(x,y)',
'from __main__ import %s as distfn; import numpy as npy; npy.random.seed(1234); d=%d; N=%d; x = npy.round(20*npy.random.rand(N,d)-10); y = npy.round(20*npy.random.rand(N,d)-10);' %(f,d,N)))
result=[]
for t in timers:
try:
result += [t.timeit(trials)]
except:
result += [npy.inf]
idx = npy.argsort(result)
bestt = result[idx[0]]
for i in idx:
print "Func",fns[i],": \t", result[i], '\t%0.3f'%(bestt/result[i])
if __name__=="__main__":
main()
_______________________________________________
Numpy-discussion mailing list
[email protected]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
