向量化为另一个数组中的每个元素查找数组中最接近的值答案

【问题标题】：Vectorize finding closest value in an array for each element in another array向量化为另一个数组中的每个元素查找数组中最接近的值
【发布时间】：2014-01-13 19:48:15
【问题描述】：

输入

known_array : numpy 数组；仅由标量值组成； shape: (m, 1)

test_array : numpy 数组；仅由标量值组成； shape: (n, 1)

输出

indices : numpy 数组； shape: (n, 1);对于test_array 中的每个值，查找known_array 中最接近值的索引

residual : numpy 数组； shape: (n, 1);对于test_array 中的每个值，找出与known_array 中最接近的值的差

示例

In [17]: known_array = np.array([random.randint(-30,30) for i in range(5)])

In [18]: known_array
Out[18]: array([-24, -18, -13, -30,  29])

In [19]: test_array = np.array([random.randint(-10,10) for i in range(10)])

In [20]: test_array
Out[20]: array([-6,  4, -6,  4,  8, -4,  8, -6,  2,  8])

示例实现（未完全矢量化）

def find_nearest(known_array, value):
    idx = (np.abs(known_array - value)).argmin()
    diff = known_array[idx] - value
    return [idx, -diff]

In [22]: indices = np.zeros(len(test_array))

In [23]: residual = np.zeros(len(test_array))

In [24]: for i in range(len(test_array)):
   ....:     [indices[i], residual[i]] =  find_nearest(known_array, test_array[i])
   ....:     

In [25]: indices
Out[25]: array([ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.])

In [26]: residual
Out[26]: array([  7.,  17.,   7.,  17.,  21.,   9.,  21.,   7.,  15.,  21.])

加快这项任务的最佳方法是什么？ Cython 是一种选择，但是，我总是希望能够删除 for 循环并让代码保持纯 NumPy。

注意：咨询了以下 Stack Overflow 问题

更新

我做了一些小的基准来比较非矢量化和矢量化解决方案（接受的答案）。

In [48]: [indices1, residual1] = find_nearest_vectorized(known_array, test_array)

In [53]: [indices2, residual2] = find_nearest_non_vectorized(known_array, test_array)

In [54]: indices1==indices2
Out[54]: array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True],   dtype=bool)

In [55]: residual1==residual2
Out[55]: array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)

In [56]: %timeit [indices2, residual2] = find_nearest_non_vectorized(known_array, test_array)
10000 loops, best of 3: 173 µs per loop

In [57]: %timeit [indices1, residual1] = find_nearest_vectorized(known_array, test_array)
100000 loops, best of 3: 16.8 µs per loop

约 10 倍 加速！

澄清

known_array 未排序。

我运行了下面@cyborg 回答中给出的基准。

案例 1：如果 known_array 已排序

known_array = np.arange(0,1000)
test_array = np.random.randint(0, 100, 10000)
print('Speedups:')
base_time = time_f('base')
for func_name in ['diffs', 'searchsorted1', 'searchsorted2']:
    print func_name + ' is x%.1f faster than base.' % (base_time / time_f(func_name))
    assert np.allclose(base(known_array, test_array), eval(func_name+'(known_array, test_array)'))

Speedups:
diffs is x0.4 faster than base.
searchsorted1 is x81.3 faster than base.
searchsorted2 is x107.6 faster than base.

首先，对于大型数组，diffs 方法实际上速度较慢，它还占用了大量 RAM，当我在实际数据上运行时我的系统挂起。

情况2：当known_array未排序时；代表实际场景

known_array = np.random.randint(0,100,100)
test_array = np.random.randint(0, 100, 100)

Speedups:
diffs is x8.9 faster than base.
AssertionError                            Traceback (most recent call last)
<ipython-input-26-3170078c217a> in <module>()
      5 for func_name in ['diffs', 'searchsorted1', 'searchsorted2']:
      6     print func_name + ' is x%.1f faster than base.' % (base_time /  time_f(func_name))
----> 7     assert np.allclose(base(known_array, test_array),  eval(func_name+'(known_array, test_array)'))

AssertionError: 


searchsorted1 is x14.8 faster than base.

我还必须评论说，这种方法也应该是内存效率的。否则我的 8 GB RAM 是不够的。在基本情况下，这很容易就足够了。

【问题讨论】：

你的数据是否排序没关系； HYRY 发布的方法处理了这种情况，并且具有线性而不是 diff 方法的二次内存性能；他的答案应该被标记为正确的

标签： python algorithm numpy vectorization cython

【解决方案1】：

如果数组很大，应该使用searchsorted:

import numpy as np
np.random.seed(0)
known_array = np.random.rand(1000)
test_array = np.random.rand(400)

%%time
differences = (test_array.reshape(1,-1) - known_array.reshape(-1,1))
indices = np.abs(differences).argmin(axis=0)
residual = np.diagonal(differences[indices,])

输出：

CPU times: user 11 ms, sys: 15 ms, total: 26 ms
Wall time: 26.4 ms

searchsorted版本：

%%time

