是否有可能优化这个动态规划代码

Question

对于200,000个浮点数的数据集，这段代码花费了半个多小时。

import numpy as np
try:
    import progressbar
    pbar = progressbar.ProgressBar(widgets=[progressbar.Percentage(),
        progressbar.Counter('%5d'), progressbar.Bar(), progressbar.ETA()]) 
except:
    pbar = list
block_length = np.loadtxt('bb.txt.gz') # get data file from http://filebin.ca/29LbYfKnsKqJ/bb.txt.gz (2MB, 200000 float numbers)
N = len(block_length) - 1
# arrays to store the best configuration
best = np.zeros(N, dtype=float)
last = np.zeros(N, dtype=int)
log = np.log
# Start with first data cell; add one cell at each iteration
for R in pbar(range(N)):
    # Compute fit_vec : fitness of putative last block (end at R)
    #fit_vec = fitfunc.fitness(
    T_k = block_length[:R + 1] - block_length[R + 1]
    #N_k = np.cumsum(x[:R + 1][::-1])[::-1]
    N_k = np.arange(R + 1, 0, -1)
    fit_vec = N_k * (log(N_k) - log(T_k))
    prior = 4 - log(73.53 * 0.05 * ((R+1) ** -0.478))
    A_R = fit_vec - prior #fitfunc.prior(R + 1, N)
    A_R[1:] += best[:R]
    i_max = np.argmax(A_R)
    last[R] = i_max
    best[R] = A_R[i_max]
# Now find changepoints by iteratively peeling off the last block
change_points = np.zeros(N, dtype=int)
i_cp = N
ind = N
while True:
    i_cp -= 1
    change_points[i_cp] = ind
    if ind == 0:
        break
    ind = last[ind - 1]
    change_points = change_points[i_cp:]
print edges[change_points] # show result

第一个循环非常慢，因为数组的长度在每次迭代时都是R，即增加，导致复杂度为N^2。

是否有办法进一步优化这段代码，例如通过预计算?我也很高兴使用其他编程语言的解决方案。

Answer 1

我可以复制A_R(直到fit-prior步骤)作为上三角形NxN矩阵，具有:

def trilog(n):
    nn = n[:-1,None]-n[None,1:]
    nn[np.tril_indices_from(nn,-1)]=1
    return nn
T_k = trilog(block_length)
N_k = trilog(-np.arange(N+1))
fit_vec = N_k * (np.log(N_k) - np.log(T_k))
R = np.arange(N)+1
prior = 4 - log(73.53 * 0.05 * (R ** -0.478))
A_R = fit_vec - prior
A_R = np.triu(A_R,0)
print(A_R)

我还没有通过计算逻辑和应用best。

我只对小数组这样做过。对于你的完整问题，对应的矩阵对于我的内存来说太大了。

B=np.ones((200000,200000),float)

因此，仅从内存考虑，您可能会坚持使用for R in range(N)迭代。