如何对这个大型数组计算进行矢量化和加速？

Question

我目前正在尝试计算10.000 x 10.000数组值中所有子方格总和的总和.例如,如果我的数组是:

1 1 1 
2 2 2
3 3 3 

我希望结果如下:

1+1+1+2+2+2+3+3+3                        [sum of squares of size 1]
+(1+1+2+2)+(1+1+2+2)+(2+2+3+3)+(2+2+3+3) [sum of squares of size 2]
+(1+1+1+2+2+2+3+3+3)                     [sum of squares of size 3]
________________________________________
68

所以,作为第一次尝试,我写了一个非常简单的python代码来做到这一点.因为它在O(k ^ 2.n ^ 2)中(n是大数组的大小,k是我们得到的子方形的大小),处理非常长.我在O(n ^ 2)中写了另一个算法来加速它:

def getSum(tab,size):
    n = len(tab)
    tmp = numpy.zeros((n,n))

    for i in xrange(0,n):
        sum = 0
        for j in xrange(0,size):
            sum += tab[j][i]
        tmp[0][i] = sum

        for j in xrange(1,n-size+1):
            sum += (tab[j+size-1][i] - tab[j-1][i])
            tmp[j][i] = sum

    finalsum = 0
    for i in xrange(0,n-size+1):
        sum = 0 
        for j in xrange(0,size):
            sum += tmp[i][j]
        finalsum += sum

        for j in xrange(1,n-size+1):
            finalsum += (tmp[i][j+size-1] - tmp[i][j-1])

return finalsum

所以这段代码工作正常.给定一个数组和子方格的大小,它将返回所有这些子方格中的值的总和.我基本上迭代了子方块的大小来获得所有可能的值.

问题是,对于大型阵列来说,这又是一件好事(对于10.000 x 10.000阵列,超过20天).我用Google搜索并了解到我可以使用numpy对迭代进行矢量化.但是,在我的情况下,我无法弄清楚如何制作它......

如果有人可以帮助我加快算法速度,或者给我一些关于这个主题的好文档,我会很高兴的!

谢谢 !

Answer 1

根据@Divakar的优秀想法,我建议使用积分图像来加速卷积.如果矩阵非常大,则必须多次卷积(每个内核大小一次).使用积分图像(也称为求和区域表)可以非常有效地计算几个卷积(或正方形内的和的评估).

一旦积分图像M计算,区域内的所有值的总和(x0, y0) - (x1, y1)可以与计算只是 4 aritmetic计算,窗口的大小无关(图片来自维基百科):

M[x1, y1] - M[x1, y0] - M[x0, y1] + M[x0, y0]

这可以很容易地在numpy中进行矢量化.可以使用计算积分图像cumsum.以下示例:

tab = np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]])
M = tab.cumsum(0).cumsum(1) # Create integral images
M = np.pad(M, ((1,0), (1,0)), mode='constant') # pad it with a row and column of zeros

M用一行和一列零填充以处理第一行(其中x0 = 0或y0 = 0).

然后,给定窗口大小W,W可以有效地计算每个窗口大小的总和,并使用numpy完全向量化:

all_sums = M[W:, W:] - M[:-W, W:] - M[W:, :-W] + M[:-W, :-W]

注意,上面的矢量化操作计算每个窗口的总和,即矩阵的每个A,B,C和D. 然后计算所有窗口的总和

total = all_sums.sum()

请注意,对于N不同大小,不同于卷积,整数图像只需计算一次,因此,代码可以非常有效地编写为:

def get_all_sums(A):
    M = A.cumsum(0).cumsum(1)
    M = np.pad(M, ((1,0), (1,0)), mode='constant')

    total = 0
    for W in range(1, A.shape[0] + 1):
        tmp = M[W:, W:] + M[:-W, :-W] - M[:-W, W:] - M[W:, :-W]
        total += tmp.sum()

    return total

示例的输出:

>>> get_all_sums(tab)
68

一些时序将卷积与具有不同大小矩阵的积分图像进行比较.getAllSums对于Divakar的卷积方法,以及get_all_sums上述基于积分图像的方法:

>>> R1 = np.random.randn(10, 10)
>>> R2 = np.random.randn(100, 100)

1)使用R110x10矩阵:

>>> %time getAllSums(R1)
CPU times: user 353 µs, sys: 9 µs, total: 362 µs
Wall time: 335 µs
2393.5912717342017

>>> %time get_all_sums(R1)
CPU times: user 243 µs, sys: 0 ns, total: 243 µs
Wall time: 248 µs
2393.5912717342012

2)使用R2100x100矩阵:

>>> %time getAllSums(R2)
CPU times: user 698 ms, sys: 0 ns, total: 698 ms
Wall time: 701 ms
176299803.29826894

>>> %time get_all_sums(R2)
CPU times: user 2.51 ms, sys: 0 ns, total: 2.51 ms
Wall time: 2.47 ms
176299803.29826882

请注意,对于足够大的矩阵,使用积分图像比卷积快300倍.

Answer 2

这些滑动求和最适合计算为 2D 卷积求和，并且可以使用 有效计算scipy's convolve2d。因此，对于特定的大小，您可以获得总和，如下所示 -

def getSum(tab,size):
    # Define kernel and perform convolution to get such sliding windowed summations
    kernel = np.ones((size,size),dtype=tab.dtype)
    return convolve2d(tab, kernel, mode='valid').sum()

为了获得所有大小的总和，我认为在内存和性能效率方面最好的方法是使用循环来循环所有可能的大小。因此，为了得到最终的总和，你需要 -

def getAllSums(tab):
    finalSum = 0
    for i in range(tab.shape[0]):
        finalSum += getSum(tab,i+1)
    return finalSum

样本运行 -

In [51]: tab
Out[51]: 
array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

In [52]: getSum(tab,1) # sum of squares of size 1
Out[52]: 18

In [53]: getSum(tab,2) # sum of squares of size 2
Out[53]: 32

In [54]: getSum(tab,3) # sum of squares of size 3
Out[54]: 18

In [55]: getAllSums(tab) # sum of squares of all sizes
Out[55]: 68