关于Dijkstra算法和其他的一些图算法解决最短路径问题的方法的Python实现

1,333 阅读13分钟
原文链接: zhuanlan.zhihu.com

这篇文章其实主要想说的是如何解决最短路径的问题

其实最短路径问题,我们在生活中都在不知不觉的使用。比如我们在上网的时候,互联网传输采用了各种各样的数据包路由方法。这些路由算法都在幕后工作。

还有一些图算法的寻路操作,比如游戏中让游戏角色自动寻路。

其实寻找最短路径一般都是其他非图类算法中的一个重要的子程序。

下面开始说一下最短路径问题

最短路径问题可以分为以下几种形式。例如,在有向图和无向图中找出最短路径,二者最重要的区别在于起点和目的地。

这个问题有很多种形式,比如从一个节点到其他所有节点的最短路径的源节点唯一的问题,从一个节点到另一个节点的最短路径的一对一问题,还有从其他所有节点到一个节点的最短路径的目标节点唯一问题,还有从所有节点到其他节点的最短路径的所有节点两两组合问题。

进去正文之前,我们先说一下松弛化和逐步改进的概念,具体的例子不再举了,大家知道是一种求解方法就好,这些方法会通过逐步接近的方式来获得相关问题的最佳解法。

然后大家看一下这个代码:

# 松弛技术
inf = float('inf')


def relax(W, u, v, D, P):
    d = D.get(u, inf) + W[u][v]
    if d < D.get(v, inf):
        D[v], P[v] = d, u
        return True

把图表示为字典的字典,用字典D存放距离值估计。P作为前导节点字典。前导指针构成了最短路径树。然后用公共代码的relax分解出松弛技术。D中不存在的项目都可以看做无穷大。也在主算法对它们初始化定义为无穷大。

这段代码,我们通过u来观察是否可以缩短路径,从而改进目前已知的到达v点的最短路径。不是最短路径的话,不再理会。是最短路径的话,就记住经过了哪些节点。

然后我们看第一个比较正式的算法,Bellman-Ford算法。这个算法适用于任意有向或者无向图的单源最短路径算法。如果图包含负环,算法将会输出信息,然后放弃查找。

# Bellman-Ford算法
def bellman_ford(G, s):
    D, P = dict(s=0), dict()
    for rnd in G:
        changed = False
        for u in G:
            for v in G[u]:
                if relax(G, u, v, D, P):
                    changed = True
        if not changed:
            break
    else:
        raise ValueError('negative cycle')
    return D, P

这个算法包含了检查变化的charged,所以不需要多次迭代,程序可以提前终止。哈可以通过判断多余的迭代是否带来变化,来检测负环是否存在。

然后我们看一个权重图:

权重图改写为字典的字典,可以用如下代码表示:

a, b, c, d, e, f, g, h = range(8)
G = {
    a: {b: 2, c: 1, d: 3, e: 9, f: 4},
    b: {c: 4, e: 3},
    c: {d: 8},
    d: {e: 7},
    e: {f: 5},
    f: {c: 2, g: 2, h: 2},
    g: {f: 1, h: 6},
    h: {f: 9, g: 8}
}

算法具体的实现,可以用调试器,增加两条打印命令来显示松弛操作的边和分配给D的值。

如果测试几次,就会发现,g、h和f的距离值估计不断建议,到了第8轮还有减一的情况发生,这就提示了有负环的存在。

解决方法可以用del G[f][g]去掉(f, g)边。这样就不会出现负环的问题。

然后我们给单源DAG图最短路径问题加一个约束条件,环可以有,但边的权值不能为负值

然后我们引入这个Dijkstra算法,这个算法的结构和Prim算法很类似,都使用优先级队列进行遍历。在Dijkstra算法中,优先级是距离值估计。然后我们看具体的算法实现:

# Dijkstra算法
from heapq import heappush, heappop


def dijkstra(G, s):
    D, P, Q, S = {s: 0}, {}, [(0, s)], set()
    while Q:
        _, u = heappop(Q)
        if u in S:
            continue
        S.add(u)
        for v in G[u]:
            relax(G, u, v, D, P)
            heappush(Q, (D[v], v))
    return D, P

这个算法的运行时间是Q((m+n)lgn),其中m为边数,n为节点数。

和Dijkstra算法关系比较密切的还有广度优先搜索(BFS)算法。我们可以考虑一下边权值为正整数的情况,将权值为w的一条边替换为w-1条无权边,连成一条由虚拟节点组成的路径,比如下图:

BFS总会找到正确答案,但是效率会变得很低。这时候我们会发现BFS解决问题的方式和Dijkstra算法十分类似:在各条边上华为的时间与边权值成正比,所以会按照距离大小从起始节点依序到达各节点。

然后下面说多对多问题

这里引入Johnson算法,算法动机其实很简单,解决稀疏图矩阵所有节点对之间的最短路径问题,对各个节点相同的情况下使用Dijkstra算法。但是Dijkstra算法不允许负权边的存在。对于单源最短路径的问题,除了改用Bellman-Ford算法,并没有其他办法。

