Geek Answers Handbook

Euclidean norm using numexpr

I need to rewrite this code with numexpr, calculating the Euclidean normal matrix data matrix [rows x cols] and vector [1 x cols].

d = ((data-vec)**2).sum(axis=1)

How can I do that? Perhaps there is another faster method?

The problem I'm using hdf5 and the data received from it. For example, this code gives an error: objects are not aligned.

 #naive numpy solution, can be parallel? def test_bruteforce_knn(): h5f = tables.open_file(fileName) t0= time.time() d = np.empty((rows*batches,)) for i in range(batches): d[i*rows:(i+1)*rows] = ((h5f.root.carray[i*rows:(i+1)*rows]-vec)**2).sum(axis=1) print (time.time()-t0) ndx = d.argsort() print ndx[:k] h5f.close() #using some tricks (don't work error: objects are not aligned ) def test_bruteforce_knn(): h5f = tables.open_file(fileName) t0= time.time() d = np.empty((rows*batches,)) for i in range(batches): d[i*rows:(i+1)*rows] = (np.einsum('ij,ij->i', h5f.root.carray[i*rows:(i+1)*rows], h5f.root.carray[i*rows:(i+1)*rows]) + np.dot(vec, vec) -2 * np.dot(h5f.root.carray[i*rows:(i+1)*rows], vec)) print (time.time()-t0) ndx = d.argsort() print ndx[:k] h5f.close()

Using numexpr: it seems that numexpr does not understand h5f.root.carray [i * rows: (i + 1) * rows] should it be reassigned?

 import numexpr as ne def test_bruteforce_knn(): h5f = tables.open_file(fileName) t0= time.time() d = np.empty((rows*batches,)) for i in range(batches): ne.evaluate("sum((h5f.root.carray[i*rows:(i+1)*rows] - vec) ** 2, axis=1)") print (time.time()-t0) ndx = d.argsort() print ndx[:k] h5f.close()
+1
source share
1 answer

There's a potentially quick way (for very large arrays) using only NumPy, which is used in scikit-learn:

 def squared_row_norms(X): # From http://stackoverflow.com/q/19094441/166749 return np.einsum('ij,ij->i', X, X) def squared_euclidean_distances(data, vec): data2 = squared_row_norms(data) vec2 = squared_row_norms(vec) d = np.dot(data, vec.T).ravel() d *= -2 d += data2 d += vec2 return d

This is based on the fact that (x - y) ² = x² + y² - 2xy, even for vectors.

Test:

 >>> data = np.random.randn(10, 40) >>> vec = np.random.randn(1, 40) >>> ((data - vec) ** 2).sum(axis=1) array([ 96.75712686, 69.45894306, 100.71998244, 80.97797154, 84.8832107 , 82.28910021, 67.48309433, 81.94813371, 64.68162331, 77.43265692]) >>> squared_euclidean_distances(data, vec) array([ 96.75712686, 69.45894306, 100.71998244, 80.97797154, 84.8832107 , 82.28910021, 67.48309433, 81.94813371, 64.68162331, 77.43265692]) >>> from sklearn.metrics.pairwise import euclidean_distances >>> euclidean_distances(data, vec, squared=True).ravel() array([ 96.75712686, 69.45894306, 100.71998244, 80.97797154, 84.8832107 , 82.28910021, 67.48309433, 81.94813371, 64.68162331, 77.43265692])

Profile:

 >>> data = np.random.randn(1000, 40) >>> vec = np.random.randn(1, 40) >>> %timeit ((data - vec)**2).sum(axis=1) 10000 loops, best of 3: 114 us per loop >>> %timeit squared_euclidean_distances(data, vec) 10000 loops, best of 3: 52.5 us per loop

Using numexpr is also possible, but it does not seem to give acceleration for 1000 points (and by 10000 it is not much better):

 >>> %timeit ne.evaluate("sum((data - vec) ** 2, axis=1)") 10000 loops, best of 3: 142 us per loop
+4
source

More articles:


All Articles