如何使用NumPy计算欧几里德距离？

我在3D空间中有两点：

a = (ax, ay, az)
b = (bx, by, bz)

我想计算它们之间的距离：

dist = sqrt((ax-bx)^2 + (ay-by)^2 + (az-bz)^2)

我如何用NumPy做到这一点？我有：

import numpy
a = numpy.array((ax, ay, az))
b = numpy.array((bx, by, bz))
1275063 次浏览
小开

这种解题方法的另一个实例：

def dist(x,y):
return numpy.sqrt(numpy.sum((x-y)**2))


a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))
dist_a_b = dist(a,b)
小开
最佳答案

使用numpy.linalg.norm：

dist = numpy.linalg.norm(a-b)

这之所以有效，是因为欧几里得距离是l2范数，而numpy.linalg.norm中ord参数的默认值是2。 有关更多理论，请参阅数据挖掘简介:

在此处输入图片描述

小开

一个漂亮的单行：

dist = numpy.linalg.norm(a-b)

但是，如果速度是一个问题，我建议在您的机器上进行实验。我发现在我的机器上使用math库的sqrt和**运算符作为正方形要比单行NumPy解决方案快得多。

我使用这个简单的程序运行我的测试：

#!/usr/bin/python
import math
import numpy
from random import uniform


def fastest_calc_dist(p1,p2):
return math.sqrt((p2[0] - p1[0]) ** 2 +
(p2[1] - p1[1]) ** 2 +
(p2[2] - p1[2]) ** 2)


def math_calc_dist(p1,p2):
return math.sqrt(math.pow((p2[0] - p1[0]), 2) +
math.pow((p2[1] - p1[1]), 2) +
math.pow((p2[2] - p1[2]), 2))


def numpy_calc_dist(p1,p2):
return numpy.linalg.norm(numpy.array(p1)-numpy.array(p2))


TOTAL_LOCATIONS = 1000


p1 = dict()
p2 = dict()
for i in range(0, TOTAL_LOCATIONS):
p1[i] = (uniform(0,1000),uniform(0,1000),uniform(0,1000))
p2[i] = (uniform(0,1000),uniform(0,1000),uniform(0,1000))


total_dist = 0
for i in range(0, TOTAL_LOCATIONS):
for j in range(0, TOTAL_LOCATIONS):
dist = fastest_calc_dist(p1[i], p2[j]) #change this line for testing
total_dist += dist


print total_dist

在我的机器上，math_calc_dist比numpy_calc_dist快得多：1.5秒对23.5秒。

为了获得fastest_calc_dist和math_calc_dist之间的可测量差异，我必须将TOTAL_LOCATIONS提高到6000。然后fastest_calc_dist需要约50秒，而math_calc_dist需要约60秒。

您还可以尝试numpy.sqrt和numpy.square，尽管它们都比我机器上的math替代品慢。

我的测试使用Python 2.6.6运行。

小开

你可以只减去向量，然后内积。

以你为榜样，

a = numpy.array((xa, ya, za))
b = numpy.array((xb, yb, zb))


tmp = a - b
sum_squared = numpy.dot(tmp.T, tmp)
result = numpy.sqrt(sum_squared)
小开

它可以像下面这样完成。我不知道它有多快，但它没有使用NumPy。

from math import sqrt
a = (1, 2, 3) # Data point 1
b = (4, 5, 6) # Data point 2
print sqrt(sum( (a - b)**2 for a, b in zip(a, b)))
小开

我在matplotlib.mlab中找到了一个“dist”函数，但我认为它不够方便。

我把它贴在这里只是为了参考。

import numpy as np
import matplotlib as plt


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


# Distance between a and b
dis = plt.mlab.dist(a, b)
小开

使用scipy.spatial.distance.euclidean：

from scipy.spatial import distance
a = (1, 2, 3)
b = (4, 5, 6)
dst = distance.euclidean(a, b)
小开

下面是Python中欧几里德距离的一些简明代码，给出了Python中表示为列表的两点。

def distance(v1,v2):
return sum([(x-y)**2 for (x,y) in zip(v1,v2)])**(0.5)
小开

我喜欢np.dot（点积）： 

a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))


distance = (np.dot(a-b,a-b))**.5
小开

定义了a和b，您还可以使用：

distance = np.sqrt(np.sum((a-b)**2))
小开

计算多维空间的欧几里德距离：

 import math


x = [1, 2, 6]
y = [-2, 3, 2]


dist = math.sqrt(sum([(xi-yi)**2 for xi,yi in zip(x, y)]))
5.0990195135927845
小开

我想用各种性能笔记来阐述简单的答案。np.linalg.norm可能会做的比你需要的更多：

dist = numpy.linalg.norm(a-b)

