如何用NumPy计算欧几里得距离？
- 2024-12-04 08:56:00
- admin 原创
- 152
问题描述：
我在三维空间中有两个点：
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))
解决方案 1：
使用numpy.linalg.norm：
dist = numpy.linalg.norm(a-b)
这是因为欧氏距离是l2 范数ord，并且中的参数的默认值为numpy.linalg.norm2。有关更多理论，请参阅数据挖掘简介：
解决方案 2：
使用scipy.spatial.distance.euclidean：
from scipy.spatial import distance
a = (1, 2, 3)
b = (4, 5, 6)
dst = distance.euclidean(a, b)
解决方案 3：
对于那些有兴趣同时计算多个距离的人，我已经使用perfplot（我的一个小项目）做了一些比较。
第一个建议是组织数据，使数组具有维度（显然是 C 连续的）。如果在连续的第一个维度上进行添加，则速度会更快，并且使用with 、with或 也(3, n)无关紧要sqrt-sum`axis=0linalg.normaxis=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")
解决方案 4：
我想通过各种性能说明来阐述这个简单的答案。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”进行范围查找，对“norm”进行范围查找，仅仅调用该函数的开销就相当于几十条 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，它是必须调用的第二个函数...)
解决方案 5：
此解题方法的另一个例子：
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)
解决方案 6：
从开始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
解决方案 7：
可以像下面这样完成。我不知道它有多快，但它没有使用 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)))
解决方案 8：
一句漂亮的话：
dist = numpy.linalg.norm(a-b)
但是，如果速度是个问题，我建议在你的机器上进行实验。我发现在我的机器上使用带有平方math运算符的库比一行 NumPy 解决方案要快得多。sqrt`**`
我使用这个简单的程序运行了测试：
#!/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 运行的。
解决方案 9：
我在 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)
解决方案 10：
您可以先减去向量，然后再计算内积。
按照你的例子，
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)
解决方案 11：
我喜欢np.dot（点积）：
a = numpy.array((xa,ya,za))
b = numpy.array((xb,yb,zb))
distance = (np.dot(a-b,a-b))**.5
解决方案 12：
自 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
解决方案 13：
使用 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)))
解决方案 14：
有了它们，a正如b您定义的那样，您还可以使用：
distance = np.sqrt(np.sum((a-b)**2))
解决方案 15：
给定两个点（以 Python 中的列表表示），下面是 Python 中欧几里得距离的一些简洁代码。
def distance(v1,v2):
return sum([(x-y)**2 for (x,y) in zip(v1,v2)])**(0.5)
解决方案 16：
import math
dist = math.hypot(math.hypot(xa-xb, ya-yb), za-zb)
解决方案 17：
计算多维空间的欧几里得距离：
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
解决方案 18：
其他答案适用于浮点数，但不能正确计算容易发生上溢和下溢的整数数据类型的距离。请注意，甚至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 中计算距离，如下所示：
np.linalg.norm(np.maximum(x, y) - np.minimum(x, y))
对于有符号整数类型，您可以先转换为浮点数：
np.linalg.norm(x.astype("float") - y.astype("float"))
具体来说，对于图像数据，可以使用 opencv 的 norm 方法：
import cv2
cv2.norm(x, y, cv2.NORM_L2)
解决方案 19：
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)
解决方案 20：
您可以轻松使用公式
distance = np.sqrt(np.sum(np.square(a-b)))
它实际上只不过是使用勾股定理来计算距离，通过将 Δx、Δy 和 Δz 的平方相加，然后对结果求根。
解决方案 21：
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)))
解决方案 22：
如果您想要更明确的内容，您可以轻松地写出如下公式：
np.sqrt(np.sum((a-b)**2))
即使数组包含 10_000_000 个元素，在我的计算机上仍需 0.1 秒才能运行。
解决方案 23：
1. SciPy 的矢量化cdist()欧几里德距离矩阵
@Nico Schlömer的基准测试表明，scipy 的euclidean()函数比其 numpy 对应函数慢得多。原因是它旨在处理一对点，而不是一个点数组；因此不是矢量化的。此外，他的基准测试使用代码来查找等长数组之间的欧几里得距离。
如果您需要从两个输入集合计算每对点之间的欧几里得距离矩阵，那么还有另一个 SciPy 函数，cdist()它比 numpy 快得多。
考虑以下示例，其中a包含 3 个点和b包含 2 个点。SciPy计算中的每个点到中的每个点cdist()之间的欧几里得距离，因此在此示例中，它将返回一个 3x2 矩阵。a`b`
import numpy as np
from scipy.spatial import distance
a = [(1, 2, 3), (3, 4, 5), (2, 3, 6)]
b = [(1, 2, 3), (4, 5, 6)]
dsts1 = distance.cdist(a, b)
# array([[0. , 5.19615242],
# [3.46410162, 1.73205081],
# [3.31662479, 2.82842712]])
如果我们有一组点，并且想要找到与除自身之外的每个点的最近距离，那么这种方法就特别有用；一种常见的用例是自然语言处理。例如，计算集合中每对点之间的欧几里得距离distance.cdist(a, a)就可以了。由于点到自身的距离为 0，因此该矩阵的对角线将全部为零。
可以使用广播通过 numpy 专用方法执行相同任务。我们只需向其中一个数组添加另一个维度即可。
# using `linalg.norm`
dsts2 = np.linalg.norm(np.array(a)[:, None] - b, axis=-1)
# using a `sqrt` + `sum` + `square`
dsts3 = np.sqrt(np.sum((np.array(a)[:, None] - b)**2, axis=-1))
# equality check
np.allclose(dsts1, dsts2) and np.allclose(dsts1, dsts3) # True
如前所述，cdist()它比 numpy 快得多。以下 perfplot 显示了这一点。1
2. Scikit-learn 的euclidean_distances()
Scikit-learn 是一个相当大的库，因此除非您不将其用于其他用途，否则仅将其导入用于欧几里得距离计算是没有意义的，但为了完整性，它还具有euclidean_distances()和可用于计算欧几里得距离的方法。它还有其他方便的成对距离计算方法值得一试。paired_distances()`pairwise_distances()`
scikit-learn 方法的一个有用之处在于它可以按原样处理稀疏矩阵，而 scipy/numpy 需要将稀疏矩阵转换为数组才能执行计算，因此根据数据的大小，scikit-learn 的方法可能是唯一运行的函数。
举个例子：
from scipy import sparse
from sklearn.metrics import pairwise
A = sparse.random(1_000_000, 3)
b = [(1, 2, 3), (4, 5, 6)]
dsts = pairwise.euclidean_distances(A, b)
1用于生成 perfplot 的代码：
import numpy as np
from scipy.spatial import distance
import perfplot
import matplotlib.pyplot as plt
def sqrt_sum(arr):
return np.sqrt(np.sum((arr[:, None] - arr) ** 2, axis=-1))
def linalg_norm(arr):
return np.linalg.norm(arr[:, None] - arr, axis=-1)
def scipy_cdist(arr):
return distance.cdist(arr, arr)
perfplot.plot(
setup=lambda n: np.random.rand(n, 3),
n_range=[2 ** k for k in range(14)],
kernels=[sqrt_sum, scipy_cdist, linalg_norm],
title="Euclidean distance between arrays of 3-D points",
xlabel="len(x), len(y)",
equality_check=np.allclose
);
解决方案 24：
首先求两个矩阵的差。然后，使用 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
解决方案 25：
使用 NumPy 或一般 Python 来实现这一点的最佳方法是什么？我有：
最好的方法就是最安全，也是最快的
我建议使用 hypot 来获得可靠的结果，因为与编写自己的 sqroot 计算器相比，下溢和溢出的可能性非常小
让我们看看 math.hypot、np.hypot 与 vanillanp.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
参考
解决方案 26：
对于大量距离，我能想到的最快解决方案是使用 numexpr。在我的计算机上，它比使用 numpy einsum 更快：
import numexpr as ne
import numpy as np
a_min_b=a-b
np.sqrt(ne.evaluate("sum((a_min_b)**2,axis=1)"))