然后我们可以这样设想,增加一个新的节点s,它到所有现在的节点的边权值为0.然后对从s出发的情况运行Bellman-Ford算法。这样可以计算出s到图中每个节点的距离。我们称为h(v)。然后我们可以用h调节各边的权值。可以如下定义:w'(u, v) = w(u, v) + h(u) - h(v)。

这样定义可以保证新的权值非负值。然后也不会引入新的干扰项。

然后再用Dijkstra算法发现最短路径,再逆向变换所有的路径长度。具体的Johnson算法实现如下:

# 求解稀疏矩阵图
from copy import deepcopy


def johnson(G):
    G = deepcopy(G)
    s = object()
    G[s] = {v: 0 for v in G}
    h, _ = bellman_ford(G, s)
    del G[s]
    for u in G:
        for v in G[u]:
            G[u][v] += h[u] - h[v]
    D, P = dict(), dict()
    for u in G:
        D[u], P[u] = dijkstra(G, u)
        for v in G:
            D[u][v] += h[v] - h[u]
    return D, P

Johnson算法效率也还不错,不过它的的运行时间是Dijkstra算法运行时间的n倍。

然后开始说下一个问题,求解所有节点对之间最短距离的方法

这个算法名字叫做Floyd-Warshall算法。基于动态编程的原理,其实Dijkstra算法也是基于动态编程。

Floyd-Warshall基于缓存式递归分解,实现的过程一般都具有迭代性。

我们需要寻找一组递归相关的子问题。

我们随意对节点排序,并限制允许用于构成最短路径的中间节点的数量,即前k个。

这里直接用三个参数对子问题进行参数化:

  1. 起始节点
  2. 终止节点
  3. 允许经过的最大节点编号

然后设节点u到节点v的最短路径的长度为d(u, v, k):

d(u, v, k) = min(d(u, v, k-1), d(u, k, k-1) + d(k, v, k-1))。

这个和背包问题一样,要考虑是否包括节点k。不包括,就是d(u, v, k-1)。包括就必须使用到达k的最短路径d(u, k, k-1)和从k出来的最短路径d(k, v, k-1)。

然后下面看代码:

# Floyd-Warshall算法的缓存式递归实现
def rec_floyd_warshall(G):
    @memo
    def d(u, v, k):
        if k == 0:
            return G[u][k]
        return min(d(u, v, k-1), d(u, k, k-1) + d(k, v, k-1))
    return {(u, v): d(u, v, len(G)) for u in G for v in G}


# 记忆体化的装饰器的函数
from functools import wraps


def memo(func):
    cache = dict()

    @wraps(func)
    def wrap(*args):
        if args not in cache:
            cache[args] = func(*args)
        return cache[args]
    return wrap

然后我们可以尝试一下迭代的版本:

三个参数,就需要三个for循环。下面直接上代码:

# Floyd-Warshall算法,仅考虑距离
def floyd_warshall(G):
    D = deepcopy(G)
    for k in G:
        for u in G:
            for v in G:
                D[u][v] = min(D[u][v], D[u][k] + D[k][v])
    return D

D为当前的距离图,先前的距离图为C。全程只用一个距离图换算之后,公式由:

D[u][v] = min(D[u][v], C[u][k] + C[k][v])。

变换为:

D[u][v] = min(D[u][v], D[u][k] + D[k][v])

如果再加入一个P矩阵的话,P[u][v]将被替换为P[k][v]。代码就会变为:

# Floyd-Warshall算法
def floyd_warshall(G):
    D, P = deepcopy(G), dict()
    for u in G:
        for v in G:
            if u == v or G[u][v] == inf:
                P[u, v] = None
            else:
                P[u, v] = u
    for k in G:
        for u in G:
            for v in G:
                shortcut = D[u][k] + D[k][v]
                if shortcut < D[u][v]:
                    D[u][v] = shortcut
                    P[u, v] = P[k, v]
    return D, P

这个地方应为shortcut < D[u][v],而不是shortcut <= D[u][v],因为在某些情况,最后一步为D[v][v],这时候将导致P[u, v] = None。

然后说下一个问题,说一下“中途相遇的问题”

Dijkstra算法子问题的解,如果节点为s和t的话,将由s到t在图上不断扩散。但是如果从起点和终点同时出发,展开遍历,这样就会减少很多工作量。具体抽象出来就像下面的图:

把原来的稍作修改就可以变为,Dijkstra算法的双向图版本。可以变成一个子解生成器,让我们尽可能提取更多的子解。

这样必须抛开距离表,仅仅依靠优先队列中保存的距离值。于是放上示例代码:

# Dijkstra算法作为解决方案生成器的实现
from heapq import heappush, heappop


def idijkstra(G, s):
    Q, S = [(0, s)], set()
    while Q:
        d, u = heappop(Q)
        if u in S:
            continue
        S.add(u)
        yield u, d
        for v in G[u]:
            heappush(Q, (d+G[u][v], v))

没有引入前导节点信息,但是可以通过向堆元组中添加前导节点来扩展这个解决方案。获取距离表,只需调用dict(idijkstra(G, s))。

但是从两个节点同时出发,可能会遇到下面的问题:

