用python/numpy加速动态编程

Question

我有一个2D成本矩阵M，可能是400x400，我正试图计算通过它的最佳路径。因此，我的函数如下：

M[i,j] = M[i,j] + min(M[i-1,j-1],M[i-1,j]+P1,M[i,j-1]+P1)

显然是递归的。P1是加性常数。我的代码或多或少是这样的：

^{pr2}$

现在我知道在Numpy中循环是个糟糕的主意，对于像计算初始成本矩阵这样的事情，我已经找到了缩短时间的捷径。然而，由于我需要对整个矩阵进行潜在的评估，所以我不确定其他方法。在我的机器上，每次调用大约需要3秒，并且必须应用于大约300个这样的成本矩阵。我不知道这个时间从何而来，因为profiling说200000次对min的调用只需要0.1s——也许是内存访问？在

有没有一种方法可以同时做到这一点呢？我想可能有，但对我来说，似乎每次迭代都是相互依赖的，除非有更聪明的方法来记忆东西。在

这个问题有相似之处：Can I avoid Python loop overhead on dynamic programming with numpy?

如果有必要，我很乐意切换到C语言，但是我喜欢Python在快速测试方面的灵活性，以及缺少文件IO的faff。在我的脑海中，像下面这样的代码是否可能会快得多？在

#define P1 10
void optimalcost(double** costin, double** costout){
    /* 
        We assume that costout is initially
        filled with costin's values.
    */
    float a,b,c,prevcost;

    for(i=0;i<400;i++){
        for(j=0;j<400;j++){
            a = prevcost+P1;
            b = costout[i][j-1]+P1;
            c = costout[i-1][j-1];
            costout[i][j] += min(prevcost,min(b,c));
            prevcost = costout[i][j];
        }
    }
}

return;

更新：

我在Mac上，我不想安装一个全新的Python工具链，所以我使用了Homebrew。在

> brew install llvm --rtti
> LLVM_CONFIG_PATH=/usr/local/opt/llvm/bin/llvm-config pip install llvmpy
> pip install numba

新的“编号”代码：

from numba import autojit, jit
import time
import numpy as np

@autojit
def cost(left, right):
    height,width = left.shape
    cost = np.zeros((height,width,width))

    for row in range(height):
        for x in range(width):
            for y in range(width):
                cost[row,x,y] = abs(left[row,x]-right[row,y])

    return cost

@autojit
def optimalcosts(initcost):
    costs = zeros_like(initcost)
    for row in range(height):
        costs[row,:,:] = optimalcost(initcost[row])
    return costs

@autojit
def optimalcost(cost):
    width1,width2 = cost.shape
    P1=10
    prevcost = 0.0
    M = np.array(cost)
    for i in range(1,width1):
        for j in range(1,width2):
            M[i,j] += min(M[i-1,j-1],prevcost+P1,M[i,j-1]+P1)
            prevcost = M[i,j]
    return M

prob_size = 400
left = np.random.rand(prob_size,prob_size)
right = np.random.rand(prob_size,prob_size)

print '---------- Numba Time ----------'
t =  time.time()
c = cost(left,right)
optimalcost(c[100])
print time.time()-t

print '---------- Native python Time --'
t =  time.time()
c = cost.py_func(left,right)
optimalcost.py_func(c[100])
print time.time()-t

用Python编写代码是很有趣的，它是如此的非Python。注意对于任何对编写Numba代码感兴趣的人，您需要在代码中显式地表示循环。以前，我有一个整洁的纽比一个班轮

abs(left[row,:][:,newaxis] - right[row,:])

计算成本。那花了他7秒钟的时间。正确地写出循环可以得到0.5s

将其与本机Python代码进行比较是不公平的，因为Numpy可以很快做到这一点，但是：

Numba编译：0.509318113327s

母语：172.70626092s

我对这些数字和转换是多么的简单给我留下了深刻的印象。在

Answer 1

如果您不难切换到Python的Anaconda发行版，那么您可以尝试使用Numba，对于这种特殊的简单动态算法，这可能会在不让您离开Python的情况下提供大量的加速。在

Answer 2

Numpy通常不擅长迭代作业（尽管它确实有一些常用的迭代函数，如np.cumsum、np.cumprod、np.linalg.*等）。但对于上面的简单任务，如寻找最短路径（或最低能量路径），您可以通过思考同时可以计算的内容来将问题矢量化（同时尽量避免复制：

假设我们正在“行”方向（即水平方向）寻找最短路径，我们可以首先创建算法输入：

# The problem, 300 400*400 matrices
# Create infinitely high boundary so that we dont need to handle indexing "-1"
a = np.random.rand(300, 400, 402).astype('f')
a[:,:,::a.shape[2]-1] = np.inf

然后准备一些我们稍后使用的实用程序数组（创建需要固定的时间）：

^{pr2}$

最后进行计算并计时：

%%timeit
for i in np.arange(a.shape[1]-1):
    A[i].min(2, T)
    B[i] += T

我的（超旧笔记本电脑）的计时结果是1.78秒，这已经快过3分钟了。我相信你可以通过优化内存布局和对齐（以某种方式）来改进更多（同时坚持numpy）。或者，您可以简单地使用multiprocessing.Pool。它很容易使用，而且这个问题很容易分解成更小的问题（通过在批处理轴上划分）。在