Numpy大步招数抱怨"阵列太大",为什么？

Question

在numpy(1.8)中,我想将这个计算从Python循环中转移到更多numpy-ish以获得更好的性能:

(width, height) = base.shape
(toolw, toolh) = tool.shape
for i in range(0, width-toolw):
    for j in range(0, height-toolh):
        zdiff[i,j] = (tool - base[i:i+toolw, j:j+toolh]).min()

base是一个~2000x2000阵列,tool是一个25x25阵列.(背景背景:基础和工具是高度图,我试图找出最接近基础的工具移动方法.)

我试图使用一个跨越式的技巧,从这开始:

base_view = np.lib.stride_tricks.as_strided(base, shape=(2000, 2000, 25, 25), 
                                            strides=(base.strides * 2))

这将base_view[10,20]是从(10,20)的左上角的基数开始的25x25值的数组.

然而,这与"阵列太大"失败了.从价值测试来看,当数组的潜在大小(例如2000*2000*25*25*8)超过2 ^ 32-ish并且触发溢出检查时,它会报告此问题,从而将所有维度相乘.(我使用的是32位Python安装).

我觉得我错过了一些东西 - 为什么当步幅值明显起作用时,它不会让我创造这种"跨步观点"？有没有办法强迫这个？

更一般地说,有没有办法优化我的循环？

更新:确切错误:

ValueError                                Traceback (most recent call last)
<ipython-input-14-313b3d6c74fa> in <module>()
----> 1 newa = np.lib.stride_tricks.as_strided(base, shape=(1000, 1000, 25, 25), strides=(base.strides * 2))

C:\Python27\lib\site-packages\numpy\lib\stride_tricks.pyc in as_strided(x, shape, strides)
     28     if strides is not None:
     29         interface['strides'] = tuple(strides)
---> 30     array = np.asarray(DummyArray(interface, base=x))
     31     # Make sure dtype is correct in case of custom dtype
     32     array.dtype = x.dtype

C:\Python27\lib\site-packages\numpy\core\numeric.pyc in asarray(a, dtype, order)
    458 
    459     """
--> 460     return array(a, dtype, copy=False, order=order)
    461 
    462 def asanyarray(a, dtype=None, order=None):

ValueError: array is too big.

Answer 1

我无法真正帮助您使用跨步方法，但确实有一种比您的原始代码更快的方法。它在tool数组上循环而不是在base数组上循环，这意味着，虽然没有完全矢量化，但更多的工作被推给了 numpy。

请注意，在您的原始代码中，我更改了范围并切换了宽度和高度，因为我认为这就是您的意图。

import numpy as np

height, width = 500, 500
toolh, toolw = 6, 6

base = np.random.rand(height, width)
tool = np.random.rand(toolh, toolw)

m, n = height-toolh+1, width-toolw+1

def height_diff_old(base, tool):
    zdiff = np.empty((m, n))
    for i in range(m):
        for j in range(n):
            zdiff[i, j] = (tool - base[i:i+toolh, j:j+toolw]).min()
    return zdiff

def height_diff_new(base, tool):
    zdiff = np.empty((m, n))
    zdiff.fill(np.inf)
    for i in range(toolh):
        for j in range(toolw):
            diff_ij = tool[i, j] - base[i:i+m, j:j+n]
            np.minimum(zdiff, diff_ij, out=zdiff)
    return zdiff

当然，您希望计算实际函数中的高度和宽度，但为了测试，将它们作为全局变量会更容易。

对于给定的数组大小，原始代码在 7.38 秒内运行，而新代码在我的系统上只需要 206 毫秒。我认为新代码对于您的数组大小也更快，但我不确定它的扩展效果如何:)

您可能感兴趣或不感兴趣的其他替代方案是使用 Numba 或 Cython，它们在许多情况下应该比您想到的任何“矢量化”numpy 代码更快。