欧几里得距离的高效精确计算

13

在查阅了一些在线资料(12numpyscipyscikitmath),我发现了几种 Python 中计算欧氏距离(Euclidean Distance)的方法:

# 1
numpy.linalg.norm(a-b)

# 2
distance.euclidean(vector1, vector2)

# 3
sklearn.metrics.pairwise.euclidean_distances  

# 4
sqrt((xa-xb)^2 + (ya-yb)^2 + (za-zb)^2)

# 5
dist = [(a - b)**2 for a, b in zip(vector1, vector2)]
dist = math.sqrt(sum(dist))

# 6
math.hypot(x, y)

我想知道哪个方法(上述任何一种或我没有找到的其他方法)在效率和精度方面被认为是最好的。如果有人知道任何讨论这个主题的资源,那就太好了。
我感兴趣的上下文是计算数对之间的欧几里德距离,例如 (52, 106, 35, 12)(33, 153, 75, 10) 之间的距离。

2
不要忘记内置的 math.hypot()。您可以使用 timeit 模块轻松测试速度。 - martineau
1
@martineau 很好的建议,我不知道存在这样的内置方法!(编辑了我的问题以包含它) - user6167676
math.hypot() 的一个可能的注意点是它只能处理二维向量,而你提到的许多其他函数可以处理三维或更高维度的向量。另一方面,如果你只需要处理二维向量,非通用内置函数可能会有速度优势。 - martineau
@martineau 有趣的限制条件,虽然对于我的情况可能是理想的。可能是一个天真的问题:当计算(52, 106, 35, 12)(33, 153, 75, 10)之间的欧几里得距离时,这两个是4D向量吗? - user6167676
不,我的意思是由以下三组端点定义的三个二维向量:(52,33)和(106,153)之间的向量、(106,153)和(35,75)之间的向量,以及(35,75)和(12,10)之间的向量。也许你应该编辑你的问题并展示所需的结果。 - martineau
显示剩余4条评论
6个回答

17

先说结论:

通过使用timeit进行效率测试的结果,我们可以得出关于效率方面的结论:

Method5 (zip, math.sqrt) > Method1 (numpy.linalg.norm) > Method2 (scipy.spatial.distance) > Method3 (sklearn.metrics.pairwise.euclidean_distances )

对于Method4,由于不适用于一般情况并且通常等同于Method5,因此我并没有真正测试它。

令人惊讶的是,Method5实际上是最快的。而使用numpyMethod1,由于它在C语言中得到了很好的优化,所以像我们所预期的那样成为了第二快的方法。

对于scipy.spatial.distance,如果你直接查看函数定义,你会发现它实际上是使用了numpy.linalg.norm,只不过它会在实际运行numpy.linalg.norm之前对两个输入向量进行验证。这就是它略慢于numpy.linalg.norm的原因。

最后对于sklearn,根据文档:

  

与其他计算距离的方法相比,这种公式有两个优点。首先,在处理稀疏数据时具有计算效率。其次,如果一个参数变化而另一个参数保持不变,则可以预先计算dot(x, x)和/或dot(y,y)。   然而,这并不是计算距离的最精确方法,并且该函数返回的距离矩阵可能不会完全对称,因为需要。

由于在你的问题中,你想使用一组固定的数据,因此该实现的优势没有得到体现。由于性能和精度之间的权衡,它也给出了所有方法中最差的精度。

关于精度方面Method5=Metho1=Method2>Method3

效率测试脚本:

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

# 1
def eudis1(v1, v2):
    return np.linalg.norm(v1-v2)

# 2
def eudis2(v1, v2):
    return distance.euclidean(v1, v2)

# 3
def eudis3(v1, v2):
    return euclidean_distances(v1, v2)

# 5
def eudis5(v1, v2):
    dist = [(a - b)**2 for a, b in zip(v1, v2)]
    dist = math.sqrt(sum(dist))
    return dist

