识别具有最小欧几里得距离的点答案

【问题标题】：Identifying points with the smallest Euclidean distance识别具有最小欧几里得距离的点
【发布时间】：2011-07-04 10:16:32
【问题描述】：

我有一个 n 维点的集合，我想找出最接近的 2 个。我能想到的最好的二维是：

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]

但这对于大型数组来说太慢了。我可以对其应用什么样的优化？

【问题讨论】：

@Ηλίας：你大概有多少分？请注意，可能有超过 2 个点（甚至所有点）的集合具有相同的距离（但不准确的计算可能无法反映这一点，因此最终您需要能够设置阈值 trh，其中距离差低于 trh 是认为相等）。您不想找出离给定点最近的点吗？
@eat 这是我正在构建的层次结构集群，我需要找到两个最近的质心。通常少于一千点，但我需要看看它可以扩展多少。在我的情况下，舍入错误不会那么重要。

标签： python algorithm numpy nearest-neighbor euclidean-distance

【解决方案1】：

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

【讨论】：

从单个整数中获取这些坐标的最简单方法是什么？
@Ηλίας 如果 distances 包含上述 pdist 调用的结果，则可以使用 np.unravel_index(np.argmin(distances), distances.shape)。
使用这种方法在 O(N^2) 时间内找到最接近的配对让我感到很痛苦，因为分而治之 O(N log N) 解决方案实际上是我的第一个算法我在学校的算法课上学过。但这实现起来要容易得多，而且对于足够小的集合也能正常工作。

【解决方案2】：

关于这个问题有一个完整的维基百科页面，请参阅：http://en.wikipedia.org/wiki/Closest_pair_of_points

执行摘要：您可以使用递归分治算法实现 O(n log n)（在上面的 Wiki 页面上进行了概述）。

【讨论】：

整洁！我很高兴我在写之前点击了刷新：“显然复杂度是 O(n^2)”;o)
太棒了。如果要连续添加点，并且要更新最小距离对，那么维护 Delaunay 三角剖分结构是有效的。

【解决方案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)：

【讨论】：

实际上可能在 2D 中工作。你有没有计时？然而，这种方法在更高的昏暗中失败得很惨。谢谢
@eat：你为什么说它“失败得很惨”？ 3D 比 2D 中的相同 N 慢 4-5 倍。但是任何方法（除了天真的蛮力方法）都会看到 D 的放缓。
好吧，尝试在 123D 中进行 Delaunay 三角剖分是没有意义的！所以这不会解决 OP 的问题（除非他的 nD 是 2 或 3）。不要误会，我真的很高兴scipy 能够如此快速地执行 Delaunay 三角剖分。请使用pdist 为 n= 2...123 进行一些计时，你会看到的。谢谢
@eat：我错过了 OP 想要一个通用的 N-D 解决方案这一事实，我的印象是它严格来说是 2D。我有点“桥梁和隧道”，有时认为 3D 不仅是“高维”，而且是最高的！

【解决方案4】：

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

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

输出 N*(N-1)/2 个唯一对（因为 r_ij == r_ji）。然后，您可以搜索最小值并避免代码中的整个循环混乱。

【讨论】：

【解决方案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]))

有了更多的点，您需要能够以某种方式利用聚类的层次结构。

【讨论】：

【解决方案6】：

与仅执行嵌套循环并跟踪最短的对相比，它的速度有多快？我认为创建一个巨大的交叉阵列可能会伤害你。如果你只做二维点，即使 O(n^2) 仍然很快。

【讨论】：

它有帮助，但对于大型矩阵会很快退化

【解决方案7】：

对于小型数据集，可接受的答案是可以的，但它的执行时间缩放为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 获得相同的解决方案，方法是将sklearn 导入替换为以下 Scipy 导入。但请注意，cKDTree 与 BallTree 不同，它不能很好地扩展到高维度

from scipy.spatial import cKDTree as tree

【讨论】：