使用numpy矩阵加速计算

Question

我有两个矩阵。一个是大矩阵，里面全是0和1，大小是3000行乘2000列；另一个是小矩阵，大小是20行乘20列，同样也是全0和1。我正在做一些类似这样的操作：

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)位置，检查大矩阵和小矩阵周围的所有像素，看看它们是否是1。如果是1，我就把这个值放到一个新的矩阵中。我在做一种碰撞检测。

Answer 1

如果你的数据是每个字节存8个比特，每个整数存32个比特，并且你可以把你的小矩阵缩小到20x16的大小，
那么可以试试下面的方法，这里是针对单行的情况。
(当20x16区域内的x,y位置周围有任意一个比特是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

Answer 2

你的示例代码看起来有点不太对劲，但你描述的问题听起来像是你想在一个大数组上进行一个小数组的二维卷积。其实在scipy.signal这个包里有一个叫做convolve2d的函数，正好可以做到这一点。你只需要用 convolve2d(bigMatrix, smallMatrix) 就能得到结果。不过可惜的是，scipy的这个实现对布尔数组没有特别的处理，所以完整的卷积过程会比较慢。下面有一个函数可以利用数组里只有0和1的特点：

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秒。运行时间还会根据小数组里1的数量而变化，每个非零元素大约需要20毫秒。

Answer 3

我想到了一些优化的方法——

你现在用的是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”的组合。）

我真的不太明白你具体想要做什么，所以无法给你更具体的代码。