将嵌套循环计算转换为 Numpy 以加速

我的Python程序的一部分包含以下代码段，其中一个新的网格 是根据旧网格中找到的数据计算的。

网格是二维浮点数列表。该代码使用了三个 for 循环：

for t in xrange(0, t, step):
    for h in xrange(1, height-1):
        for w in xrange(1, width-1):
            new_gr[h][w] = gr[h][w] + gr[h][w-1] + gr[h-1][w] + t * gr[h+1][w-1]-2 * (gr[h][w-1] + t * gr[h-1][w])
    gr = new_gr

return gr

对于大网格和大量时间来说，代码非常慢t.

我尝试使用 Numpy 通过替换内部循环来加速此代码 和：

J = np.arange(1, width-1)
new_gr[h][J] = gr[h][J] + gr[h][J-1] ...

但产生的结果（数组中的浮点数）大约比 他们的列表计算对应物。

• 使用以下方法将浮点数列表转换为 Numpy 浮点数数组时，预计会出现什么精度损失np.array(pylist)然后进行计算？

• 我应该如何将三重 for 循环转换为漂亮且快速的 Numpy 代码？ （或者还有其他显着加快代码速度的建议吗？）

If gr是浮点列表，如果您希望使用 NumPy 进行矢量化，第一步是转换gr到 NumPy 数组np.array() http://docs.scipy.org/doc/numpy/reference/generated/numpy.array.html.

接下来，我假设你有new_gr用形状的零点初始化(height,width)。在最里面的两个循环中执行的计算基本上表示2D convolution。所以，你可以使用signal.convolve2d http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.signal.convolve2d.html具有适当的kernel。决定kernel，我们需要查看缩放因子并做出3 x 3将它们从内核中取出并对它们取反以模拟我们在每次迭代中所做的计算。因此，您将拥有一个矢量化解决方案，其中两个最里面的循环被删除以获得更好的性能，如下所示 -

import numpy as np
from scipy import signal

# Get the scaling factors and negate them to get kernel
kernel = -np.array([[0,1-2*t,0],[-1,1,0,],[t,0,0]])

# Initialize output array and run 2D convolution and set values into it
out = np.zeros((height,width))
out[1:-1,1:-1] = signal.convolve2d(gr, kernel, mode='same')[1:-1,:-2]

验证输出和运行时测试

定义函数：

def org_app(gr,t):
    new_gr = np.zeros((height,width))
    for h in xrange(1, height-1):
        for w in xrange(1, width-1):
            new_gr[h][w] = gr[h][w] + gr[h][w-1] + gr[h-1][w] + t * gr[h+1][w-1]-2 * (gr[h][w-1] + t * gr[h-1][w]) 
    return new_gr

def proposed_app(gr,t):
    kernel = -np.array([[0,1-2*t,0],[-1,1,0,],[t,0,0]])
    out = np.zeros((height,width))
    out[1:-1,1:-1] = signal.convolve2d(gr, kernel, mode='same')[1:-1,:-2]
    return out

Verify -

In [244]: # Inputs
     ...: gr = np.random.rand(40,50)
     ...: height,width = gr.shape
     ...: t = 1
     ...: 

In [245]: np.allclose(org_app(gr,t),proposed_app(gr,t))
Out[245]: True

Timings -

In [246]: # Inputs
     ...: gr = np.random.rand(400,500)
     ...: height,width = gr.shape
     ...: t = 1
     ...: 

In [247]: %timeit org_app(gr,t)
1 loops, best of 3: 2.13 s per loop

In [248]: %timeit proposed_app(gr,t)
10 loops, best of 3: 19.4 ms per loop