首先-此函数旨在处理列表并返回所有值，例如比较从pA到点集sP的距离：

sP = set(points)
pA = point
distances = np.linalg.norm(sP - pA, ord=2, axis=1.)  # 'distances' is a list

记住几件事：

  • Python函数调用是昂贵的。
  • [常规]Python不缓存名称查找。

所以

def distance(pointA, pointB):
dist = np.linalg.norm(pointA - pointB)
return dist

并不像看起来那么无辜。

>>> dis.dis(distance)
2           0 LOAD_GLOBAL              0 (np)
2 LOAD_ATTR                1 (linalg)
4 LOAD_ATTR                2 (norm)
6 LOAD_FAST                0 (pointA)
8 LOAD_FAST                1 (pointB)
10 BINARY_SUBTRACT
12 CALL_FUNCTION            1
14 STORE_FAST               2 (dist)


3          16 LOAD_FAST                2 (dist)
18 RETURN_VALUE

首先-每次我们调用它时，我们都必须对“np”进行全局查找，对“linalg”进行范围查找，对“规范”进行范围查找，并且仅呼吁的开销函数可以相当于数十条python指令。

最后，我们浪费了两个操作来存储结果并重新加载它以返回…

改进的第一步：加快查找速度，跳过存储

def distance(pointA, pointB, _norm=np.linalg.norm):
return _norm(pointA - pointB)

我们得到了更精简的：

>>> dis.dis(distance)
2           0 LOAD_FAST                2 (_norm)
2 LOAD_FAST                0 (pointA)
4 LOAD_FAST                1 (pointB)
6 BINARY_SUBTRACT
8 CALL_FUNCTION            1
10 RETURN_VALUE

不过，函数调用开销仍然需要一些工作。你需要做基准测试来确定你自己是否可以更好地进行数学运算：

def distance(pointA, pointB):
return (
((pointA.x - pointB.x) ** 2) +
((pointA.y - pointB.y) ** 2) +
((pointA.z - pointB.z) ** 2)
) ** 0.5  # fast sqrt

在某些平台上，**0.5比math.sqrt更快。您的里程可能会有所不同。

****高级性能说明。

你为什么要计算距离？如果唯一的目的是显示它，

 print("The target is %.2fm away" % (distance(a, b)))

但是，如果您正在比较距离，进行范围检查等，我想添加一些有用的性能观察结果。

让我们来看两种情况：按距离排序或剔除列表中满足范围约束的项。

# Ultra naive implementations. Hold onto your hat.


def sort_things_by_distance(origin, things):
return things.sort(key=lambda thing: distance(origin, thing))


def in_range(origin, range, things):
things_in_range = []
for thing in things:
if distance(origin, thing) <= range:
things_in_range.append(thing)

我们需要记住的第一件事是，我们使用毕达哥拉斯来计算距离（dist = sqrt(x^2 + y^2 + z^2)），所以我们进行了很多sqrt调用。数学101：

dist = root ( x^2 + y^2 + z^2 )
:.
dist^2 = x^2 + y^2 + z^2
and
sq(N) < sq(M) iff M > N
and
sq(N) > sq(M) iff N > M
and
sq(N) = sq(M) iff N == M

简而言之：直到我们实际需要以X为单位的距离而不是X^2，我们可以消除计算中最难的部分。

# Still naive, but much faster.


def distance_sq(left, right):
""" Returns the square of the distance between left and right. """
return (
((left.x - right.x) ** 2) +
((left.y - right.y) ** 2) +
((left.z - right.z) ** 2)
)


def sort_things_by_distance(origin, things):
return things.sort(key=lambda thing: distance_sq(origin, thing))


def in_range(origin, range, things):
things_in_range = []


# Remember that sqrt(N)**2 == N, so if we square
# range, we don't need to root the distances.
range_sq = range**2


for thing in things:
if distance_sq(origin, thing) <= range_sq:
things_in_range.append(thing)

太好了，这两个函数都不再做任何昂贵的平方根。这样会快得多，但在你进一步讨论之前，先检查一下自己：为什么sort_things_by_distance两次都需要一个“天真”的免责声明？答案在最底部（*a1）。

我们可以通过将其转换为生成器来改进in_range：

def in_range(origin, range, things):
range_sq = range**2
yield from (thing for thing in things
if distance_sq(origin, thing) <= range_sq)

如果你正在做这样的事情，这尤其有好处：

if any(in_range(origin, max_dist, things)):
...

但是如果你接下来要做的事情需要一段距离，

for nearby in in_range(origin, walking_distance, hotdog_stands):
print("%s %.2fm" % (nearby.name, distance(origin, nearby)))

考虑产生元组：

