在Python中查找两个列表/数组中最近的项

Question

我有两个numpy数组，x和y，里面存的是浮点数。对于x中的每一个值，我想找到y中最接近的元素，而且不能重复使用y中的元素。最终的结果应该是x中元素的索引和y中元素的索引之间的一一对应关系。这里有一种不好的方法，它依赖于排序。这个方法会把每个配对的元素从列表中移除。如果不排序，这种方法就不好用了，因为配对的结果会受到原始输入数组顺序的影响。

def min_i(values):
    min_index, min_value = min(enumerate(values),
                               key=operator.itemgetter(1))
    return min_index, min_value

# unsorted elements
unsorted_x = randn(10)*10
unsorted_y = randn(10)*10

# sort lists
x = sort(unsorted_x)
y = sort(unsorted_y)

pairs = []
indx_to_search = range(len(y))

for x_indx, x_item in enumerate(x):
    if len(indx_to_search) == 0:
        print "ran out of items to match..."
        break
    # until match is found look for closest item
    possible_values = y[indx_to_search]
    nearest_indx, nearest_item = min_i(possible_values)
    orig_indx = indx_to_search[nearest_indx]
    # remove it
    indx_to_search.remove(orig_indx)
    pairs.append((x_indx, orig_indx))
print "paired items: "
for k,v in pairs:
    print x[k], " paired with ", y[v]

我更希望在不先排序的情况下完成这个任务，但如果它们已经排序了，我想得到原始未排序列表unsorted_x和unsorted_y中的索引。用numpy、scipy、Python或者pandas，最好的方法是什么呢？谢谢。

补充说明：为了澄清，我并不是想在所有元素中找到最佳匹配（比如说最小化距离的总和），而是想为每个元素找到最佳匹配，偶尔牺牲其他元素也是可以的。我假设y通常比x大很多，和上面的例子相反，因此对于x中的每个值，在y中通常会有很多非常好的匹配，我只想高效地找到其中一个。

有人能给我展示一下如何用scipy的kd树来实现这个吗？文档上相关的信息很少。

kdtree = scipy.spatial.cKDTree([x,y])
kdtree.query([-3]*10) # ?? unsure about what query takes as arg

Answer 1

编辑 2 使用 KDTree 的方法如果能选择一个合适的邻居数量，可以让每个数组里的项目都有唯一的邻居，这样效果会很好。用下面的代码：

def nearest_neighbors_kd_tree(x, y, k) :
    x, y = map(np.asarray, (x, y))
    tree =scipy.spatial.cKDTree(y[:, None])    
    ordered_neighbors = tree.query(x[:, None], k)[1]
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    nearest_neighbor.fill(-1)
    used_y = set()
    for j, neigh_j in enumerate(ordered_neighbors) :
        for k in neigh_j :
            if k not in used_y :
                nearest_neighbor[j] = k
                used_y.add(k)
                break
    return nearest_neighbor

再加上一个包含 n=1000 个点的样本，我得到：

In [9]: np.any(nearest_neighbors_kd_tree(x, y, 12) == -1)
Out[9]: True

In [10]: np.any(nearest_neighbors_kd_tree(x, y, 13) == -1)
Out[10]: False

所以最佳的邻居数量是 k=13，然后时间是：

In [11]: %timeit nearest_neighbors_kd_tree(x, y, 13)
100 loops, best of 3: 9.26 ms per loop

但在最糟糕的情况下，你可能需要 k=1000，那样的话：

In [12]: %timeit nearest_neighbors_kd_tree(x, y, 1000)
1 loops, best of 3: 424 ms per loop

这比其他选项要慢：

In [13]: %timeit nearest_neighbors(x, y)
10 loops, best of 3: 60 ms per loop

In [14]: %timeit nearest_neighbors_sorted(x, y)
10 loops, best of 3: 47.4 ms per loop

---

编辑 在搜索之前对数组进行排序，对于超过1000个项目的数组来说是很有帮助的：

def nearest_neighbors_sorted(x, y) :
    x, y = map(np.asarray, (x, y))
    y_idx = np.argsort(y)
    y = y[y_idx]
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    for j, xj in enumerate(x) :
        idx = np.searchsorted(y, xj)
        if idx == len(y) or idx != 0 and y[idx] - xj > xj - y[idx-1] :
            idx -= 1
        nearest_neighbor[j] = y_idx[idx]
        y = np.delete(y, idx)
        y_idx = np.delete(y_idx, idx)
    return nearest_neighbor

对于一个有10000个元素的数组：

In [2]: %timeit nearest_neighbors_sorted(x, y)
1 loops, best of 3: 557 ms per loop

In [3]: %timeit nearest_neighbors(x, y)
1 loops, best of 3: 1.53 s per loop

对于较小的数组，效果稍微差一些。

---

你需要遍历所有项目来实现你的 贪心 最近邻算法，至少要去掉重复项。考虑到这一点，这是我能想到的最快的方法：

def nearest_neighbors(x, y) :
    x, y = map(np.asarray, (x, y))
    y = y.copy()
    y_idx = np.arange(len(y))
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    for j, xj in enumerate(x) :
        idx = np.argmin(np.abs(y - xj))
        nearest_neighbor[j] = y_idx[idx]
        y = np.delete(y, idx)
        y_idx = np.delete(y_idx, idx)

    return nearest_neighbor

现在用：

n = 1000
x = np.random.rand(n)
y = np.random.rand(2*n)

我得到：

In [11]: %timeit nearest_neighbors(x, y)
10 loops, best of 3: 52.4 ms per loop