I'm looking for resources on fast, numerically stable pairwise euclidean distance algorithms. In particular, suppose $A \in \mathbb{R}^{M \times D}$ and $B \in \mathbb{R}^{N \times D}$ are two sets of row vectors. I would like to compute the matrix,

$$X \in \mathbb{R}^{M \times N}, \quad X_{i,j} = \|A_i - B_j\|_2^2.$$

So far I have found two methods:

Method 1 - The simple approach is to loop over each vector $A_i$ and each vector $B_j$ and populate $X$ with the corresponding sum of squares.

Method 2 - Take advantage of the fact that,

$$\|A_i - B_j\|_2^2 = \langle A_i - B_j, A_i - B_j \rangle = \|A_i\|_2^2 + \|B_j\|_2^2 - 2 \langle A_i, B_j \rangle.$$

Taking advantage of efficient code for computing each of these three terms, Method 2 is almost an order of magnitude faster. However, it is less numerically stable than Method 1. For example, Method 2 can output negative distances.

Are there alternate approaches to computing pairwise distances which are faster than Method 1 but guarantee (at the very least) all non-negative distances? have better numerical stability?

Code Demonstration - Below is an example comparison written in Python. Scipy's cdist() function is, effectively, an implementation of Method 1 whereas cdist_fast() below is an implementation of Method 2:

# experiment.py
import numpy as np
import time
from scipy.spatial.distance import cdist

def cdist_fast(XA, XB):
    XA_norm = np.sum(XA**2, axis=1)
    XB_norm = np.sum(XB**2, axis=1)
    XA_XB_T = np.dot(XA, XB.T)
    distances = XA_norm.reshape(-1,1) + XB_norm - 2*XA_XB_T
    return distances

def main():
    M,N = 5000, 128
    XA = np.random.randn(M,N)

    t = time.time()
    distances_cdist = cdist(XA, XA, metric='sqeuclidean')
    time_cdist = time.time() - t

    t = time.time()
    distances_cdist_fast = cdist_fast(XA, XA)
    time_cdist_fast = time.time() - t

    print(f'time_cdist = {time_cdist:.3f} s')
    print(f'time_cdist_fast = {time_cdist_fast:.3f} s')

    # check validity of results
    assert np.allclose(distances_cdist, distances_cdist_fast)

    # check that the results are non-negative
        assert (distances_cdist >= 0.0).all()
    except AssertionError:
        print('Numerical instability in cdist()')

        assert (distances_cdist_fast >= 0.0).all()
    except AssertionError:
        print('Numerical instability in cdist_fast()')

if __name__ == '__main__':

Script output on a 3.1 GHz Intel Core i7:

$ python experiment.py
time_cdist = 3.457 s
time_cdist_fast = 0.625 s
Numerical instability in cdist_fast()
  • $\begingroup$ Do you have a bound on the values A and B could take? $\endgroup$ Commented Oct 21, 2018 at 17:59
  • $\begingroup$ @TolgaBirdal - Not really. For now, let's just assume that A and B consist of either normally or uniformly distributed data at some more-or-less arbitrary scale. $\endgroup$ Commented Oct 22, 2018 at 18:12
  • $\begingroup$ Does scipy.spatial.cdist use multiple cores on your machine, e.g. with export OMP_NUM_THREADS=4 VECLIB_MAXIMUM_THREADS=4 ? $\endgroup$
    – denis
    Commented Nov 20, 2018 at 12:18

2 Answers 2


I would strongly advise against using the method 2: $$\|A_i - B_j\|_2^2 = \langle A_i - B_j, A_i - B_j \rangle = \|A_i\|_2^2 + \|B_j\|_2^2 - 2 \langle A_i, B_j \rangle.$$ whether you are uing the absolute value on top of it or not. Because you are going to get your positive (but wrong) distances; however, you are still lacking any numerical stability.

Now, for 5000 vectors in $\mathbb R^{128}$, your (scipy) stable version is only ~6 times slower than the fast one. I don't think you are able to find a numerically stable version that is faster than provided by scipy.cdist. I would make sure that it is actually a bottleneck in your code. Especially since method 1 naturally sums only positive entries, an ideal and cheap scenario to not deal with the stability!

You definitely can try using Kahan or Neumaier summation, but I would be surprised that if they are employed (note, Neumaier has to be applied for each entry individually), the aforementioned method 2 would be any faster than scipy.cdist, most likely, at least tens time slower.

However, if the positiveness of the distance is the only thing you care, you certainly can use an absolute value on top of your cdist_fast(.,.) function.

  • $\begingroup$ How do you distinguish round-off errors from "wrong" distances? Which knowledge do you lose when you replace the round-off error by its absolute value? $\endgroup$
    – user7440
    Commented Oct 17, 2018 at 16:42
  • $\begingroup$ Thank you for your response. I know that Numpy already implements pairwise summation which would at least reduce the floating point error in the calculation of the first two terms in "Method 2". It might be worth looking at these other two approaches. As for the slowdown, I admit that I haven't fully explored the problem size space. Depending on my application cdist() may be faster than I think. As for positiveness, the numerical error can sometimes produce multiple zero-distances even when the set XA consists of all distinct points. $\endgroup$ Commented Oct 17, 2018 at 17:57
  • $\begingroup$ @user7440 I don't mind having an absolute value on top of the $||A_i||_2+...$ expression. With or without absolute value, it's just not a stable way to do the calculation. So, you don't lose any knowledge, as far as I can see. $\endgroup$
    – Anton Menshov
    Commented Oct 18, 2018 at 9:37

The computed values are negative only when $A_i = B_j + O(\varepsilon)$, where $\varepsilon$ is the machine precision. The negative values are round-off errors.

If you return the absolute value of your computation, $\left| ||A_i||_2^2 + ||B_j||_2^2 - 2 <A_i,B_j> \right|$, you get the speed-up and positive numbers.

  • $\begingroup$ Regarding non-negativity, see my comment in Anton's answer. It's not just that the results are sometimes negative but they are also sometimes equal up to float32 precision even if the points are distinct. Unfortunately, the problem I'm considering is at too large of a scale to copy the data to higher precision formats. $\endgroup$ Commented Oct 17, 2018 at 17:59
  • $\begingroup$ If your vectors were normalized, you have cases where the result is float32 precision, even if the points are distinct? $\endgroup$
    – user7440
    Commented Oct 17, 2018 at 20:24

Your Answer

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

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