计算numpy数组之间的距离

Question

我有两个数组，第一个np.array是源点，第二个np.array是我需要计算的所有距离。

例如：

 import numpy as np
 from_array = np.array([(0,1), (1,1), ..., (x,y)])
 to_array = np.array([(5,1), (3,1), ..., (x,y)])

我需要做的是取from_array的第一个条目，并计算from_array [0]到to_array中所有点之间的距离，然后保留最大距离。

所以我可以用暴力方法解决：

 def get_distances(from_array, to_array):
     results = []
     distances = []
     for pt in from_array:
        for to in to_array:
            results.append(calc_dist(pt, to))
        distances.append(results)
     return distances

但是这种方法很慢，我正在寻找一种优化的计算方式，因为我可能会有成千上万个点需要计算。

最终目标是计算豪斯多夫距离。

fhd = np.mean(np.min(SomeDistanceArray,axis=0))
rhd = np.mean(np.min(SomeDistanceArray,axis=1))
print (max(fhd, rhd))

我希望仅使用numpy完成这个任务。我的距离可以是欧几里得距离或平方欧几里得距离。

因此，我正在寻求一种优化的方法来计算两个np.arrays的欧几里得距离方法。需要注意的是，数组1可能比数组2有更多的行。这意味着2D数组（x，y）的长度可能会将10行与30行进行比较。

- code base 5000

---

看看Scipy的cdist。 - Divakar

我感觉就像昨天已经回答了这个问题。 - Daniel F

Scipy 版本 0.19.0（截至2017年1月的开发版本）具有有向Hausdorff距离的实现，您可能需要了解一下。它使用高效算法（参见论文），其运行时间接近O(m)，而不是O(m*n)。 - weiji14

Answer 1

这是一个基于NumPy的方法，使用np.einsum -

subs = from_array[:,None] - to_array
sq_eucliean_dist = np.einsum('ijk,ijk->ij',subs,subs)
eucliean_dist = np.sqrt(sq_eucliean_dist)

注意: 如果您接下来要计算 np.mean(np.min(SomeDistanceArray,axis=0))，则可以跳过对 eucliean_dist 的计算，直接使用 sq_eucliean_dist 作为 SomeDistanceArray，因为计算平方根会非常昂贵。

---

np.einsum('ijk, ijk->ij', subs, subs)是什么意思？它对同一个数组subs执行逐元素乘法，即平方，然后沿着最后一个轴进行求和约简，因此在该约简过程中丢失它。那么，为什么不明确地进行平方和求和呢？好处是np.einsum可以一步完成平方和求和，从而提高性能效率。因此，如果from_array是(N x 2)数组，to_array是(M x 2)数组，则np.einsum的输出将是作为二维数组(N x M)的平方欧几里得距离。关于字符串符号本身的更多信息将涉及更长的讨论，其中一些可以在this post和早期发布的官方文档链接中找到。

- Divakar

---

你能解释一下einsum函数吗？我不熟悉ijk、ikj这种符号。>ij - code base 5000

@josh1234 在上面添加了一些注释。 - Divakar

Answer 2

仅使用Numpy，因此不使用scipy.spatial.distance.cdist 首先，不要使用元组，而要使用2xN和2xM数组。然后进行广播。

np.linalg.norm(from_array[:,:,None]-to_array[:,None,:], axis=0)

如果您使用的是旧版本的numpy，没有可向量化的linalg.norm（即您正在使用Abaqus），请执行以下操作：

np.sum((from_array[:,:,None]-to_array[:,None,:])**2, axis=0).__pow__(0.5)

- Daniel F

---

是的，只需确保公共维度在前面，比如2x30和2x50，这样你就可以像上面那样广播。 - Daniel F

抱歉，我的意思是，我需要旋转我的数组？ - code base 5000

转置可以通过获取 T 属性来完成，例如 x.T 是 x 的转置。 - Mad Physicist

为什么cdist不起作用？这只是问题的限制吗？ - 3pitt

@DanielF 无法相信这不是被接受的答案。 - DuttaA