python - 使用每个矩阵行的可变范围向量化 numpy.argmin 搜索

Question

有没有办法去掉下面代码中的循环并用向量化操作代替？

给定一个数据矩阵，对于每一行，我想找到适合单独数组中定义的范围(每行)的最小值的索引。

这是一个例子:

import numpy as np
np.random.seed(10)

# Values of interest, for this example a random 6 x 100 matrix
data = np.random.random((6,100))

# For each row, define an inclusive min/max range
ranges = np.array([[0.3, 0.4],
                   [0.35, 0.5],
                   [0.45, 0.6],
                   [0.52, 0.65],
                   [0.6,  0.8],
                   [0.75,  0.92]])


# For each row, find the index of the minimum value that fits inside the given range
result = np.zeros(6).astype(np.int)
for i in xrange(6):
    ind = np.where((ranges[i][0] <= data[i]) & (data[i] <= ranges[i][1]))[0]
    result[i] = ind[np.argmin(data[i,ind])]

print result
# Result: [35  8 22  8 34 78]

print data[np.arange(6),result]
# Result: [ 0.30070006  0.35065639  0.45784951  0.52885388  0.61393513  0.75449247]

Answer 1

方法#1:使用 broadcasting 和np.minimum.reduceat -

mask = (ranges[:,None,0] <= data) & (data <= ranges[:,None,1])
r,c = np.nonzero(mask)
cut_idx = np.unique(r, return_index=1)[1]
out = np.minimum.reduceat(data[mask], cut_idx)

避免的改进np.nonzero并计算cut_idx直接来自mask :

cut_idx = np.concatenate(( [0], np.count_nonzero(mask[:-1],1).cumsum() ))

方法#2:使用broadcasting并用 NaNs 填充无效位置然后使用 np.nanargmin -

mask = (ranges[:,None,0] <= data) & (data <= ranges[:,None,1])
result = np.nanargmin(np.where(mask, data, np.nan), axis=1)
out = data[np.arange(6),result]

方法#3:如果您的迭代次数不够(就像示例中具有 6 次迭代的循环一样)，您可能希望坚持使用循环以提高内存效率，但请使用更高效masking用 bool 数组代替 -

out = np.zeros(6)
for i in xrange(6):
    mask_i = (ranges[i,0] <= data[i]) & (data[i] <= ranges[i,1])
    out[i] = np.min(data[i,mask_i])

方法#4:这里还有一种可能的循环解决方案。这个想法是对每行数据进行排序。然后，在np.searchsorted的帮助下，使用每行的两个范围限制来决定开始和停止索引。 。此外，我们将使用这些索引进行切片，然后获取最小值。以这种方式切片的好处是，我们将使用 views因此，无论是在内存还是性能方面都非常高效。

实现看起来像这样 -

out = np.zeros(6)
sdata = np.sort(data, axis=1)
for i in xrange(6):
    start = np.searchsorted(sdata[i], ranges[i,0])
    stop = np.searchsorted(sdata[i], ranges[i,1], 'right')    
    out[i] = np.min(sdata[i,start:stop])

此外，我们可以在执行 vectorized searchsorted 后以矢量化方式获取这些开始、停止索引。 .

基于 @Daniel F 的建议对于我们处理给定data范围内的范围的情况，我们可以简单地使用起始索引 -

out[i] = sdata[i, start]