在Python中为NumPy计算欧几里得距离

Question

我刚开始学习Python，所以这个问题可能看起来很琐碎。但是，我没有找到一个类似于我的情况。我有一个由20个节点的坐标矩阵。我想要计算这组节点中所有节点对之间的欧几里得距离，并将它们存储在一个成对矩阵中。例如，如果我有20个节点，我想要最终的结果是一个(20,20)大小的矩阵，其中包含每对节点之间的欧几里得距离值。我尝试使用for循环遍历坐标集合中的每个元素，并计算欧几里得距离，代码如下：

ncoord=numpy.matrix('3225   318;2387    989;1228    2335;57      1569;2288  8138;3514   2350;7936   314;9888    4683;6901   1834;7515   8231;709   3701;1321    8881;2290   2350;5687   5034;760    9868;2378   7521;9025   5385;4819   5943;2917   9418;3928   9770')
n=20 
c=numpy.zeros((n,n))
for i in range(0,n):
    for j in range(i+1,n):
        c[i][j]=math.sqrt((ncoord[i][0]-ncoord[j][0])**2+(ncoord[i][1]-ncoord[j][1])**2)

然而，我遇到了一个“input must be a square array”的错误。我想知道这里发生了什么。谢谢。

- Fairy

---

请[编辑]您的问题，包括ncoord的定义。感谢您提高问题的参考价值并使其更易回答！ - Nathan Tuggy

你的n是多少？for j in range(i+1,n-1)将会执行j=i+1, i+2, ..., n-2。我猜你想让这两个范围都增加到n，而不是n-1。 - MarkG

@MarkG 是的，我有20个节点（n=20），我希望两个索引都可以到达n。我尝试使用n而不是n-1，但是我得到了相同的错误。我可以很容易地在MATLAB中编写此代码，但我必须使用Python。Python中的索引方式不同，所以我可能错了。 - Fairy

那么你的两个for循环都应该到n：for i in range(0,n):和for j in range(i+1,n):。如果这不是你的错误，那么你需要展示更多的代码。 - MarkG

@MarkG 是的，这不是我的错误。我的代码就是我在主要问题中提到的那样。我没有更多的东西。 - Fairy

