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

Question

У меня есть две матрицы. Оба заполнены нулями и единицами. Один большой (3000 х 2000 элементов), а другой - меньше (20 х 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

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

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

Answer 1

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

Я не могу точно понять, что вы ищете - но, с одной стороны, если плотность вашей большой матрицы "1" низкая, вы, безусловно, можете использовать python's "any" функция на срезах 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 ».

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

Answer 2

Ваш пример кода не имеет смысла, но описание вашей проблемы звучит так, как будто вы пытаетесь выполнить двумерную свертку небольшого битрейра над большим битаррей. В пакете scipy.signal есть функция convolve2d, которая делает именно это. Просто сделайте 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 на int, и вы можете уменьшить ваш smallMatrix до 20x16,
затем попробуйте следующее здесь для одной строки.
( newMatrix [x, y] = 1 , когда любой бит 20x16 вокруг x, y равен 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