如何计算Dijkstra算法中的跳跃次数？

Question

在numpy中，计算dijkstra算法使用的跳跃次数的最快方法是什么？我有一个10000x10000的连接矩阵，使用scipy.sparse.csgraph.dijkstra来计算填充的距离矩阵和前驱矩阵。我的简单做法如下：

import numpy as np
from numpy.random import rand
from scipy.sparse.csgraph import dijkstra

def dijkway(dijkpredmat, i, j):
    """calculate the path between two nodes in a dijkstra matrix"""
    wayarr = []
    while (i != j) & (j >= 0):
        wayarr.append(j)
        j = dijkpredmat[i,j]
    return np.array(wayarr)

def jumpvec(pmat,node):
    """calculate number of jumps from one node to all others"""
    jumps = np.zeros(len(pmat))
    jumps[node] = -999
    while 1:
        try:
            rvec = np.nonzero(jumps==0)[0]
            r = rvec.min()
            dway = dijkway(pmat, node, r)
            jumps[dway] = np.arange(len(dway),0,-1)
        except ValueError:
            break
    return jumps

#Create a matrix
mat = (rand(500,500)*20)
mat[(rand(50000)*500).astype(int), (rand(50000)*500).astype(int)] = np.nan
dmat,pmat = dijkstra(mat,return_predecessors=True)

timeit jumpvec(pmat,300)
In [25]: 10 loops, best of 3: 51.5 ms per loop

每个节点大约需要50毫秒是可以接受的，但当距离矩阵扩展到10000个节点时，时间增加到每个节点大约需要2秒。jumpvec也需要执行10000次……

Answer 1

这里有一个（渐进最优的）O(n)算法。

首先，创建一个未访问的顶点集合，最开始这个集合里包含除了源顶点以外的所有顶点。然后，把源顶点到它自己的跳跃向量的值初始化为0。接下来，只要这个集合不为空，就取出一个元素v。利用前驱矩阵，收集v以及每一个后续的祖先，直到找到一个已经访问过的顶点为止。然后，反向遍历这个列表，把每个节点w的跳跃向量值设置为它父节点的值加1，然后把w从集合中移除。

Answer 2

以下代码在我的电脑上可以加速4倍，之所以更快，是因为：

• 使用 ndarray.item() 从数组中获取值。
• 使用集合对象来保存未处理的索引。
• 在循环中不创建 numpy.arange()。

Python代码：

def dijkway2(dijkpredmat, i, j):
    wayarr = []
    while (i != j) & (j >= 0):
        wayarr.append(j)
        j = dijkpredmat.item(i,j)
    return wayarr

def jumpvec2(pmat,node):
    jumps = np.zeros(len(pmat))
    jumps[node] = -999
    todo = set()
    for i in range(len(pmat)):
        if i != node:
            todo.add(i)

    indexs = np.arange(len(pmat), 0, -1)
    while todo:
        r = todo.pop()
        dway = dijkway2(pmat, node, r)
        jumps[dway] = indexs[-len(dway):]
        todo -= set(dway)
    return jumps

为了进一步加速，你可以使用cython：

import numpy as np
cimport numpy as np
import cython

@cython.wraparound(False)
@cython.boundscheck(False)
cpdef dijkway3(int[:, ::1] m, int i, int j):
    cdef list wayarr = []
    while (i != j) & (j >= 0):
        wayarr.append(j)
        j = m[i,j]
    return wayarr

@cython.wraparound(False)
@cython.boundscheck(False)
def jumpvec3(int[:, ::1] pmat, int node):
    cdef np.ndarray jumps
    cdef int[::1] jumps_buf
    cdef int i, j, r, n
    cdef list dway
    jumps = np.zeros(len(pmat), int)
    jumps_buf = jumps
    jumps[node] = -999

    for i in range(len(jumps)):
        if jumps_buf[i] != 0:
            continue
        r = i
        dway = dijkway3(pmat, node, r)
        n = len(dway)
        for j in range(n):
            jumps_buf[<int>dway[j]] = n - j
    return jumps

这是我的测试，使用cython的版本快了80倍：

%timeit jumpvec3(pmat,1)
%timeit jumpvec2(pmat, 1)
%timeit jumpvec(pmat, 1)

输出结果：

1000 loops, best of 3: 138 µs per loop
100 loops, best of 3: 2.81 ms per loop
100 loops, best of 3: 10.8 ms per loop