首次相遇的节点(高亮的节点)并不一定位于最短路径(高亮的边)上。

其实这个问题只需要对终止条件做限定就好了,只要看他们已经走了多远,就是目前获得的最新的距离,这个是不能减小的。二者综合至少和我们目前发现的最佳路径相等,这样就不可能再找到更好的方案。

如果G是无向图的话,并且任意节点u满足G[u][u] = 0,则可以使用下面的代码:

# Dijkstra算法的双向图版本
from itertools import cycle


def bidir_dijkstra(G, s, t):
    Ds, Dt = dict(), dict()
    forw, back = idijkstra(G, s), idijkstra(G, t)
    dirs = (Ds, Dt, forw), (Dt, Ds, back)
    try:
        for D, other, step in cycle(dirs):
            v, d = next(step)
            D[v] = d
            if v in other:
                break
    except StopAsyncIteration:
        return inf
    m = inf
    for u in Ds:
        for v in G[u]:
            if not v in Dt:
                continue
            m = min(m, Ds[u] + G[u][v] + Dt[v])
    return m

当然,你也可以很轻松的拓展这段代码。

下面我们会引入A*算法

如果真的知道哪一条边可以离目标更近,我们就可以用贪婪算法解决问题。直接沿着最短的路径移动,不用理会其他的旁支路线。

A*算法有点像人工智能中启发式搜索的概念。而不是像DFS和BFS那样子盲目搜索。也不像Dijkstra算法一样对未来走向一无所知。

因为A*加入了一个潜在势函数,也可以叫做启发式函数h(v)。

就像上面介绍的Johnson算法一样,我们可以定义修正后的边权:

w'(u, v) = w(u, v) - h(u) + h(v)

然后你会发现,这样子调整之后,我们可以奖励正确的边,惩罚不正确的边。然后给各个边的权值加上了启发式番薯。这个算法将沿导致剩余距离下降的方向发展。

A*算法相当于针对修正图的Dijkstra算法,h可行的话,算法就是正确的。

h(s)是一个常数,所以我们只用将h(v)加入现有的优先级的队列中。这个和是我们对从s到t的最佳估计。如果w'(u, v)可行,那么h(v)也会是d(v, t)的下界。

然后这只直接给出代码,里面调用了上面代码的idijkstra函数:

# A*算法
from heapq import heappush, heappop

inf = float('inf')


def a_star(G, s, t, h):
    P, Q = dict(), [(h(s), None, s)]
    while Q:
        d, p, u = heappop(Q)
        if u in P:
            continue
        P[u] = p
        if u == t:
            return d - h(t), P
        for v in G[u]:
            w = G[u][v] - h(u) + h(v)
            heappush(Q, (d + w, u, v))
    return inf, None

A*算法的优势就在于启发式函数,至于具体的启发式函数是什么取决于要解决的问题。

另外,A*算法也可以搜索解空间。可以解决魔方问题或者单词梯问题。这里说一下后一个问题。

单词梯是从第一个起始单词开始构建的,比如 lead,然后以另外一个单词作为结尾,比如gold。单词梯每一步搭建都要用到实际单词,从一个单词推进到下一个单词,只能更换一个字母。比如这个题目的一个解法就可以通过 load 和 goad 这两个单词,到达 lead 和 gold。如果每个单词看做一个节点,我们可以加上边。可能没有必要这样建立,不过这样假设,便可以做出如下代码:

# 单词梯路径的隐式图
from string import ascii_letters as chars


def variants(wd, words):
    wasl = list(wd)
    for i, c in enumerate(wasl):
        for oc in chars:
            if c == oc:
                continue
            wasl[i] = oc
            ow = ''.join(wasl)
            if ow in words:
                yield  ow
        words[i] = c

class WordSpace:

    def __init__(self, words):
        self.words = words
        self.M = dict()

    def __getitem__(self, wd):
        if wd not in self.M:
            self.M[wd] = dict.fromkeys(self.variants(wd, self.words), 1)
        return self.M[wd]

    def heuristic(self, u, v):
        return sum(a!=b for a, b in zip(u, v))

    def ladder(self, s, t, h=None):
        if h is None:
            def h(v):
                return self.heuristic(v, t)
        _, p = a_star(self, s, t, h)
        if p is None:
            return [s, None, t]
        u, p = t, []
        while u is not None:
            p.append(u)
            u = P[u]
        p.reverse()
        return p

WordSpace可以作为加权图使用,配合a_star使用。这里的启发式函数只是统计单词不同的字符位数。G是WordSpace的对象,G['lead'] 就是字典,其他单词是键值。各个边的权值是1。

文章写到了这里,差不多就说完了。

A*算法一般在人工智能的书里面才会涉及到。其他的一些算法在一般的算法书中也能找得到。

设计新的算法的时候,比如Dijkstra算法的双向版本和A*的启发式结合在一起,很可能因为一些陷阱而导致算法无效。

以上。谢谢大家关注。

今天立冬,祝大家冬天愉快,所有的不开心伴随着秋天的过去都过去了。新的季节来了,希望大家都有好心情。

天气寒冷,大家注意保暖。