def in_range_with_dist_sq(origin, range, things):
range_sq = range**2
for thing in things:
dist_sq = distance_sq(origin, thing)
if dist_sq <= range_sq: yield (thing, dist_sq)

如果您可能会进行范围检查（“找到靠近X且在Y的Nm范围内的东西”，这可能特别有用，因为您不必再次计算距离）。

但是，如果我们正在搜索一个非常大的things列表，并且我们预计其中很多不值得考虑呢？

实际上有一个非常简单的优化：

def in_range_all_the_things(origin, range, things):
range_sq = range**2
for thing in things:
dist_sq = (origin.x - thing.x) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.y - thing.y) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.z - thing.z) ** 2
if dist_sq <= range_sq:
yield thing

这是否有用将取决于“事物”的大小。

def in_range_all_the_things(origin, range, things):
range_sq = range**2
if len(things) >= 4096:
for thing in things:
dist_sq = (origin.x - thing.x) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.y - thing.y) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.z - thing.z) ** 2
if dist_sq <= range_sq:
yield thing
elif len(things) > 32:
for things in things:
dist_sq = (origin.x - thing.x) ** 2
if dist_sq <= range_sq:
dist_sq += (origin.y - thing.y) ** 2 + (origin.z - thing.z) ** 2
if dist_sq <= range_sq:
yield thing
else:
... just calculate distance and range-check it ...

再一次，考虑屈服dist_sq。我们的热狗示例变成：

# Chaining generators
info = in_range_with_dist_sq(origin, walking_distance, hotdog_stands)
info = (stand, dist_sq**0.5 for stand, dist_sq in info)
for stand, dist in info:
print("%s %.2fm" % (stand, dist))

（*a1：sort_things_by_distance的排序键对每个项目都distance_sq，而那个看起来无辜的键是一个lambda，这是必须调用的第二个函数…）

小开

对于任何对一次计算多个距离感兴趣的人，我使用灌流图（我的一个小项目）做了一些比较。

第一个建议是组织您的数据，使数组具有维度(3, n)（并且显然是C-连续的）。如果添加发生在连续的第一维，速度会更快，如果您将sqrt-sum与axis=0、linalg.norm与axis=0或

a_min_b = a - b
numpy.sqrt(numpy.einsum('ij,ij->j', a_min_b, a_min_b))

这是，略微差一点，最快的变体。（这实际上也适用于一行。）

您在第二个轴上求和的变体axis=1都要慢得多。

在此处输入图片描述

重现情节的代码：

import numpy
import perfplot
from scipy.spatial import distance




def linalg_norm(data):
a, b = data[0]
return numpy.linalg.norm(a - b, axis=1)




def linalg_norm_T(data):
a, b = data[1]
return numpy.linalg.norm(a - b, axis=0)




def sqrt_sum(data):
a, b = data[0]
return numpy.sqrt(numpy.sum((a - b) ** 2, axis=1))




def sqrt_sum_T(data):
a, b = data[1]
return numpy.sqrt(numpy.sum((a - b) ** 2, axis=0))




def scipy_distance(data):
a, b = data[0]
return list(map(distance.euclidean, a, b))




def sqrt_einsum(data):
a, b = data[0]
a_min_b = a - b
return numpy.sqrt(numpy.einsum("ij,ij->i", a_min_b, a_min_b))




def sqrt_einsum_T(data):
a, b = data[1]
a_min_b = a - b
return numpy.sqrt(numpy.einsum("ij,ij->j", a_min_b, a_min_b))




def setup(n):
a = numpy.random.rand(n, 3)
b = numpy.random.rand(n, 3)
out0 = numpy.array([a, b])
out1 = numpy.array([a.T, b.T])
return out0, out1




b = perfplot.bench(
setup=setup,
n_range=[2 ** k for k in range(22)],
kernels=[
linalg_norm,
linalg_norm_T,
scipy_distance,
sqrt_sum,
sqrt_sum_T,
sqrt_einsum,
sqrt_einsum_T,
],
xlabel="len(x), len(y)",
)
b.save("norm.png")
小开 
import numpy as np
from scipy.spatial import distance
input_arr = np.array([[0,3,0],[2,0,0],[0,1,3],[0,1,2],[-1,0,1],[1,1,1]])
test_case = np.array([0,0,0])
dst=[]
for i in range(0,6):
temp = distance.euclidean(test_case,input_arr[i])
dst.append(temp)
print(dst)
小开 
import math


dist = math.hypot(math.hypot(xa-xb, ya-yb), za-zb)
小开

你可以很容易地使用公式

distance = np.sqrt(np.sum(np.square(a-b)))

它实际上只是使用毕达哥拉斯定理来计算距离，通过添加Δx，Δy和Δz的平方并生根结果。

小开

