4
\$\begingroup\$
  • 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[0].shape[0], y.T[0].shape[0]))
    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?

\$\endgroup\$
2
  • 2
    \$\begingroup\$ 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. \$\endgroup\$
    – greybeard
    Oct 28, 2021 at 1:40
  • 3
    \$\begingroup\$ 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? \$\endgroup\$
    – greybeard
    Oct 28, 2021 at 2:04

3 Answers 3

3
\$\begingroup\$

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 * ([email protected]) + np.sum(B ** 2,axis=1 )[np.newaxis,:]
distance_matrix = np.sqrt(squared_matrix)

Check:

np.allclose(distance_matrix,euclidean_distance(A,B))
>>> True
\$\endgroup\$
1
  • 1
    \$\begingroup\$ You can probably write this with np.einsum in a more succinct way but einsum makes my head hurt. \$\endgroup\$
    – kubatucka
    Oct 29, 2021 at 12:46
2
\$\begingroup\$

sklearn's library Euclidean distance calculation was tested to perform better than scipy & numpy due to used vectorisation implementation here . Codes were compared:

import numpy as np
from scipy.spatial.distance import cdist
from sklearn.metrics.pairwise import euclidean_distances

vec = np.array([ [ 0, 1 ], [ 1, 0 ],
            [ 1, 2 ], [ 2, 1 ] ])  
######################## 1
scipy_cdist = cdist(vec, vec, metric='euclidean')
print(scipy_cdist)

######################## 2
sklearn_dist = euclidean_distances(vec, vec)
print(sklearn_dist)

########################  3  
#  method consuming memory.
x = np.sum(vec**2, axis=1)[:, np.newaxis]
y = np.sum(vec**2,    axis=1)
xy = np.dot(vec, vec.T)
dist = np.sqrt(x + y - 2*xy)
print(dist)

######################## 4
# Without pre-allocating memory
dist = []
for i in range(len(vec)):
    dist.append(((vec- vec[i])**2).sum(axis=1)**0.5)
print("Without pre-allocating memory:\n", np.array(dist, dtype= float))

# pre-allocating memory
D = np.empty((len(vec),len(vec)))
for i in range(len(vec)):
    D[i, :] = ((vec-vec[i])**2).sum(axis=1)**0.5
print("pre-allocating memory:\n", D)

plus preallocation of memory can influence the performance... can see by link testing_code to compare all these 4 approaches

\$\endgroup\$
4
  • 1
    \$\begingroup\$ This is missing norm(). Also, for most intents, anything linked doesn't exist. You need to paste times and analyses here. \$\endgroup\$
    – Reinderien
    Jun 30 at 17:51
  • \$\begingroup\$ why norm() is missing? Euclidean is 2-norm. I corrected the code to its simpliest version to be reproducible, think that adding more complicated multidimensional array & timing will not become too hard for unbelievers. The fourth approach with memory-preallocation I'm not going to test myself - i don't need it, as I believe to the link (I also checked) -- but this can be the advice for topic-starter if he still wants to save memory - or you can invest your own time if you doubt about the conclusion I've already done from the article - sklearn performs the best \$\endgroup\$
    – JeeyCi
    Jun 30 at 18:20
  • \$\begingroup\$ BTW, scipy's debugging messages seems more explaining than sklearn's ones, therefore the choice of the library is developer's choice anyway \$\endgroup\$
    – JeeyCi
    Jun 30 at 19:10
  • \$\begingroup\$ added memory_preallocation from link's example ... but still it is a little worth \$\endgroup\$
    – JeeyCi
    Jul 1 at 16:11
1
\$\begingroup\$

Ideas for 2:

  • preallocate res
  • don't use one_res, just a comprehensions
\$\endgroup\$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.