Блог

допустим, у меня есть 2D массив, т. е.:

 In [136]: ary          
array([[6, 7, 9],
       [0, 2, 5],
       [3, 3, 4],
       [2, 2, 8],
       [3, 4, 9],
       [0, 5, 7],
       [2, 4, 9],
       [3, 5, 7],
       [7, 8, 8],
       [0, 2, 3]])
 

Я хочу быстро вычислить перекрытие с 1D вектором.

Я почти могу сделать это с помощью (8 мс на большом массиве):

  (ary == match)     # .sum(axis=1).argsort()[::-1]
 

Проблема в том, что он совпадает только в том случае, если совпадают и позиция, и значение.

 match == [6,5,4]                                                                                                                                                      
 
array([[ True, False, False],
       [False, False, False],
       [False, False,  True],
       [False, False, False],
       [False, False, False],
       [False,  True, False],
       [False, False, False],
       [False,  True, False],
       [False, False, False],
       [False, False, False]])
 

Например, 5 во 2-м столбце 1d vec не совпало с 5 в 3-м столбце во 2-й строке.

Он работает с .isin()

   np.isin(ary,match,assume_unique=True).sum(axis=1).argsort()[::-1][:5]
 

но это медленно на больших массивах (200000,10) ~20 мс

Помогите мне расширить первый случай, чтобы он мог сопоставить значение в любой позиции вектора 1D со строкой.

ожидаемый результат-индексы строк, упорядоченные по КОЛИЧЕСТВУ ПЕРЕКРЫТИЙ, позволяет использовать [2,4,5], потому что в нем больше совпадений:

 In [147]: np.isin(ary,[2,5,4],assume_unique=True)                                                                                                                             
Out[147]: 
array([[False, False, False],
       [False,  True,  True],
       [False, False,  True],
       [ True,  True, False],
       [False,  True, False],
       [False,  True, False],
       [ True,  True, False],
       [False,  True, False],
       [False, False, False],
       [False,  True, False]])
 

Перекрытие :

 In [149]: np.isin(ary,[2,5,4],assume_unique=True).sum(axis=1)                                                                                                                 
Out[149]: array([0, 2, 1, 2, 1, 1, 2, 1, 0, 1])
 

Порядок перекрытия :

 In [148]: np.isin(ary,[2,5,4],assume_unique=True).sum(axis=1).argsort()[::-1]                                                                                                 
Out[148]: array([6, 3, 1, 9, 7, 5, 4, 2, 8, 0])
 

Смотрите строки : 6,3,1 имеют перекрытие:2 вот почему они первые

Варианты:

 #could be from 1000,10,10 to 2000,100,20 ..    
def data(cells=2000,seg=100,items=10): 
    ary = np.random.randint(0,cells,(cells*seg,items))
    rnd = np.random.randint(0,cells*seg)
    return ary, ary[rnd]

def best2(match,ary): #~20ms, (200000,10)
    return np.isin(ary,match,assume_unique=True).sum(axis=1).argsort()[::-1][:5]                                                                                                     

def best3(match,ary): #Corralien  ~20ms
    return np.logical_or.reduce(np.ravel(ary) == match[:, None], axis=0).reshape(ary.shape).sum(axis=1).argsort()[::-1][:5]
 

Можно ли ускорить это, если использовать numba cuda ИЛИ cupy на GPU ?

Основная проблема всех подходов так быстро заключается в том, что они создают огромный временный массив, в то время как, наконец, важны только 5 элементов. Numba можно использовать для вычисления массивов на лету (с эффективными циклами JIT-компиляции), избегая некоторых временных массивов. Более того, полная сортировка не требуется, так как необходимо извлечь только 5 лучших элементов. Вместо этого можно использовать раздел. Можно даже использовать более быстрый подход, поскольку важны только 5 выбранных элементов, а не остальные. Вот полученный код:

 @nb.njit('int32[::1](int32[::1], int32[:,::1])')
