通过识别欧几里得距离最小的点来解决问题

Question

我有一组n维点，想找出其中距离最近的两个点。对于二维情况，我能想到的最好方法是：

from numpy import *
myArr = array( [[1, 2],
                [3, 4],
                [5, 6],
                [7, 8]] )

n = myArr.shape[0]
cross = [[sum( ( myArr[i] - myArr[j] ) ** 2 ), i, j]
         for i in xrange( n )
         for j in xrange( n )
         if i != j
         ]

print min( cross )

这提供了

[8, 0, 1]

但是对于大型数组来说，这种方法速度太慢了。我可以应用什么样的优化方案呢？

相关问题：

---

计算两个不同Numpy数组中点之间的欧几里得距离，而不是同一数组内的点

- Ηλίας

---

@Ηλίας：你大概有多少个点？请注意，可能存在一组距离相同的点（甚至所有点），但是不准确的计算可能无法反映这一点，因此最终需要能够设置一个阈值trh，其中距离差小于trh被视为相等。你不想找到给定点的最近点吗？ - eat

@eat 我正在构建一个层次聚类，需要找到两个最接近的质心。通常少于一千个数据点，但我需要看看它能扩展多少。在我的情况下，舍入误差不是很重要。 - Ηλίας

Answer 1

尝试使用scipy.spatial.distance.pdist(myArr)。这将为您提供压缩的距离矩阵。您可以在其上使用argmin并找到最小值的索引。这可以转换为成对信息。

- tkerwin

---

从那个单一的整数中获取这些坐标的最简单方法是什么？ - Ηλίας

@Ηλίας 如果 distances 包含上面 pdist 调用的结果，您可以使用 np.unravel_index(np.argmin(distances), distances.shape)。 - sffc

使用这种O(N^2)时间复杂度的方法来寻找最近对让我感到不适，因为在我的算法课上，分治的O(N log N)解法实际上是我学习的第一个算法。但这种方法实现起来就容易多了，而且对于足够小的数据集来说也运行得很好。 - sffc

Answer 2

这个问题在维基百科上有专门的页面，链接如下：http://en.wikipedia.org/wiki/Closest_pair_of_points。简单总结一下：你可以使用递归分治算法（在上述维基页面中有概述）来实现O(n log n)的时间复杂度。

- payne

---

太好了！我很高兴在写“显然复杂度是O(n^2)”之前刷新了一下；o) - das_weezul

很好。如果要逐个添加点，并更新最小距离对，则维护Delaunay三角剖分结构是高效的。 - Alexandre C.

Answer 3

你可以利用最新版本的SciPy(v0.9)的Delaunay三角剖分工具。你可以确保最接近的两个点将成为三角剖分中简单形式的边缘，这比对每个组合进行操作的子集要小得多。

下面是更新后的通用N-D代码：

import numpy
from scipy import spatial

def closest_pts(pts):
    # set up the triangluataion
    # let Delaunay do the heavy lifting
    mesh = spatial.Delaunay(pts)

    # TODO: eliminate reduncant edges (numpy.unique?)
    edges = numpy.vstack((mesh.vertices[:,:dim], mesh.vertices[:,-dim:]))

    # the rest is easy
    x = mesh.points[edges[:,0]]
    y = mesh.points[edges[:,1]]

    dists = numpy.sum((x-y)**2, 1)
    idx = numpy.argmin(dists)

    return edges[idx]
    #print 'distance: ', dists[idx]
    #print 'coords:\n', pts[closest_verts]

dim = 3
N = 1000*dim
pts = numpy.random.random(N).reshape(N/dim, dim)

似乎是O(n)的复杂度：

- Paul

---

可能在二维中确实有效。你做过任何计时吗？然而，这种方法在更高的维度中失败得很惨。谢谢。 - eat

@eat：你为什么说它“失败得很惨”? 3D 比同样的 N 在 2D 中慢 4-5 倍。但是任何方法（除了天真的暴力方法）在 D 中都会看到减速。 - Paul

嗯，在123D中尝试进行Delaunay三角剖分有点毫无意义！因此，这不会解决OP的问题（除非他的nD是2或3）。别误会，我实际上非常高兴scipy能够如此快速地执行Delaunay三角剖分。请使用pdist进行n = 2…123的一些计时，你会看到的。谢谢。 - eat

@eat：我错过了OP想要一个通用的N-D解决方案这一事实，我以为它严格是2D。有时候我有点“桥隧”思维，认为3D不仅是“高维”的，而且是最高的！ - Paul

Answer 4

有一个Scipy函数pdist可以以相当高效的方式获取数组中点之间的成对距离：

http://docs.scipy.org/doc/scipy/reference/spatial.distance.html

该函数输出N*(N-1)/2个唯一对（因为r_ij == r_ji）。您可以搜索最小值并避免在代码中编写整个循环。

- JoshAdel

Answer 5

也许你可以沿着这些线路继续进行：

In []: from scipy.spatial.distance import pdist as pd, squareform as sf
In []: m= 1234
In []: n= 123
In []: p= randn(m, n)
In []: d= sf(pd(p))
In []: a= arange(m)
In []: d[a, a]= d.max()
In []: where(d< d.min()+ 1e-9)
Out[]: (array([701, 730]), array([730, 701]))

如果你有更多的数据点，你需要能够以某种方式利用聚类的层次结构。

- eat

Answer 6

对于小数据集，接受的答案是可以的，但其执行时间随着 n**2 的增长而增加。然而，正如 @payne 指出的那样，最优解可以实现 n*log(n) 的计算时间缩放。

可以使用 sklearn.neighbors.BallTree 获得此最佳解决方案，具体操作如下。

import matplotlib.pyplot as plt
import numpy as np
from sklearn.neighbors import BallTree as tree

n = 10
dim = 2
xy = np.random.uniform(size=[n, dim])

# This solution is optimal when xy is very large
res = tree(xy)
dist, ids = res.query(xy, 2)
mindist = dist[:, 1]  # second nearest neighbour
minid = np.argmin(mindist)

plt.plot(*xy.T, 'o')
plt.plot(*xy[ids[minid]].T, '-o')

这个过程适用于非常大的xy值集合，甚至适用于大维度dim（尽管示例说明了dim=2的情况）。生成的输出如下所示

使用scipy.spatial.cKDTree也可以得到相同的解决方案，只需用以下Scipy导入替换sklearn导入即可。但需要注意的是，cKDTree与BallTree不同，在高维情况下不具有良好的可扩展性。

from scipy.spatial import cKDTree as tree

- divenex

Answer 7

相比于只做嵌套循环并跟踪最短对，它有多快？我认为创建一个巨大的交叉数组可能是在伤害你。即使你只对二维点进行操作，O(n^2)仍然相当迅速。

- Anthony -GISCOE-

---

它有所帮助，但对于大矩阵很快就会退化。 - Ηλίας