用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”时陈述，你的速度和你一样慢。

我无法确切地知道你在寻找什么 - 但是一方面，如果你的大矩阵“1”的密度很低，你肯定可以使用python的“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）循环。 （比如，在最后一个“if”中插入一个“break”语句，并为“b”循环插入另一个if..break对。

我真的无法确切地知道你的意图是什么 - 所以我不能给你更多的具体代码。

Answer 2

您的示例代码毫无意义，但您对问题的描述听起来像是在尝试对大比特阵列进行小比特的二维卷积。在scipy.signal包中有一个convolve2d函数就是这样做的。只需执行 convolve2d（bigMatrix，smallMatrix）即可获得结果。不幸的是，scipy实现没有布尔数组的特殊情况，因此完全卷积相当慢。这是一个利用数组只包含1和0的事实的函数：

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

在我的机器上，对于3000x2000乘20x20卷积，它在不到9秒的时间内运行。运行时间取决于较小数组中的数量，每个非零元素为20ms。

Answer 3

如果你的位真正打包每个字节8个/每个int 32个， 你可以将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