dis1 = (52, 106, 35, 12)
dis2 = (33, 153, 75, 10)
v1, v2 = np.array(dis1), np.array(dis2)

import timeit

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

wrappered1 = wrapper(eudis1, v1, v2)
wrappered2 = wrapper(eudis2, v1, v2)
wrappered3 = wrapper(eudis3, v1, v2)
wrappered5 = wrapper(eudis5, v1, v2)
t1 = timeit.repeat(wrappered1, repeat=3, number=100000)
t2 = timeit.repeat(wrappered2, repeat=3, number=100000)
t3 = timeit.repeat(wrappered3, repeat=3, number=100000)
t5 = timeit.repeat(wrappered5, repeat=3, number=100000)

print('\n')
print('t1: ', sum(t1)/len(t1))
print('t2: ', sum(t2)/len(t2))
print('t3: ', sum(t3)/len(t3))
print('t5: ', sum(t5)/len(t5))

效率测试输出:

t1:  0.654838958307
t2:  1.53977598714
t3:  6.7898791732
t5:  0.422228400305

精准度测试脚本 & 结果:

In [8]: eudis1(v1,v2)
Out[8]: 64.60650122085238

In [9]: eudis2(v1,v2)
Out[9]: 64.60650122085238

In [10]: eudis3(v1,v2)
Out[10]: array([[ 64.60650122]])

In [11]: eudis5(v1,v2)
Out[11]: 64.60650122085238

1
请添加内置的 math.hypot()。(顺便提一下,OP正在使用Python 3)。 - martineau
@MaThMaX 好东西!正如 @martineau 所建议的那样,如果你能添加内置的 math.hypot(),那将是令人惊叹的。特别是因为我以前从未使用/听说过它。 - user6167676
当计算小向量之间的距离时,性能效率为Method5(zip,math.sqrt)> Method1(numpy.linalg.norm)。然而,当我测试大于128的向量大小时,Method1 > Method5。 - RyanLiu
1
关于sklearn和文档:计算优势只在更大的距离矩阵中显示出来。基准测试实际上测试了单个距离值,但是如果您有数千个点并且想要计算它们之间的成对距离并将结果存储在矩阵中怎么办?这就是sklearn变得更加优越(以损失精度为代价)的场景-正如文档中所示。 - no_use123

6

这并不完全回答问题,但可能值得提到的是,如果你对欧几里得距离本身不感兴趣,而只想比较欧几里得距离之间的大小关系,那么平方根是单调函数,也就是说,当且仅当 x < y 时,x**(1/2) < y**(1/2)。

因此,如果你不需要显式距离,而仅仅想知道 vector1 的欧几里得距离是否更接近一个向量列表(称为 vectorlist),你可以避免计算昂贵(精度和时间)的平方根,而用类似下面的方法:

