# How can I optimize this pairwise Euclidean distance search algorithm?

• x and y are two-dimensional arrays with dimensions (AxN) and (BxN), i.e. they have the same number of columns.

• I need to get a matrix of Euclidean distances between each pair of rows from x and y.

I have already written four different methods that give equally good results. But the first three methods are not fast enough, and the fourth method consumes too much memory.

1:

from scipy.spatial import distance

def euclidean_distance(x, y):
return distance.cdist(x, y)


2:

import numpy as np
import math

def euclidean_distance(x, y):
res = []
one_res = []
for i in range(len(x)):
for j in range(len(y)):
one_res.append(math.sqrt(sum(x[i] ** 2) + sum(y[j] ** 2) - 2 *
np.dot(x[i], np.transpose(y[j]))))
res.append(one_res)
one_res = []
return res


3:

import numpy as np

def euclidean_distance(x, y):
distances = np.zeros((x.T.shape, y.T.shape))
for x, y in zip(x.T, y.T):
distances = np.subtract.outer(x, y) ** 2 + distances
return np.sqrt(distances)


4:

import numpy as np

def euclidean_distance(x, y):
x = np.expand_dims(x, axis=1)
y = np.expand_dims(y, axis=0)
return np.sqrt(((x - y) ** 2).sum(-1))


It looks like the second way could still be improved, and it also shows the best timing. Can you help me with the optimization, please?

• Welcome to Code Review@SE. Time allowing, (re?)visit How to get the best value out of Code Review. As per How do I ask a Good Question?, the title of a CR question (and, preferably, the code itself) should state what the code is to accomplish. Oct 28, 2021 at 1:40
• It may help the non-mathheads if you added why&how 2-4 compute the distances. Please add how the result is going to be used - e.g, for comparisons, you don't need sqrt (cdist(x, y, 'sqeuclidean')). Do you need all values at once/as a matrix? Otherwise, a generator may help. What are typical dimensions of x and y? Oct 28, 2021 at 2:04

### Performance:

A way to use broadcasting but avoid allocating the 3d array is to separate the distance formula into $(x - y) ** 2 = x^2 + y^2 - 2xy$.

squared_matrix = np.sum(A ** 2,axis=1 )[:,np.newaxis] - 2 * (A@B.T) + np.sum(B ** 2,axis=1 )[np.newaxis,:]
distance_matrix = np.sqrt(squared_matrix)


Check:

np.allclose(distance_matrix,euclidean_distance(A,B))
>>> True

• You can probably write this with np.einsum in a more succinct way but einsum makes my head hurt. Oct 29, 2021 at 12:46

Ideas for 2:

• preallocate res
• don't use one_res, just a comprehensions