我可以避免使用 numpy 进行动态编程的 Python 循环开销吗？

Question

我需要帮助解决以下问题的 Pythonic 循环开销：我正在编写一个计算像素流算法的函数，该算法是 2D Numpy 数组上的经典动态编程算法。它需要：

1) 像这样访问数组的所有元素至少一次：

for x in range(xsize):
    for y in range(ysize):
         updateDistance(x,y)

2) 可能会根据看起来像这样的元素的邻居的值遵循元素路径

while len(workingList) > 0:
   x,y = workingList.pop()
   #if any neighbors of x,y need calculation, push x,y and neighbors on workingList 
   #else, calculate flow on pixels as a sum of flow on neighboring pixels

不幸的是，即使对 updateDistance 的调用是 pass，我似乎在 #1 上也得到了很多 Pythonic 循环开销。我认为这是一个足够经典的算法，在 Python 中必须有一个好的方法可以避免一些循环开销。我还担心如果我能修复 #1，我会在 #2 的循环中被破坏。

关于快速遍历二维 numpy 数组中的元素并可能更新该数组中的元素链有什么建议吗？

编辑：清除 #2 的更多细节

似乎我可以对第一个循环进行矢量化，也许可以通过矢量化 np.meshgrid 调用。

部分循环有点复杂，但这里有一个简化的版本。我担心循环和对相邻元素的索引：

#A is a 2d cost matrix
workingList = [(x,y)]
while len(workingList) > 0:
   x,y = workingList.pop()
   neighborsToCalculate = []
   for n in neighborsThatNeedCalculation(x,y): #indexes A to check neighbors of (x,y)
      neighborsToCalculate.append(n)
   if len(neighborstToCalculate) != 0:
       workingList.append((x,y))
       workingList.extend(neighborsToCalculate)
   else:
       for xn,yn in neighbors(x,y):
          A[x,y] += 1+A[xn,yn]

这是一个经典的广度优先搜索问题。如果这可以并行化，那就太好了。它可能不能以目前的形式出现，因为它遵循一条路径，但我很乐意提供建议。

Answer 1

如果您在算法中使用 python 循环，您不会从 numpy 获得任何速度提升。您需要并行化您的问题。

在图像处理中，并行化意味着对所有像素使用相同的函数，即使用内核。在 numpy 中，而不是这样做：

for x in range(xsize):
    for y in range(ysize):
         img1[y, x] = img2[y, x] + img3[y, x]

你会的：

img1 = img2 + img3 # add 2 images pixelwise

以便循环发生在 c 中。您拥有每个像素长度未知的邻居列表这一事实使您的问题难以以这种方式并行化。您应该重新处理您的问题（您能更具体地了解您的算法吗？），或者使用另一种语言，例如 cython。

编辑：

如果不更改算法，您将无法从 Numpy 中受益。 Numpy 允许您执行线性代数运算。执行任意操作时，您无法避免使用此库循环开销。

为了优化这一点，您可以考虑：

• 切换到另一种语言，如 cython（专门用于 python 扩展）以消除循环成本

• 优化你的算法：如果你可以只使用线性代数运算（这取决于neighborsThatNeedCalculation 函数）得到相同的结果，你可以使用 numpy，但你需要制定一个新的架构。

• 使用 MapReduce 等并行化技术。使用 python，您可以使用一个工作池（在多处理模块中提供），如果您切换到另一种语言，您将获得更快的速度提升，因为 python 会有其他瓶颈。

如果您想要一些易于设置和集成的东西，并且您只需要具有类似 c 的性能，如果您不能重新设计您的算法，我强烈建议您使用 cython。

Answer 2

对于第一部分，您可以使用numpy.vectorize，但只有在无法使用数组操作来实现updateDistance 的功能时才应该这样做。这是一个例子：

import numpy as np    
updateDistance = np.vectorize(lambda x: x + 1) # my updateDistance increments

实际上，如果这是您尝试执行的操作，只需执行 a + 1。 因此，如果我们获取一个数组并应用 updateDistance：

>>> a = np.ones((3,3))
>>> updateDistance(a)
array([[ 2.,  2.,  2.],
       [ 2.,  2.,  2.],
       [ 2.,  2.,  2.]])

至于第二部分，我认为我对细节的理解不足以提出更好的选择。听起来您需要反复查看最近的邻居，所以我怀疑您至少可以改进if-else 中的内容。

---

更新：第一部分的时间安排。

注意：这些时间是在我的机器上完成的，没有尝试标准化环境。

循环时间通过以下方式生成：

python -mtimeit 'import numpy as np' 'n = 100' 'a = np.ones((n, n))' 'b = np.zeros((n, n))' 'for x in range(n): ' '    for y in range(n):' '        b[x,y] = a[x,y] + 1'

np.vectorize 时间是通过以下方式生成的：

python -mtimeit 'import numpy as np' 'n = 100' 'a = np.ones((n, n))' 'updateDistance = np.vectorize(lambda x: x + 1)' 'b = updateDistance(a)'

在这两种情况下，n = 100 都会导致一个 100 x 100 的数组。根据需要替换100。

Array size    Loop version    np.vectorize version    np.vectorize speed up
100 x 100     20.2 msec        2.6 msec               7.77x
200 x 200     81.8 msec       10.4 msec               7.87x
400 x 400     325 msec        42.6 msec               7.63x

最后，将np.vectorize 示例与简单地使用数组操作进行比较，您可以这样做：

python -mtimeit 'import numpy as np' 'n = 100' 'a = np.ones((n, n))' 'a += 1'

在我的机器上，这会生成以下内容。

Array size    Array operation version    Speed up over np.vectorize version
100 x 100     23.6 usec                  110.2x
200 x 200     79.7 usec                  130.5x
400 x 400     286 usec                   149.0x

总之，使用np.vectorize而不是循环有一个优势，但是如果可能的话，使用数组操作来实现updateDistance的功能会有更大的动力。

Answer 3

您应该考虑使用 C-extension/Cython。 如果您继续使用 Python，则可以通过替换来实现一项重大改进：

for xn,yn in neighbors(x,y):
      A[x,y] += 1+A[xn,yn]

与：

n = neighbors(x,y)
A[x,y] += len(n)+sum(A[n])

neighbors 应该返回索引，而不是下标。