你为什么选择使用矩阵而不是数组？你考虑过将输入数据解析成np.array([ [x0,y0],[x1,y1] ... [xn,yn] ],dtype=([('x','<f8'),('y','<f8')])的形式吗？请参考此链接：https://dev59.com/I3M_5IYBdhLWcg3wXyHZ?rq=1 - user1121588

Answer 1

针对这种情况，有比嵌套for循环更快的替代方法。我将向您展示两种不同的方法 - 第一种是更通用的方法，将介绍广播和矢量化的概念；第二种使用了一个更方便的scipy库函数。

---

一般的方法是使用广播和向量化。

---

我建议你首先做的一件事是改用np.array而不是np.matrix。数组更受欢迎，有许多原因，最重要的是它们可以有>2个维度，并且它们使逐元素乘法变得不那么笨拙。

import numpy as np

ncoord = np.array(ncoord)

使用数组，我们可以通过插入新的单例维度并广播减法来消除嵌套的for循环：

# indexing with None (or np.newaxis) inserts a new dimension of size 1
print(ncoord[:, :, None].shape)
# (20, 2, 1)

# by making the 'inner' dimensions equal to 1, i.e. (20, 2, 1) - (1, 2, 20),
# the subtraction is 'broadcast' over every pair of rows in ncoord
xydiff = ncoord[:, :, None] - ncoord[:, :, None].T

print(xydiff.shape)
# (20, 2, 20)

这相当于使用嵌套的for循环遍历每一对行，但速度要快得多！

xydiff2 = np.zeros((20, 2, 20), dtype=xydiff.dtype)
for ii in range(20):
    for jj in range(20):
        for kk in range(2):
            xydiff[ii, kk, jj] = ncoords[ii, kk] - ncoords[jj, kk]

# check that these give the same result
print(np.all(xydiff == xydiff2))
# True

我们可以使用向量化操作来完成剩下的工作:

# we square the differences and sum over the 'middle' axis, equivalent to
# computing (x_i - x_j) ** 2 + (y_i - y_j) ** 2
ssdiff = (xydiff * xydiff).sum(1)

# finally we take the square root
D = np.sqrt(ssdiff)

整个事情可以像这样一行完成：

D = np.sqrt(((ncoord[:, :, None] - ncoord[:, :, None].T) ** 2).sum(1))

---

2. 使用 pdist 的懒惰方式

---

原来已经有一个快速便捷的函数可以计算所有成对距离：scipy.spatial.distance.pdist。

from scipy.spatial.distance import pdist, squareform

d = pdist(ncoord)

# pdist just returns the upper triangle of the pairwise distance matrix. to get
# the whole (20, 20) array we can use squareform:

print(d.shape)
# (190,)

D2 = squareform(d)
print(D2.shape)
# (20, 20)

# check that the two methods are equivalent
print np.all(D == D2)
# True

- ali_m

---

这个广播对我来说很神奇。我该如何获得一些关于它的直觉？ - sakimarquis

感谢这个神奇的方法，但它仍然比叉积慢得多，尽管复杂度看起来是一样的。 - Jean Paul

当我使用方法1（非scipy）计算一个大矩阵（1000 * 20000）时，我也遇到了一些内存问题，但是使用方法2（scipy）却没有这个问题。 - Jean Paul

Answer 2

for i in range(0, n):
    for j in range(i+1, n):
        c[i, j] = math.sqrt((ncoord[i, 0] - ncoord[j, 0])**2 
        + (ncoord[i, 1] - ncoord[j, 1])**2)

注意: 对于Numpy矩阵，ncoord[i, j]不同于ncoord[i][j]。这是混淆的根源。如果ncoord是一个Numpy数组，则它们将返回相同的结果。

对于Numpy矩阵，ncoord[i]返回ncoord的第i行，它本身是一个形状为1 x 2的Numpy矩阵对象。因此，ncoord[i][j]实际上是指：获取ncoord的第i行，然后获取该1 x 2矩阵的j列。当 j > 0时，这就是索引问题所在。

关于您对c[i][j]赋值“work”的评论，它不应该。至少在我的Numpy 1.9.1中，如果您的索引i和j迭代到n，那么它不应该工作。

另外，请记得将矩阵c的转置添加到本身。

建议使用Numpy数组而不是矩阵。请参阅此帖子。

如果您的坐标存储为Numpy数组，则可以计算成对距离：

from scipy.spatial.distance import pdist

pairwise_distances = pdist(ncoord, metric="euclidean", p=2)

或者简单地说

pairwise_distances = pdist(ncoord)

默认的度量方式是“欧几里得距离”，默认的“p”为2。

在下面的评论中，我错误地提到了pdist的结果是n x n矩阵。 要得到n x n矩阵，您需要执行以下操作：

from scipy.spatial.distance import pdist, squareform

pairwise_distances = squareform(pdist(ncoord))

或者

from scipy.spatial.distance import cdist

pairwise_distances = cdist(ncoord, ncoord)

- lightalchemist

---

我确实做了这件事，但没有放在这里。我的代码的最后一行是：c[j][i]=c[i][j]。 - Fairy

谢谢。但我仍然不理解ncoord[i,j]和ncoord[i][j]之间的区别。 - Fairy

@Fairy，我已经修复了关于对c元素赋值的拼写错误。它们都应该使用相同的格式c[i, j]而不是c[i][j]。 - lightalchemist

c[i][j] 工作正常。它们之间有什么区别？（c[i][j] 和 c[i,j]） - Fairy

@Fairy 当您像这样链接索引时，实际上是首先获取 c 的第 i 行，然后在该行的第一维中索引第 j 个元素。对于 np.array，c[i].shape == (n,)（即它只有一个长度为 n 的维度），因此索引第 j 个元素没有问题。但是，np.matrix 的行为不同，因为 c[i].shape == (1, n)。由于保留了单例行维度，因此当您尝试在 c[i] 的第一维中索引第 j 个元素时，索引会越界。重点是要使用 np.array，而不是 np.matrix！ - ali_m

Answer 3

我猜你想要做的是：你说你想要一个20x20的矩阵...但是你编写的那个是三角形的。

因此，我编写了一个完整的20x20的矩阵。

distances = []
for i in range(len(ncoord)):
    given_i = []
    for j in range(len(ncoord)):
        d_val = math.sqrt((ncoord[i, 0]-ncoord[j,0])**2+(ncoord[i,1]-ncoord[j,1])**2)
        given_i.append(d_val)

    distances.append(given_i)

    # distances[i][j] = distance from i to j

SciPy方式：

from scipy.spatial.distance import cdist
# Isn't scipy nice - can also use pdist... works in the same way but different recall method.
distances = cdist(ncoord, ncoord, 'euclidean')

- Abhinav Ramakrishnan

---

谢谢您的评论。我也会尝试您的方法。 - Fairy

每当你需要在NumPy中通过一个数组进行双重循环时，你就会失去NumPy一开始提供的速度优势。你应该尽可能地使用广播。然而，在某些操作中，我认为包括这个操作，你不能使用广播，因为每一步的值都取决于它们的邻居。在这些情况下，SciPy解决方案通常在C级别上进行了优化（参见cython），因此它们仍然可以更快。我希望cdist函数比双重循环要快得多。 - Adam Hughes

Answer 4

使用自己的自定义平方根和平方和并不总是安全的，它们可能会溢出或下溢。就速度而言，它们是相同的。

np.hypot(
    np.subtract.outer(x, x),
    np.subtract.outer(y, y)
)

下溢

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

参考

- eroot163pi