使用numpy矩阵加速计算

Question

我有两个矩阵。两者都用0和1填充。一个是大的(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“函数来快速检查其中是否有任何集合元素--你可以在那里获得几倍的速度：

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
        (...) 

此时，如果仍然需要对每个元素进行交互，您可以执行另一对索引来遍历步骤中的每个位置-但我认为您已经明白了这一点。

除了使用像这样的"any“用法的内部数值操作之外，当找到第一个匹配的像素时，您当然可以添加一些控制流代码来中断(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

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

Answer 3

如果你的位真的是每字节8位/每整型32位，并且你可以将你的smallMatrix减少到20x16，

然后尝试执行以下操作，这里只有一行。

(当x，y周围的20x16的任何一位是1时，newMatrix[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