【问题标题】:How to speed up the performance of array masking from the results of numpy.searchsorted in python?python - 如何从python中numpy.searchsorted的结果加快数组屏蔽的性能?
【发布时间】:2020-10-23 09:04:29
【问题描述】:

我想从 numpy.searchsorted() 的结果中生成一个掩码:

import numpy as np

# generate test examples
x = np.random.rand(1000000)
y = np.random.rand(200)

# sort x
idx = np.argsort(x)
sorted_x = np.take_along_axis(x, idx, axis=-1)

# searchsort y in x
pt = np.searchsorted(sorted_x, y)

pt 是一个数组。然后我想创建一个大小为(200, 1000000) 的布尔掩码,当它的索引为idx[0:pt[i]] 时,它的值为真,我想出了一个这样的for循环:

mask = np.zeros((200, 1000000), dtype='bool')
for i in range(200):
     mask[i, idx[0:pt[i]]] = True

有人有加快 for 循环的想法吗?

【问题讨论】:

    标签: python performance numpy


    【解决方案1】:

    方法#1

    根据从OP's comments 获取的新信息,指出只有y 是实时变化的,我们可以对x 周围的大量内容进行预处理,从而做得更好。我们将创建一个散列数组来存储阶梯掩码。对于涉及y 的部分,我们将简单地使用从searchsorted 获得的索引来索引散列数组,这将近似于最终的掩码数组。鉴于 numba 的参差不齐的性质,分配剩余布尔值的最后一步可以卸载。如果我们决定扩大y 的长度,这也应该是有益的。

    让我们看看实现。

    使用x 进行预处理:

    sidx = x.argsort()
    ssidx = x.argsort().argsort()
    
    # Choose a scale factor. 
    # 1. A small one would store more mapping info, hence faster but occupy more mem
    # 2. A big one would store less mapping info, hence slower, but memory efficient.
    scale_factor = 100
    mapar = np.arange(0,len(x),scale_factor)[:,None] > ssidx
    

    y 的剩余步骤:

    import numba as nb
    
    @nb.njit(parallel=True,fastmath=True)
    def array_masking3(out, starts, idx, sidx):
        N = len(out)
        for i in nb.prange(N):
            for j in nb.prange(starts[i], idx[i]):
                out[i,sidx[j]] = True
        return out
    
    idx = np.searchsorted(x,y,sorter=sidx)
    s0 = idx//scale_factor
    starts = s0*scale_factor
    out = mapar[s0]
    out = array_masking3(out, starts, idx, sidx)
    

    基准测试

    In [2]: x = np.random.rand(1000000)
       ...: y = np.random.rand(200)
    
    In [3]: ## Pre-processing step with "x"
       ...: sidx = x.argsort()
       ...: ssidx = x.argsort().argsort()
       ...: scale_factor = 100
       ...: mapar = np.arange(0,len(x),scale_factor)[:,None] > ssidx
    
    In [4]: %%timeit
       ...: idx = np.searchsorted(x,y,sorter=sidx)
       ...: s0 = idx//scale_factor
       ...: starts = s0*scale_factor
       ...: out = mapar[s0]
       ...: out = array_masking3(out, starts, idx, sidx)
    41 ms ± 141 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    
     # A 1/10th smaller hashing array has similar timings
    In [7]: scale_factor = 1000
       ...: mapar = np.arange(0,len(x),scale_factor)[:,None] > ssidx
    
    In [8]: %%timeit
       ...: idx = np.searchsorted(x,y,sorter=sidx)
       ...: s0 = idx//scale_factor
       ...: starts = s0*scale_factor
       ...: out = mapar[s0]
       ...: out = array_masking3(out, starts, idx, sidx)
    40.6 ms ± 196 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    
    # @silgon's soln    
    In [5]: %timeit x[np.newaxis,:] < y[:,np.newaxis]
    138 ms ± 896 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
    

    方法#2

    这很好地借鉴了OP's solution

    import numba as nb
    
    @nb.njit(parallel=True)
    def array_masking2(mask1D, mask_out, idx, pt):
        n = len(idx)
        for j in nb.prange(len(pt)):
            if mask1D[j]:
                for i in nb.prange(pt[j],n):
                    mask_out[j, idx[i]] = False
            else:
                for i in nb.prange(pt[j]):
                    mask_out[j, idx[i]] = True
        return mask_out
    
    def app2(idx, pt):
        m,n = len(pt), len(idx)      
        mask1 = pt>len(x)//2
        mask2 = np.broadcast_to(mask1[:,None], (m,n)).copy()
        return array_masking2(mask1, mask2, idx, pt)
    

    因此,想法是,我们有超过一半的索引要设置True,我们切换到设置False,而不是在将这些行预先分配为所有True。这会减少内存访问,从而显着提升性能。

    基准测试

    OP的解决方案:

    @nb.njit(parallel=True,fastmath=True)
    def array_masking(mask, idx, pt):
        for j in nb.prange(pt.shape[0]):
            for i in nb.prange(pt[j]):
                mask[j, idx[i]] = True
        return mask
    
    def app1(idx, pt):
        m,n = len(pt), len(idx)      
        mask = np.zeros((m, n), dtype='bool')
        return array_masking(mask, idx, pt)
    

    时间安排 -

    In [5]: np.random.seed(0)
       ...: x = np.random.rand(1000000)
       ...: y = np.random.rand(200)
    
    In [6]: %timeit app1(idx, pt)
    264 ms ± 8.91 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    
    In [7]: %timeit app2(idx, pt)
    165 ms ± 3.43 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
    

    【讨论】:

    • 谢谢。这是对性能的提升。你的代码在我电脑上的运行速度比我之前的快 30%。
    • @f.c.方法#1对您来说是一个可行的解决方案吗?你能给我们一些关于你倾向于获得赏金的解决方案的反馈吗?如果需要可以改进答案,可能会有所帮助。
    • 抱歉回复晚了。我已经测试了你的方法#1。到目前为止,你是赢家。谢谢!
    • 顺便说一句,由于内存问题,我将scale_factor增加到2000,在我的计算机上加速大约是2.5倍,非常棒。
    • @f.c.惊人的!感谢您的反馈。
    【解决方案2】:

    这是一个替代答案,但不确定它是否正是您所需要的。

    x = np.random.rand(1000000)
    y = np.random.rand(200)
    mask = x[np.newaxis,:] < y[:,np.newaxis]
    

    注意:我提到它可能不是您需要的,因为您指定了 numpy.searchsorted() 的需要,在这里我不使用它,但是我得到了相同的结果。如果它不能完全满足您的需求,它也可能在将来对其他人有用;)。

    时间安排(@DanielF 编辑)

    设置:

    import numpy as np
    
    # generate test examples
    x = np.random.rand(1000000)
    y = np.random.rand(200)
    
    # sort x
    idx = np.argsort(x)
    sorted_x = np.take_along_axis(x, idx, axis=-1)
    

    跑步:

    %%timeit   #  silgon
    mask = x[np.newaxis,:] < y[:,np.newaxis]
    
    166 ms ± 3.99 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
    
    
    %%timeit    # Divakar
    pt = np.searchsorted(sorted_x, y)
    mask = app2(idx, pt)
    
    316 ms ± 29 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    
    
    %%timeit   #  f.c.
    pt = np.searchsorted(sorted_x, y)
    mask = np.zeros((200, 1000000), dtype='bool')
    for i in range(200):
         mask[i, idx[0:pt[i]]] = True
         
    466 ms ± 13.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    

    【讨论】:

    • 你是对的。它是提供相同结果的蛮力搜索版本。使用 numpy.searchsorted 的原因是为了加快计算速度。
    • @f.c.与此相比,np.searchsorted 真的可以加快计算速度吗?
    • @DanielF 好问题,是的。它在我的计算机上节省了 70% 的时间。
    • 用我的时间编辑了答案。
    • 在那之前我已经运行编译它了。在那之前是~800ms
    【解决方案3】:

    我已经实现了一个 numba 版本的 for 循环,它比 python 版本稍微快一点。代码如下:

    import numba as nb
    @nb.njit(parallel=True,fastmath=True)
    def array_masking(mask, idx, pt):
        for j in nb.prange(pt.shape[0]):
            for i in nb.prange(pt[j]):
                mask[j, idx[i]] = True
    

    我正在寻求进一步的加速。有什么想法吗?

    【讨论】:

      【解决方案4】:

      为了更清楚起见,您正在寻找一种快速方法来确定小于 y[i]x 项的索引掩码。例如,如果对 x 的项目进行排序的索引是:

      np.argsort(x) = [5, 0, 2, 10, 7, 8, 9, 11, 1, 6, 3, 4]
      

      并且您知道 8 的项目小于 y[i],然后您需要从该列表的倒序中选择前 8 个项目:

      arg_inv = [1, 8, 2, 10, 11, 0, 9, 4, 5, 6, 3, 7]
      

      这个问题最简单的方法是高级索引:

      length_x, length_y = len(x), len(y)
      idx = np.argsort(x)
      arg_inv = np.argsort(idx)
      pt = np.searchsorted(x, y, sorter=idx)
      mask = np.zeros((length_y, length_x), dtype='bool')
      row, col = np.divmod(np.arange(length_x * length_y), length_x)
      mask[row, col] = arg_inv[col] < pt[row]
      return mask
      

      我还加了一个小样本的例子:

      x = [0.809 0.958 0.881 0.146 0.882 0.421 0.604]
      y = [0.119 0.981 0.775 0.254]
      np.sort(x) = [0.146 0.421 0.604 0.809 0.881 0.882 0.958]
      np.argsort(x) = [3 5 6 0 2 4 1]
      arg_inv = [3 6 4 0 5 1 2]
      pt = [0 7 3 1]
      Process of advanced indexing:
      
          row  col  arg_inv[col]  pt[row]  arg_inv[col] < pt[row]
      0     0    0             3        0                       0
      1     0    1             6        0                       0
      2     0    2             4        0                       0
      3     0    3             0        0                       0
      4     0    4             5        0                       0
      5     0    5             1        0                       0
      6     0    6             2        0                       0
      7     1    0             3        7                       1
      8     1    1             6        7                       1
      9     1    2             4        7                       1
      10    1    3             0        7                       1
      11    1    4             5        7                       1
      12    1    5             1        7                       1
      13    1    6             2        7                       1
      14    2    0             3        3                       0
      15    2    1             6        3                       0
      16    2    2             4        3                       0
      17    2    3             0        3                       1
      18    2    4             5        3                       0
      19    2    5             1        3                       1
      20    2    6             2        3                       1
      21    3    0             3        1                       0
      22    3    1             6        1                       0
      23    3    2             4        1                       0
      24    3    3             0        1                       1
      25    3    4             5        1                       0
      26    3    5             1        1                       0
      27    3    6             2        1                       0
      

      【讨论】:

      • 如果x有重复的元素似乎不起作用,例如x = np.array([0.1 0.3 0.2 0.25 0.2 0.25])
      • @f.c.感谢您的回复。我认为np.random.rand 通话是不可能的。稍后我会解决这个问题,同时,arg_inv = np.argsort(np.argsort(x)) 而不是np.unique 可能有帮助吗?
      • 我使用np.random.rand 来生成测试数据。实际上,x 确实有重复值。谢谢。
      • 我认为 np.argsort(np.argsort(x)) 应该可以工作。我会测试性能。
      • 不幸的是,您的屏蔽比 for 循环慢 7 倍,并且消耗的内存多 4 倍。我会保留for循环。但是感谢您的尝试。
      猜你喜欢
      • 1970-01-01
      • 2013-02-14
      • 2021-10-16
      • 2011-08-11
      • 1970-01-01
      • 1970-01-01
      • 2017-08-23
      • 2017-06-18
      • 1970-01-01
      相关资源
      最近更新 更多