加快阵列中所有可能对之间的距离

Question

我有几个 (~10^10) 点的 x、y、z 坐标数组（此处仅显示 5 个）

a= [[ 34.45  14.13   2.17]
    [ 32.38  24.43  23.12]
    [ 33.19   3.28  39.02]
    [ 36.34  27.17  31.61]
    [ 37.81  29.17  29.94]]

我想创建一个新数组，其中仅包含与列表中所有其他点至少有一定距离d 的那些点。我用while循环写了一段代码，

 import numpy as np
 from scipy.spatial import distance 

 d=0.1 #or some distance 
 i=0
 selected_points=[]
 while i < len(a):
          interdist=[]  
          j=i+1
          while j<len(a):
              interdist.append(distance.euclidean(a[i],a[j]))
              j+=1

          if all(dis >= d for dis in interdist):
              np.array(selected_points.append(a[i]))
          i+=1

这可行，但执行此计算需要很长时间。我在某处读到while 循环非常慢。

我想知道是否有人对如何加快计算速度有任何建议。

编辑：虽然我找到与所有其他粒子至少有一定距离的粒子的目标保持不变，但我刚刚意识到我的代码存在严重缺陷，假设我有3个粒子，我的代码执行以下操作，对于i的第一次迭代，它计算距离1->2，1->3，假设1->2小于阈值距离d，所以代码抛出远离粒子1。对于i的下一次迭代，它只做2->3，假设它发现它大于d，所以它保持粒子2，但这是错误的！因为 2 也应该与粒子 1 一起丢弃。 @svohara 的解决方案是正确的！

Answer 1

对于大数据集和低维点（例如您的 3 维数据），有时使用空间索引方法有很大的好处。低维数据的一种流行选择是 k-d 树。

策略是对数据集进行索引。然后使用相同的数据集查询索引，以返回每个点的 2 个最近邻。第一个最近的邻居总是点本身（dist=0），所以我们真的想知道下一个最近的点有多远（第二个最近的邻居）。对于那些 2-NN > 阈值的点，您就有了结果。

from scipy.spatial import cKDTree as KDTree
import numpy as np

#a is the big data as numpy array N rows by 3 cols
a = np.random.randn(10**8, 3).astype('float32')

# This will create the index, prepare to wait...
# NOTE: took 7 minutes on my mac laptop with 10^8 rand 3-d numbers
#  there are some parameters that could be tweaked for faster indexing,
#  and there are implementations (not in scipy) that can construct
#  the kd-tree using parallel computing strategies (GPUs, e.g.)
k = KDTree(a)

#ask for the 2-nearest neighbors by querying the index with the
# same points
(dists, idxs) = k.query(a, 2)
# (dists, idxs) = k.query(a, 2, n_jobs=4)  # to use more CPUs on query...

#Note: 9 minutes for query on my laptop, 2 minutes with n_jobs=6
# So less than 10 minutes total for 10^8 points.

# If the second NN is > thresh distance, then there is no other point
# in the data set closer.
thresh_d = 0.1   #some threshold, equiv to 'd' in O.P.'s code
d_slice = dists[:, 1]  #distances to second NN for each point
res = np.flatnonzero( d_slice >= thresh_d )

Answer 2

这是使用distance.pdist 的矢量化方法-

# Store number of pts (number of rows in a)
m = a.shape[0]

# Get the first of pairwise indices formed with the pairs of rows from a
# Simpler version, but a bit slow : idx1,_ = np.triu_indices(m,1)
shifts_arr = np.zeros(m*(m-1)/2,dtype=int)
shifts_arr[np.arange(m-1,1,-1).cumsum()] = 1
idx1 = shifts_arr.cumsum()

# Get the IDs of pairs of rows that are more than "d" apart and thus select 
# the rest of the rows using a boolean mask created with np.in1d for the 
# entire range of number of rows in a. Index into a to get the selected points.
selected_pts = a[~np.in1d(np.arange(m),idx1[distance.pdist(a) < d])] 

对于像10e10 这样的庞大数据集，我们可能必须根据可用的系统内存分块执行操作。

Answer 3

您的算法是二次的（10^20 次操作），如果分布几乎是随机的，这是一种线性方法。 将您的空间分割成d/sqrt(3)^3 大小的盒子。把每个点放在它的盒子里。

那么对于每个盒子，

• 如果只有一个点，您只需计算与小邻域中的点的距离。

• 否则无事可做。

Answer 4

1. 放弃追加，它一定很慢。您可以拥有一个静态的距离向量，并使用 [] 将数字放在正确的位置。

2. 使用 min 而不是 all。您只需要检查最小距离是否大于 x。

3. 实际上，您可以在发现距离小于限制的那一刻中断追加，然后您可以丢弃两个点。这样，您甚至不必节省任何距离（除非您以后需要它们）。

   4. 由于 d(a,b)=d(b,a) 您只能对以下点进行内部循环，忘记您已经计算的距离。如果您需要它们，您可以从数组中选择更快的。

从您的评论来看，如果您没有重复的观点，我相信这样做可以。

selected_points = []
for p1 in a:
    save_point = True
    for p2 in a:
        if p1!=p2 and distance.euclidean(p1,p2)<d:
            save_point = False
            break
    if save_point:
        selected_points.append(p1)

return selected_points

最后我检查了 a,b 和 b,a，因为您不应该在处理列表时修改它，但是使用一些附加变量可以更聪明。