首先找到两个矩阵的差值。然后，使用numpy的乘法命令进行元素乘法。然后，找到元素乘以新矩阵的求和。最后，找到求和的平方根。

def findEuclideanDistance(a, b):
euclidean_distance = a - b
euclidean_distance = np.sum(np.multiply(euclidean_distance, euclidean_distance))
euclidean_distance = np.sqrt(euclidean_distance)
return euclidean_distance
小开

从Python 3.8开始，math模块直接提供dist函数，它返回两点之间的欧几里得距离（以元组或坐标列表给出）：

from math import dist


dist((1, 2, 6), (-2, 3, 2)) # 5.0990195135927845

如果您正在使用列表：

dist([1, 2, 6], [-2, 3, 2]) # 5.0990195135927845
小开

从Python 3.8开始

从Python 3.8开始，math模块包含函数math.dist()。
看这里https://docs.python.org/3.8/library/math.html#math.dist。

math.dist（p1，p2）
返回两点p1和p2之间的欧几里得距离， 每个都作为坐标序列（或可迭代）给出。

import math
print( math.dist( (0,0),   (1,1)   )) # sqrt(2) -> 1.4142
print( math.dist( (0,0,0), (1,1,1) )) # sqrt(3) -> 1.7321
小开

使用Python 3.8，这很容易。

https://docs.python.org/3/library/math.html#math.dist

math.dist(p, q)

返回两点p和q之间的欧几里得距离，每个点都给定 作为坐标序列（或可迭代）。这两个点必须有 相同的维度。

大致相当于：

sqrt(sum((px - qx) ** 2.0 for px, qx in zip(p, q)))

小开 
import numpy as np
# any two python array as two points
a = [0, 0]
b = [3, 4]

您首先将列表更改为numpy数组并这样做：print(np.linalg.norm(np.array(a) - np.array(b)))。直接从python列表中的第二个方法为：print(np.linalg.norm(np.subtract(a,b)))

小开

其他答案适用于浮点数，但不能正确计算容易溢出和下限溢位的整数dtype的距离。请注意，即使scipy.distance.euclidean也有这个问题：

>>> a1 = np.array([1], dtype='uint8')
>>> a2 = np.array([2], dtype='uint8')
>>> a1 - a2
array([255], dtype=uint8)
>>> np.linalg.norm(a1 - a2)
255.0
>>> from scipy.spatial import distance
>>> distance.euclidean(a1, a2)
255.0

这很常见，因为许多图像库将图像表示为dtype="uint8"的ndarray。这意味着如果你有一张灰度图像由非常暗的灰色像素组成（假设所有像素都有颜色#000001），并且你将其与黑色图像（#000000）进行了差异，你最终可以在所有单元格中得到x-y由255组成，这注册为这两张图像彼此相距甚远。对于无符号整数类型（例如uint8），你可以安全地以numpy as计算距离：

np.linalg.norm(np.maximum(x, y) - np.minimum(x, y))

对于有符号整数类型，您可以首先强制转换为浮点数：

np.linalg.norm(x.astype("float") - y.astype("float"))

对于特定的图像数据，您可以使用opencv的规范方法：

import cv2
cv2.norm(x, y, cv2.NORM_L2)
小开

使用NumPy或一般的Python执行此操作的最佳方法是什么？我有：

最好的方法是最安全也是最快的

我建议使用hypot来获得可靠的结果，因为与编写自己的sqroot计算器相比，下限溢位和溢出的机会很小

让我们看看math.hypot，np.hypotvs香草np.sqrt(np.sum((np.array([i, j, k])) ** 2, axis=1))

i, j, k = 1e+200, 1e+200, 1e+200
math.hypot(i, j, k)
# 1.7320508075688773e+200

np.sqrt(np.sum((np.array([i, j, k])) ** 2))
# RuntimeWarning: overflow encountered in square

速度明智math.hypot看起来更好

%%timeit
math.hypot(i, j, k)
# 100 ns ± 1.05 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

%%timeit
np.sqrt(np.sum((np.array([i, j, k])) ** 2))
# 6.41 µs ± 33.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

下溢

i, j = 1e-200, 1e-200
np.sqrt(i**2+j**2)
# 0.0

溢出

i, j = 1e+200, 1e+200
np.sqrt(i**2+j**2)
# inf

没有暗流

i, j = 1e-200, 1e-200
np.hypot(i, j)
# 1.414213562373095e-200

无溢出

i, j = 1e+200, 1e+200
np.hypot(i, j)
# 1.414213562373095e+200

参考

小开

对于大量距离，我能想到的最快的解决方案是使用numexr。在我的机器上，它比使用numpy einsum更快：

import numexpr as ne
import numpy as np
np.sqrt(ne.evaluate("sum((a_min_b)**2,axis=1)"))

提交答案