两个不同Numpy数组中点的最小欧氏距离，不包括内部

Question

我有两个数组，里面存的是 x-y 坐标，我想找出一个数组中的每个点和另一个数组中所有点之间的最小欧几里得距离。注意，这两个数组的大小不一定相同。例如：

xy1=numpy.array(
[[  243,  3173],
[  525,  2997]])

xy2=numpy.array(
[[ 682, 2644],
[ 277, 2651],
[ 396, 2640]])

我现在的方法是遍历第一个数组 xy1 中的每个坐标，然后计算这个坐标和第二个数组中其他坐标之间的距离。

mindist=numpy.zeros(len(xy1))
minid=numpy.zeros(len(xy1))

for i,xy in enumerate(xy1):
    dists=numpy.sqrt(numpy.sum((xy-xy2)**2,axis=1))
    mindist[i],minid[i]=dists.min(),dists.argmin()

有没有办法可以去掉这个循环，直接在两个数组之间进行逐个元素的计算呢？我想象中可以生成一个距离矩阵，然后从中找到每一行或每一列的最小值。

还有一种看待这个问题的方法。假设我把 xy1（长度为 m）和 xy2（长度为 p）合并成一个新的数组 xy（长度为 n），并且我记录下原始数组的长度。理论上，我应该能够从这些坐标生成一个 n x n 的距离矩阵，然后从中提取出一个 m x p 的子矩阵。有没有什么高效的方法来生成这个子矩阵呢？

Answer 1

这个被接受的答案并没有完全解决问题，因为问题是要找出两组点之间的最小距离，而不是两组中每一个点之间的距离。

虽然直接解决原问题的方法确实是计算每一对点之间的距离，然后找出最小的那个，但如果我们只关心最小距离，其实不需要这样做。对于后者的问题，有一种更快的解决方案。

所有提出的解决方案的运行时间都是按 m*p = len(xy1)*len(xy2) 来计算的。对于小数据集来说，这样是可以的，但我们可以写出一个更优的解决方案，它的运行时间是 m*log(p)，这样在处理大的 xy2 数据集时能节省很多时间。

要实现这个最佳的执行时间，可以使用 scipy.spatial.KDTree，具体方法如下：

import numpy as np
from scipy import spatial

xy1 = np.array(
    [[243,  3173],
     [525,  2997]])

xy2 = np.array(
    [[682, 2644],
     [277, 2651],
     [396, 2640]])

# This solution is optimal when xy2 is very large
tree = spatial.KDTree(xy2)
mindist, minid = tree.query(xy1)
print(mindist)

# This solution by @denis is OK for small xy2
mindist = np.min(spatial.distance.cdist(xy1, xy2), axis=1)
print(mindist)

这里的 mindist 是 xy1 中每个点与 xy2 中点的最小距离。

Answer 2

要计算一个 m 行 p 列的距离矩阵，可以这样做：

>>> def distances(xy1, xy2):
...   d0 = numpy.subtract.outer(xy1[:,0], xy2[:,0])
...   d1 = numpy.subtract.outer(xy1[:,1], xy2[:,1])
...   return numpy.hypot(d0, d1)

这里的 .outer 函数会生成两个矩阵，这两个矩阵是沿着两个方向的数值差异。然后，.hypot 函数会把这些差异转换成一个形状相同的矩阵，这个矩阵里面的值就是欧几里得距离，也就是我们常说的直线距离。

Answer 3

（几个月后）
scipy.spatial.distance.cdist( X, Y ) 这个函数可以计算所有点对之间的距离，
无论是二维、三维还是更多维度的点都可以使用。
它还支持22种不同的距离计算方式，详细信息可以在这里找到。

# cdist example: (nx,dim) (ny,dim) -> (nx,ny)

from __future__ import division
import sys
import numpy as np
from scipy.spatial.distance import cdist

#...............................................................................
dim = 10
nx = 1000
ny = 100
metric = "euclidean"
seed = 1

    # change these params in sh or ipython: run this.py dim=3 ...
for arg in sys.argv[1:]:
    exec( arg )
np.random.seed(seed)
np.set_printoptions( 2, threshold=100, edgeitems=10, suppress=True )

title = "%s  dim %d  nx %d  ny %d  metric %s" % (
        __file__, dim, nx, ny, metric )
print "\n", title

#...............................................................................
X = np.random.uniform( 0, 1, size=(nx,dim) )
Y = np.random.uniform( 0, 1, size=(ny,dim) )
dist = cdist( X, Y, metric=metric )  # -> (nx, ny) distances
#...............................................................................

print "scipy.spatial.distance.cdist: X %s Y %s -> %s" % (
        X.shape, Y.shape, dist.shape )
print "dist average %.3g +- %.2g" % (dist.mean(), dist.std())
print "check: dist[0,3] %.3g == cdist( [X[0]], [Y[3]] ) %.3g" % (
        dist[0,3], cdist( [X[0]], [Y[3]] ))


# (trivia: how do pairwise distances between uniform-random points in the unit cube
# depend on the metric ? With the right scaling, not much at all:
# L1 / dim      ~ .33 +- .2/sqrt dim
# L2 / sqrt dim ~ .4 +- .2/sqrt dim
# Lmax / 2      ~ .4 +- .2/sqrt dim