在python中计算numpy的欧几里得距离答案

【问题标题】：Computing Euclidean distance for numpy in python在python中计算numpy的欧几里得距离
【发布时间】：2015-04-25 13:49:05
【问题描述】：

我是 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)

但是，我收到“输入必须是方阵”的错误 ”。我想知道是否有人知道这里发生了什么。谢谢

【问题讨论】：

请edit您的问题包括ncoord的定义。感谢您提高问题的参考价值并使其更易于回答！
你的 n 是多少？ for j in range(i+1,n-1) 会做j=i+1, i+2, ..., n-2。我猜你希望这两个范围都达到n，而不是n-1。
@MarkG 是的，我有 20 个节点（n=20），我希望两个索引都上升到 n。我尝试了 n 而不是 n-1 但我得到了同样的错误。我可以在 MATLAB 中轻松编写代码，但我必须使用 Python。 Python 中的索引是不同的，所以我可能错了。
那么你的两个 for 循环都应该上升到 n：for i in range(0,n): 和 for j in range(i+1,n): 如果这不是你的错误，那么你需要显示更多代码。
@MarkG 是的，这不是我的错误。我的代码就是我在主要问题中提到的。我什么都没有了

标签： python numpy

【解决方案1】：

为此使用嵌套的for 循环有很多更快的替代方法。我将向您展示两种不同的方法 - 第一种将是一种更通用的方法，将向您介绍广播和矢量化，第二种使用更方便的 scipy 库函数。

1。通用方式，使用广播和矢量化

我建议做的第一件事是改用np.array 而不是np.matrix。 a number of reasons 首选数组，最重要的是因为它们可以具有 >2 维，并且它们使逐元素乘法变得不那么尴尬。

import numpy as np

ncoord = np.array(ncoord)

使用数组，我们可以通过插入新的单例维度和broadcasting 的减法来消除嵌套的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

【讨论】：

这个广播对我来说很神奇。我怎样才能得到一些关于它的直觉？
感谢这个神奇的方法，但它仍然比叉积慢得多，而复杂度看起来是一样的。
在使用方法 1 计算大型矩阵 (1000 * 20000) 时，我也遇到了一些内存问题，而方法 2 (scipy) 则没有。

【解决方案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 array，那么它们将给出相同的结果。

对于一个 Numpy matrix，ncoord[i] 返回ncoord 的第 ith 行，它本身就是一个 Numpy matrix 对象，具有形状在您的情况下为 1 x 2。因此，ncoord[i][j] 实际上意味着：取ncoord 的第 i 行 和取那 1 x 2 的 第 j 行 矩阵。这就是当j > 0 时出现索引问题的地方。

关于分配给c[i][j]“工作”的cmets，它不应该。至少在我构建的 Numpy 1.9.1 中，如果您的索引 i 和 j 迭代到 n，它就不应该工作。

顺便说一句，记得将矩阵c 的转置添加到自身。

建议使用 Numpy 数组而不是矩阵。见this post。

如果您的坐标存储为 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)

【讨论】：

我确实这样做了，但没有放在这里。我的代码的最后一行是：c[j][i]=c[i][j]
谢谢它现在正在工作。然而我现在被误解了。我想当我们想在 Python 中调用矩阵的元素时，我们需要将其称为 a[][]，但您使用的是 a[,]。为什么使用第二种格式从ncoord读取数据，却通过调用c[][]等c元素将距离保存在c矩阵中？
非常感谢您提供完整的信息。我将尝试您提到的另一种方法，看看是否可以获得与结果相同的矩阵大小（我假设 pairwise_distances 是一个 n*n 矩阵）
是的 pairwise_distances 将是一个 n x n 矩阵，如果 n 是您拥有的点数。
谢谢。但是我还是不明白 ncoord[i,j] 和 ncoord[i][j] 的区别

【解决方案3】：

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

因此我编写了一个完整的 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')

【讨论】：

感谢您的评论。我也会试试你的方法。
任何时候你必须在 numpy 中对数组进行双循环，你就会失去 NumPy 首先提供的速度优势。您希望尽可能进行广播。但是，对于某些操作，我认为包括这一操作在内，您不能广播，因为每一步的值都取决于它们的邻居。在这些情况下，SciPy 解决方案通常在 c 级别进行优化（请参阅 cython），因此它们仍然可以更快。我希望 cdist 函数比双循环快得多。

【解决方案4】：

使用您自己的自定义 sqrt sum sqaures 并不总是安全的，它们可能会上溢或下溢。速度方面他们是一样的

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

Refer

【讨论】：

1。通用方式，使用广播和矢量化

2。懒惰的方式，使用pdist

下溢

溢出

没有下溢

无溢出

2。懒惰的方式，使用`pdist`