min(vectorlist, key = lambda compare: sum([(a - b)**2 for a, b in zip(vector1, compare)])


1

在接受的答案上改进基准测试后,我发现,假设您已经以numpy数组格式获得输入,则可以更好地编写method5:

import numpy as np
from numba import jit

@jit(nopython=True)
def euclidian_distance(y1, y2):
    return np.sqrt(np.sum((y1-y2)**2)) # based on pythagorean

速度测试:

euclidian_distance(y1, y2)
# 2.03 µs ± 138 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

np.linalg.norm(y1-y2)
# 17.6 µs ± 5.08 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

有趣的事实是,你可以在numpy函数中添加jit

@jit(nopython=True)
def jit_linalg(y1, y2):
    return np.linalg.norm(y1-y2)

jit_linalg(y[i],y[j])
# 2.91 µs ± 261 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

1
这里有一个关于如何仅使用numpy的示例。
import numpy as np

a = np.array([3, 0])
b = np.array([0, 4])

c = np.sqrt(np.sum(((a - b) ** 2)))
# c == 5.0

0

通常而言,尽可能使用scipynumpy的实现,因为它们是向量化的,比本地Python代码快得多。(主要原因是:C语言实现,向量化消除了循环的类型检查开销。)

(另外:我的答案没有涉及精确度,但我认为相同的原则也适用于效率方面。)

作为额外的奖励,我会提供一些关于如何分析您的代码以衡量效率的信息。如果您正在使用IPython解释器,秘密就在于使用%prun命令。

In [1]: import numpy

In [2]: from scipy.spatial import distance

In [3]: c1 = numpy.array((52, 106, 35, 12))

In [4]: c2 = numpy.array((33, 153, 75, 10))

In [5]: %prun distance.euclidean(c1, c2)
         35 function calls in 0.000 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 linalg.py:1976(norm)
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.dot}
        6    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.array}
        4    0.000    0.000    0.000    0.000 numeric.py:406(asarray)
        1    0.000    0.000    0.000    0.000 distance.py:232(euclidean)
        2    0.000    0.000    0.000    0.000 distance.py:152(_validate_vector)
        2    0.000    0.000    0.000    0.000 shape_base.py:9(atleast_1d)
        1    0.000    0.000    0.000    0.000 misc.py:11(norm)
        1    0.000    0.000    0.000    0.000 function_base.py:605(asarray_chkfinite)
        2    0.000    0.000    0.000    0.000 numeric.py:476(asanyarray)
        1    0.000    0.000    0.000    0.000 {method 'ravel' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 linalg.py:111(isComplexType)
        1    0.000    0.000    0.000    0.000 <string>:1(<module>)
        2    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.issubclass}
        4    0.000    0.000    0.000    0.000 {built-in method builtins.len}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        2    0.000    0.000    0.000    0.000 {method 'squeeze' of 'numpy.ndarray' objects}


In [6]: %prun numpy.linalg.norm(c1 - c2)
         10 function calls in 0.000 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 linalg.py:1976(norm)
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.dot}
        1    0.000    0.000    0.000    0.000 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 numeric.py:406(asarray)
        1    0.000    0.000    0.000    0.000 {method 'ravel' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 linalg.py:111(isComplexType)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.issubclass}
        1    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.array}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

%prun 的作用是告诉你一个函数调用需要多长时间才能运行,包括一些跟踪信息以找出瓶颈所在。在这种情况下,scipy.spatial.distance.euclideannumpy.linalg.norm 实现都非常快。假设你定义了一个函数 dist(vect1, vect2),你可以使用相同的 IPython 魔法调用进行分析。另外,%prun 也适用于 Jupyter 笔记本,并且你可以通过将 %%prun 放在代码单元格的第一行来分析整个代码单元格,而不仅仅是一个函数。


这个答案是错误的。对于特定的用例,如果您不需要任何处理程序,使用基本的numpy运算符有时比使用numpy函数更快。此外,还有numba可以将您的函数编译成jit。尝试我能找到的最快的方法:def euclidian_distance(y1, y2): return np.sqrt(np.sum((y1-y2)**2))``` - Muhammad Yasirroni

0

我不知道精度和速度如何与您提到的其他库相比,但是您可以使用内置的math.hypot()函数来处理二维向量:

from math import hypot

def pairwise(iterable):
    "s -> (s0, s1), (s1, s2), (s2, s3), ..."
    a, b = iter(iterable), iter(iterable)
    next(b, None)
    return zip(a, b)

a = (52, 106, 35, 12)
b = (33, 153, 75, 10)

dist = [hypot(p2[0]-p1[0], p2[1]-p1[1]) for p1, p2 in pairwise(tuple(zip(a, b)))]
print(dist)  # -> [131.59027319676787, 105.47511554864494, 68.94925670375281]

谢谢,我会尝试测试和计时。您能简要解释一下pairwise方法是做什么的吗? - user6167676
1
当然。pairwise()函数是在itertools recipes文档中展示的那个函数的略微变化。它和原来的函数一样,从传入的可迭代参数中按照函数开头的doc string所示的顺序返回值对。 - martineau

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接