def computeScore(match, ary):
    n, m = ary.shape
    assert m == match.shape[0]
    tmp = np.empty(n, dtype=np.int32)
    for i in range(n):
        s = 0
        # Count the number of matching items (with repetition)
        for j in range(m):
            # Find a match
            item = ary[i, j]
            found = False
            for k in range(m):
                found |= item == match[k]
            s  = found
        tmp[i] = s
    return tmp

def best4(match, ary):
    n, m = ary.shape
    score = computeScore(match, ary)
    bestItems = np.argpartition(score, n-5)[n-5:] # sadly not supported by Numba yet
    order = np.argsort(-score[bestItems]) # bastItems is not sorted and likely needs to be
    return bestItems[order]
 

Обратите внимание, что это best4 может привести к результатам, отличным best2 от тех, когда соответствующий балл (сохраненный в tmp ) равен между несколькими элементами. Это связано с алгоритмом сортировки, который по умолчанию нестабилен в Numpy ( kind параметр можно использовать для адаптации этого поведения). Это также верно для алгоритма разбиения, хотя Numpy, похоже, еще не обеспечивает стабильный алгоритм разбиения.

Этот код должен быть быстрее, чем другие реализации, но не с большим отрывом. Одна из проблем заключается в том, что Numba (и большинству компиляторов C/C , подобных тому, который используется для компиляции Numpy) не удается создать быстрый код, поскольку он не знает значения m во время компиляции. В результате наиболее агрессивные оптимизации (например, развертывание циклов и использование инструкций SIMD) вряд ли могут быть применены. Вы можете помочь Numba, используя утверждения или избегая условных выражений.

Более того, код может быть распараллелен с использованием нескольких потоков, чтобы сделать его намного быстрее на основных платформах. Обратите внимание, что распараллеленная версия может работать не быстрее на небольших данных и не на всех платформах, поскольку создание потоков приводит к накладным расходам, которые могут быть больше, чем фактические вычисления.

Вот результат реализации:

 @nb.njit('int32[::1](int32[::1], int32[:,::1])', parallel=True)
def computeScoreOpt(match, ary):
    n, m = ary.shape
    assert m == match.shape[0]
    assert m == 10
    tmp = np.empty(n, dtype=np.int32)
    for i in nb.prange(n):
        # Thie enable Numba to assume m=10 in the following code
        # and generate a very efficient code for this specific case.
        # The assert should be enough but the internals of Numba 
        # prevent the information to be propagatted to this portion
        # of the code when it is parallelized.
        if m != 10: continue

        s = 0
        for j in range(m):
            item = ary[i, j]
            found = False
            for k in range(m):
                found |= item == match[k]
            s  = found
        tmp[i] = s
    return tmp

def best5(match, ary):
    n, m = ary.shape
    score = computeScoreOpt(match, ary)
    bestItems = np.argpartition(score, n-5)[n-5:]
    order = np.argsort(-score[bestItems])
    return bestItems[order]
 

Вот тайминги на моей машине с примером набора данных:

 best2:                            18.2 ms
best3:                            17.8 ms
best4 (sequential -- default):    12.0 ms
best4 (parallel):                  3.1 ms
best5 (sequential):                3.2 ms
best5 (parallel -- default):       1.2 ms
 

Самая быстрая реализация в 15 раз быстрее, чем исходная эталонная реализация.

Обратите внимание, что если m значение больше примерно 30, то лучше использовать более продвинутый алгоритм, основанный на множествах. Альтернативным решением match в этом случае является сначала сортировка, а затем использование np.isin в цикле на i основе.

Использовать вещание и np.logical_or.reduce :

 # match = np.array(match)
>>> np.logical_or.reduce(np.ravel(ary) == match[:, None], axis=0) 
                 .reshape(ary.shape)

array([[ True, False, False],
       [False, False,  True],
       [False, False,  True],
       [False, False, False],
       [False,  True, False],
       [False,  True, False],
       [False,  True, False],
       [False,  True, False],
       [False, False, False],
       [False, False, False]])
 

Производительность

 match = np.array([6, 5, 4])
ary = np.random.randint(0, 10, (200000, 10))
 

 %timeit np.logical_or.reduce(np.ravel(ary) == match[:, None], axis=0).reshape(ary.shape)
7.49 ms ± 174 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)