使用 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”的密度很低，你当然可以在 bigMtarix 的切片上使用 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
        (...) 

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

除了使用像“任何”用法这样的内部数字操作之外，您当然可以添加一些控制流代码来在找到第一个匹配像素时中断 (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 x 20x20 卷积的时间不到 9 秒。运行时间取决于较小数组中的 1 数量，每个非零元素为 20 毫秒。

Answer 3

如果你的位确实被打包为每个字节 8 个/每个整数 32 个，
你可以将你的smallMatrix减少到20x16，
然后尝试以下操作，此处针对单行。
（newMatrix[x, y] = 1，当 x,y 周围 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