Ускорение вычислений с матрицами numpy

Question

У меня есть две матрицы. Оба заполнены нулями и единицами. Один из них большой (3000 x 2000 элементов), а другой - меньшие (20 x 20) элементов. Я делаю что-то наподобие:

newMatrix = (size of bigMatrix), filled with zeros 
l = (a constant) 

for y in xrange(0, len(bigMatrix[0])): 
    for x in xrange(0, len(bigMatrix)): 

     for b in xrange(0, len(smallMatrix[0])): 
      for a in xrange(0, len(smallMatrix)): 

       if (bigMatrix[x, y] == smallMatrix[x + a - l, y + b - l]): 
        newMatrix[x, y] = 1 

Который является болезненно медленным. Я что-то делаю неправильно? Есть ли умный способ сделать эту работу быстрее?

Редактирование: В основном я для каждой (x, y) в большой матрице, проверяя все пиксели как большой матрицы, так и малой матрицы вокруг (x, y), чтобы увидеть, являются ли они 1. Если они 1, то я устанавливаю это значение на newMatrix. Я делаю своеобразное обнаружение столкновений.

Answer 1

Я могу придумать пару оптимизаций там - Поскольку вы используете 4 вложенных оператора python for for, вы примерно так же медленны, как и вы.

Я не могу точно определить, что вы ищете - , но с одной стороны, если плотность вашей большой матрицы «1» низкая, вы можете, конечно, использовать «любую» функцию python на срезах bigMtarix, чтобы быстро проверить если есть какие-либо набор элементов там - вы могли бы получить увеличивается в несколько раз скорость там:

step = len(smallMatrix[0]) 
for y in xrange(0, len(bigMatrix[0], step)): 
    for x in xrange(0, len(bigMatrix), step): 
     if not any(bigMatrix[x: x+step, y: y + step]): 
      continue 
     (...) 

на данный момент, если по-прежнему должны взаимодействовать на каждом элементе, вы делаете еще одну пару индексов, чтобы ходить каждую позицию внутри шага - но я думаю, что у вас есть идея.

Помимо использования внутренних цифровых операций, таких как «любое» использование, вы могли бы, конечно, добавить некоторый код потока управления, чтобы отключить цикл (b, a) при обнаружении первого совпадающего пикселя. (Например, вставив инструкцию «break» внутри вашей последней «if» и другой пары if..break для цикла «b».

Я действительно не могу понять, в чем вы намерены, поэтому я могу 't дать вам более конкретный код.

Answer 2

Ваш пример кода не имеет смысла, но описание вашей проблемы звучит так, будто вы пытаетесь выполнить 2-мерную свертку небольшого bitarray над большим bitarray. В scipy есть функция convolve2d .signal, который делает именно это. Просто сделайте convolve2d(bigMatrix, smallMatrix), чтобы получить результат. К сожалению, у scipy-реализации нет специального случая для булевых массивов, поэтому полная свертка довольно медленная. Вот функция, которая использует тот факт, что массивы содержат только единицы и нули:

import numpy as np 

def sparse_convolve_of_bools(a, b): 
    if a.size < b.size: 
     a, b = b, a 
    offsets = zip(*np.nonzero(b)) 
    n = len(offsets) 
    dtype = np.byte if n < 128 else np.short if n < 32768 else np.int 
    result = np.zeros(np.array(a.shape) + b.shape - (1,1), dtype=dtype) 
    for o in offsets: 
     result[o[0]:o[0] + a.shape[0], o[1]:o[1] + a.shape[1]] += a 
    return result 

На моей машине он работает менее чем за 9 секунд для свертки 3000x2000 на 20x20. Время работы зависит от количества единиц в меньшем массиве, составляя 20 мс на каждый ненулевой элемент.

Answer 3

Если биты действительно упакованы 8 на байты/32 за междунар, и вы можете уменьшить smallMatrix до 20x16,
попробуйте следующее, здесь для одной строки.
(newMatrix[x, y] = 1 когда любой бит 20x16 вокруг х, у 1 ?? Что вы действительно ищете?)

python -m timeit -s ' 
""" slide 16-bit mask across 32-bit pairs bits[j], bits[j+1] """ 

import numpy as np 

bits = np.zeros(2000 // 16, np.uint16) # 2000 bits 
bits[::8] = 1 
mask = 32+16 
nhit = 16 * [0] 

def hit16(bits, mask, nhit): 
    """ 
     slide 16-bit mask across 32-bit pairs bits[j], bits[j+1] 
     bits: long np.array(uint16) 
     mask: 16 bits, int 
     out: nhit[j] += 1 where pair & mask != 0 
    """ 
    left = bits[0] 
    for b in bits[1:]: 
     pair = (left << 16) | b 
     if pair: # np idiom for non-0 words ? 
      m = mask 
      for j in range(16): 
       if pair & m: 
        nhit[j] += 1 
        # hitposition = jb*16 + j 
       m <<= 1 
     left = b 
    # if any(nhit): print "hit16:", nhit 

' \ 
' 
hit16(bits, mask, nhit) 
' 

# 15 msec per loop, bits[::4] = 1 
# 11 msec per loop, bits[::8] = 1 
# mac g4 ppc