index_sorted = np.argsort(known_array)
known_array_sorted = known_array[index_sorted]

idx1 = np.searchsorted(known_array_sorted, test_array)
idx2 = np.clip(idx1 - 1, 0, len(known_array_sorted)-1)

diff1 = known_array_sorted[idx1] - test_array
diff2 = test_array - known_array_sorted[idx2]

indices2 = index_sorted[np.where(diff1 <= diff2, idx1, idx2)]
residual2 = test_array - known_array[indices]

输出：

CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 311 µs

我们可以检查结果是否相同：

assert np.all(residual == residual2)
assert np.all(indices == indices2)

【讨论】：

我不会安静地从这两种方法中得到相同的结果
我在问题末尾添加了一个澄清部分来解释这一点
搜索排序算法对我来说效果很好，但如果test_array 的任何值大于known_array 的最大值，它就会失败。在这种情况下，np.searchsorted 将返回一个 1 太大而无法用作known_array_sorted 的索引的索引。我已经编辑了上面的答案来解决这个问题。请检查它是否适合您！

【解决方案2】：

TL;DR：使用numpy.searchsorted()。

import inspect
from timeit import timeit
import numpy as np

known_array = np.arange(-10, 10)
test_array = np.random.randint(-10, 10, 1000)
number = 1000

def base(known_array, test_array):
    def find_nearest(known_array, value):
        idx = (np.abs(known_array - value)).argmin()
        return idx
    indices = np.zeros_like(test_array, dtype=known_array.dtype)
    for i in range(len(test_array)):
        indices[i] =  find_nearest(known_array, test_array[i])
    return indices

def diffs(known_array, test_array):
    differences = (test_array.reshape(1,-1) - known_array.reshape(-1,1))
    indices = np.abs(differences).argmin(axis=0)
    return indices

def searchsorted1(known_array, test_array):
    index_sorted = np.argsort(known_array)
    known_array_sorted = known_array[index_sorted]
    idx1 = np.searchsorted(known_array_sorted, test_array)
    idx1[idx1 == len(known_array)] = len(known_array)-1
    idx2 = np.clip(idx1 - 1, 0, len(known_array_sorted)-1)
    diff1 = known_array_sorted[idx1] - test_array
    diff2 = test_array - known_array_sorted[idx2]
    indices2 = index_sorted[np.where(diff1 <= diff2, idx1, idx2)]
    return indices2

def searchsorted2(known_array, test_array):
    index_sorted = np.argsort(known_array)
    known_array_sorted = known_array[index_sorted]
    known_array_middles = known_array_sorted[1:] - np.diff(known_array_sorted.astype('f'))/2
    idx1 = np.searchsorted(known_array_middles, test_array)
    indices = index_sorted[idx1]
    return indices

def time_f(func_name):
    return timeit(func_name+"(known_array, test_array)",
        'from __main__ import known_array, test_array, ' + func_name, number=number)

print('Speedups:')
base_time = time_f('base')
for func_name in ['diffs', 'searchsorted1', 'searchsorted2']:
    print func_name + ' is x%.1f faster than base.' % (base_time / time_f(func_name))

输出：

Speedups:
diffs is x29.9 faster than base.
searchsorted1 is x37.4 faster than base.
searchsorted2 is x64.3 faster than base.

【讨论】：

能否提供一个符合我提供的测试数据的示例代码？
known_array未排序时断言失败；例如。 known_array = np.random.randint(0, 100, 100) test_array = np.random.randint(0, 100, 1000)
我在问题的末尾添加了一个澄清部分来说明这一点
searchsorted1 和 searchsorted2 均失败
我删除了断言，因为它没有定义如果有多个最接近的值应该选择 known_array 的哪个值。

【解决方案3】：

例如，您可以使用以下方法计算所有差异：

differences = (test_array.reshape(1,-1) - known_array.reshape(-1,1))

并使用argmin 和花哨的索引以及np.diagonal 来获得所需的索引和差异：

indices = np.abs(differences).argmin(axis=0)
residual = np.diagonal(differences[indices,])

所以对于

>>> known_array = np.array([-24, -18, -13, -30,  29])
>>> test_array = np.array([-6,  4, -6,  4,  8, -4,  8, -6,  2,  8])

一个得到

>>> indices
array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2])
>>> residual
array([ 7, 17,  7, 17, 21,  9, 21,  7, 15, 21])

【讨论】：

解决方案有效！我只是在做一些小的基准测试来比较性能。如果您能突出显示您对矢量化此任务的看法，那就太好了。那会更有帮助！
在问题的更新部分添加了基准。讨论两个版本的内存要求也可能会有所帮助。

【解决方案4】：

我能想到的最快的。这需要对 array 进行排序。值可以是标量或列表/数组。

def find_nearest(value, array):
    idx = np.searchsorted(array, value, side="left")
    cond = np.logical_and(idx>0, np.logical_or(idx == len(array), np.fabs(value - array[idx-1]) < np.fabs(value - array[idx])))
    return np.where(cond, array[idx-1], array[idx])

【讨论】：