全源最短路径
全源最短路径问题是指给定带权有向图,求任意两点\(u,v\)从\(u\)到\(v\)的最短路径”。可把该问题直接作为\(|V|\)个单源最短路径问题,使用Dijkstra算法(二叉堆)的代价为\(O(|V||E|lg|V|)\),Bellman-Ford算法的代价为\(O(|V|^2|E|)\),对于稠密图二者的代价可分别写成\(O(|V|^{3}lg|V|)\)与\(O(|V|^4)\)。本文将讨论专门用于解决全源最短路问题的算法,这些算法大多基于邻接矩阵讨论,且不允许存在任何负权环路的。下面给出相关定义:
1) 带权邻接矩阵:带权有向图\(G=(V,E)\)与其权重可被\(|V|\times |V|\)的矩阵\(W\)描述。其中\(W_{ij}\)表示图边\((i,j)\)的权重,并且当\(i=j\)时有\(W_{ij}=0\),当图边\((i,j)\)不存在时有\(W_{ij}=\infty\);
2) 最短路径值矩阵:用矩阵\(D\)描述最短路径值,其中\(D_{ij} = \delta(i,j)\);
3) 前驱节点矩阵:当\(i=j\)或\(i\)不可达\(j\)时,\(\Pi_{ij}=None\)。否则\(\Pi_{ij}\)为\(i\)到\(j\)某个最短路中\(j\)的前驱;
基于矩阵乘法的算法
该算法的限制条件和用到的性质为:
a) 最优子结构:若\(i-…-k-j\)为最短路,则\(i-…-k\)为最短路,则\(D_{ij}=D_{ik}+W_{kj}\);
b) 有环图:该算法不允许负权环路,故若\(i\)可达\(j\)则必然存在无环最短路,边数至多为\(|V|-1\);
然后通过矩阵\(L^{(m)}\)给出递归解的结构:
1) 设矩阵\(L_{ij}^{(m)}\)为从\(i\)到\(j\)且至多含\(m\)条边的最短路值,则\(L^{(1)}=W\);
2) 规定\(m=0\)时,若\(i=j\),则\(L_{ij}^{(0)}=0\),否则\(L_{ij}^{(0)}=\infty\);
3) \(L_{ij}^{(m)} = min\{ L_{ij}^{(m-1)},~ min\{ ~ L_{ik}^{(m-1)} + W_{kj} ~|~ k\in V ~\} ~ \} = min\{ ~ L_{ik}^{(m-1)} + W_{kj} ~|~ k\in V ~\}\);
4) 根据上述(b)可知\(D_{ij} = L_{ij}^{(|V|-1)} = L_{ij}^{(|V|)} = L_{ij}^{(|V|+1)} = … = L_{ij}^{(\infty)}\);
类比\(A = L^{(m-1)} \times W\)的矩阵乘法,分析利用\(L^{(m-1)}\)与\(W\)求\(L^{(m)}\)的流程,发现把(i)的加法换为\(min(x,y)\),乘法换为加法即得(ii),该运算(计为\(\times\))是结合的,故\(L^{(m)}= L^{(m-1)} \times W = W^m\)。最短路值矩阵\(D\)可归纳为(iii)。
i) \(A_{ij} = \sum_{k\in V } (L_{ik}^{(m-1)} \cdot W_{kj}) = add(add(add( L_{i1}^{(m-1)} \cdot W_{1j} , L_{i2}^{(m-1)} \cdot W_{2j} ), L_{i3}^{(m-1)} \cdot W_{3j}),… \);
ii) \(L_{ij}^{(m)} = min\{ ~ L_{ik}^{(m-1)} + W_{kj} ~|~ k\in V \} = min(min( L_{i1}^{(m-1)} + W_{1j} , L_{i2}^{(m-1)} + W_{2j} ),… \);
iii) 定义矩阵乘法\((AB)_{ij} = min\{ A_{ik} + B_{kj} |k\in V \}\),给定\(n\times n\)的带权邻接矩阵\(W\),则\(D=W^{n-1}\);
考虑基于朴素的矩阵乘法与“重复平方”策略将上述分析转化为算法实现:
1) 对于两个\(n\times n\)方阵间的乘法,朴素矩阵乘法算法需要\(O(n^3)\)的代价,故对于\(n\times n\)方阵,计算\(A^{n-1}\)需要\(O(n^4)\)代价,故类比朴素算法计算上述“\(D=W^{n-1}\)”同样为\(O(n^4)=O(|V|^4)\)代价;
2) 可通过二分思路对矩阵幂次做“重复平方”优化,其流程为迭代的计算\((((W^2)^2)^2)^2…\),第\(i\)次迭代求得\(W^{2^i}\),每次迭代完若\(2^i \geq |V|-1 \)则该次迭代值即为解,无需严格要求\(2^i = |V|-1 \)是因为上述(3)给出的性质;
3) 由于每次迭代可看作求上次迭代结果的平方,故\(i\)次迭代需\(i\)次矩阵乘法,总乘法次数\(k\)有上界\(2^{k+1} \leq |V|-1 \),即\(k = O(lg|V|)\),“重复平方”算法的最坏代价为\(O(|V|^3lg|V|)\);
最后讨论前驱节点矩阵,设\(\Pi_{ij}^{(m)}\)为从\(i\)到\(j\)且至多含\(m\)条边的某最短路中\(j\)的前驱,通过从左到右逐步“相乘”的方式求\(L_{ij}^{(m)}\),即朴素乘法\((((W)W)W)…\)的计算方法,随这个流程记录取\(min\)时的具体顶点可求得\(\Pi_{ij}^{(m)}\)并整理得出\(\Pi_{ij}\)。但若用“重复平方”策略,计算\((L^{(k)})^2 \)对应\(\Pi^{(2k)}\)无明确意义,故上述算法未构造前驱节点矩阵\(\Pi\)。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | ''' W = [ 0,2,3, 4,0,float("inf"), 7,8,0, ] ''' def extendMultiply(A,B): n = int(len(A) ** 0.5) C = [] for i in range(n): for j in range(n): C.append(min( [ A[k+i*n] + B[j+k*n] for k in range(n) ] )) return C def fasterShorestPath(W): n = int(len(W) ** 0.5) m = 1 while m < n - 1: W = extendMultiply(W,W) m = 2 * m return W |
Floyd-Warshall算法
该算法的限制条件和用到的性质为:
a) 顶点的表示:用邻接矩阵的行标(列标)表示顶点,比如对于\(n\times n\)邻接矩阵,其顶点集为\(\{1,2…,n\}\);
b) 中间节点:简单路径\(p=1->2->…(n-1)->n\)的中间节点为集合\(\{2,…,{n-1}\}\);
c) 最优子结构:设\(V=\{1,2,…,n\}\)。其有子集\(\{1,2,…,k\}\),其中\(k<n\)。对于任意节点对\((i,j)\),考虑其所有中间节点取自\(\{1,2,…,k\}\)的路径,\(p\)为其中最短路。若\(k\)不是\(p\)的中间节点,则\(p\)的中间节点都取自\(\{1,2,…,k-1\}\),故中间节点取自\(\{1,2,…,k-1\}\)的最短路也是中间节点取自\(\{1,2,…,k\}\)的最短路。若\(k\)是\(p\)的中间节点,则可将\(p\)分解为两个子最短路\(p_1 = i -…- k\)与\(p_2 = k-…-j\)。那么\(p_1\)与\(p_2\)的中间节点都取自\(\{1,2,…,k-1\}\);
d) 有环图:该算法同样不允许负权环路;
然后通过矩阵\(D^{(m)}\)给出最短路值矩阵\(D\)的递归结构:
1) 设\(D_{ij}^{(k)}\)为从\(i\)到\(j\)所有中间节点都取自\(\{1,2…k\}\)的最短路的权重;
2) \(D_{ij}^{(0)}\)表示无中间节点的情况,所以规定\(D_{ij}^{(0)}=W_{ij}\);
3) \(k\)的其他情况根据最优子结构,\(D_{ij}^{(k)} = min\{~ D_{ij}^{(k-1)} , D_{ik}^{(k-1)} + D_{kj}^{(k-1)} ~\}\);
4) 设\(n\)为图顶点数,则最终\(D_{ij} = \delta(i,j) = D_{ij}^{(n)}\);
然后通过矩阵\(\Pi^{(m)}\)给出前驱节点矩阵\(\Pi\)的递归结构:
1) 设\(\Pi^{(k)}\)为从\(i\)到\(j\)的某条所有中间节点都来自\(\{1,2…k\}\)的最短路中\(j\)的前驱;
2) \(\Pi_{ij}^{(0)}\)表示无中间节点的情况,当\(i=j\)或\(W_{ij}=\infty\)时\(\Pi_{ij}^{(0)}=None\),否则\(W_{ij}=j\);
3) \(k\)的其他情况,若\(D_{ij}^{(k-1)} \leq D_{ik}^{(k-1)} + D_{kj}^{(k-1)} \)有\(\Pi_{ij}^{(k)} = \Pi_{ij}^{(k-1)} \);
4) \(k\)的其他情况,若\(D_{ij}^{(k-1)} > D_{ik}^{(k-1)} + D_{kj}^{(k-1)} \)有\(\Pi_{ij}^{(k)} = \Pi_{kj}^{(k-1)} \);
最后通过自底向上的迭代次序实现该算法,可以发现该算法属于DP算法,其代价为\(O(|V|^3)\);
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 | ''' W = [ 0,2,3, 4,0,float("inf"), 7,8,0, ] ''' def initMatrix(W): n = int(len(W) ** 0.5) PI = [] for i in range(n): for j in range(n): if i==j or W[j+n*i]==float('inf'): PI.append(None) else: PI.append(j) kPI = [PI[:] for i in range(n+1)] kD = [W[:] for i in range(n+1)] return kD,kPI def floyd_warshall(W): n = int( len(W) ** 0.5 ) D,P = initMatrix(W) for k in range(1,n+1): for i in range(n): for j in range(n): x = D[k-1][j+i*n] y = D[k-1][(k-1)+i*n] + D[k-1][j+(k-1)*n] if x <= y: D[k][j+i*n] = x P[k][j+i*n] = P[k-1][j+i*n] else: D[k][j+i*n] = y P[k][j+i*n] = P[k-1][j+(k-1)*n] return D[-1],P[-1] print(floyd_warshall([ 0,1,4, 1,0,2, 1,5,0, ])) |
Johnson算法
该算法的限制条件和用到的性质为:
a) 有环图:该算法同样不允许负权环路,但其可以检测到负权环路;
b) “重新赋予权重”:若图中无负权边,则对每个节点执行Dijkstra算法即可。若图中存在负权边,但无负权环路,则需要为每个边计算一组新的非负权重值,然后执行Dijkstra算法即可。其中新权重函数\(w’\)必须满足如下性质:
b.i) 原权重\(w\)下的任意节点对间的任意最短路\(p\),新权重\(w’\)必须使\(p\)仍为最短路;
b.ii) 任意边的新权重\(w’\)是非负值;
c) “重新赋予权重”引理:\(w'(u,v)=w(u,v)+h(u)-h(v)\),\(h(x)\)为任意函数。设\(p\)为\(u\)到\(v\)的任意路径。
c.i) \(p\)为\(w\)下的最短路径 \(<=> \) \(p\)为\(w’\)下的最短路径;
证明:对\(p\)求权重和,错位相减得\(\sum_{x \in p}w'(x)=(\sum_{x \in p}w(x))+h(u)-h(v)=(\sum_{x \in p}w(x))+C\)。故从\(u\)到\(v\)的任意路径的总权重,其与原权重的差是固定的常数,原权重最短路必为新权重最短路。
c.ii) \(G\)在\(w\)下无负权环路 \(<=> \) \(G\)在\(w’\)下无负权环路;
证明:若\(p\)为环路,\(\sum_{x \in p}w'(x)=(\sum_{x \in p}w(x))+h(u)-h(u)=\sum_{x \in p}w(x)\)。
给出算法的流程如下,其中步骤(1)(2)的方法在之前差分约束的文章中讨论过。
1) 在\(G\)中加入新顶点\(s\),并初始化\(s\)有到达所有其他顶点的0权重出边,把更改后的图计为\(G’\);
2) 在\(G’\)的\(s\)运行Bellman-Ford算法,求得每个\(\delta(s,v)\)并判断原图\(G\)有无负权环路(差分约束文章里讨论过);
3) 设\(h(v)=\delta(s,v)< \infty\),利用三角不等式有\(w'(u,v)=w(u,v)+h(u)-h(v) \geq 0\);
4) 然后以\(w’\)为权重,对\(G\)的每个顶点执行Dijkstra算法;
5) 原权重与新权重的最短路值满足\(\delta(s,v) + h(s) – h(v) = \delta'(s,v)\)(利用和证明引理相同的方法可证);
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 | class Graph: def __init__(self): self.d = dict() def addVertex(self,i): self.d[i] = dict() def addEdge(self,i1,i2,w): self.d[i1][i2] = w def getWeight(self,i1,i2): return self.d[i1][i2] def setWeight(self,i1,i2,w): self.addEdge(i1,i2,w) def getVertices(self,i): return self.d[i].keys() def vertices(self): return self.d.keys() class MinHeap: @staticmethod def parent(i):return 0 if i==0 else (i//2-1 if i%2==0 else i//2) @staticmethod def left(idx):return 2 * idx + 1 @staticmethod def right(idx):return 2 * idx + 2 def _heapify(self, idx, lastIdx): arr = self.arr keys = self.keys indexs = self.indexs l = self.left(idx) r = self.right(idx) x = idx if l <= lastIdx and keys[arr[l]] < keys[arr[idx]]: x = l if r <= lastIdx and keys[arr[r]] < keys[arr[x]]: x = r if x != idx: arr[x], arr[idx] = arr[idx], arr[x] indexs[arr[x]], indexs[arr[idx]] = indexs[arr[idx]], indexs[arr[x]] self._heapify(x,lastIdx) # 数组lst存放元组( 优先级值(key), 对象 ) def __init__(self,lst): self.arr = [] self.keys = {} self.indexs = {} for tp in lst: self.arr.append(tp[1]) self.indexs[tp[1]] = len(self.arr)-1 self.keys[tp[1]] = tp[0] def contain(self,k): return self.keys.get(k)!=None def pop(self): arr = self.arr keys = self.keys indexs = self.indexs if len(arr) == 0: return None if len(arr) == 1: k = arr.pop() keys.pop(k) indexs.pop(k) return k last = len(arr) - 1 arr[last], arr[0] = arr[0], arr[last] indexs[arr[last]], indexs[arr[0]] = indexs[arr[0]], indexs[arr[last]] self._heapify(0,len(arr)-2) k = arr.pop() keys.pop(k) indexs.pop(k) return k # 把堆中的对象k的优先级降为key并维持堆性质 def decrease(self,k,key): arr = self.arr keys = self.keys indexs = self.indexs i = self.indexs[k] keys[arr[i]] = key p = self.parent(i) while i > 0 and keys[arr[p]] > keys[arr[i]]: arr[i], arr[p] = arr[p], arr[i] indexs[arr[i]], indexs[arr[p]] = indexs[arr[p]], indexs[arr[i]] i = p p = self.parent(i) def dijkstra_relax(d,pi,graph,heap,u,v): l = d[u] + graph.getWeight(u,v) if l < d[v]: if heap.contain(v): heap.decrease(v,l) d[v] = l pi[v] = u def dijkstra(graph,s): d,pi = {},{} for v in graph.vertices(): if v!=s: d[v] = float('inf') pi[v] = None heap = [(v,k) for k,v in d.items()] d[s] = 0 heap.insert(0,(0,s)) heap = MinHeap(heap) while len(heap.arr)>0: u = heap.pop() for v in graph.getVertices(u): dijkstra_relax(d,pi,graph,heap,u,v) return d,pi def bellman_ford_relax(d,pi,graph,u,v): l = d[u] + graph.getWeight(u,v) if l < d[v]: d[v] = l pi[v] = u def bellman_ford(graph,s): d = {} pi = {} for v in graph.vertices(): d[v] = float('inf') pi[v] = None d[s] = 0 times = len(d) while times > 0: times -= 1 for u in d.keys(): for v in graph.getVertices(u): bellman_ford_relax(d,pi,graph,u,v) for u in d.keys(): for v in graph.getVertices(u): if d[u] + graph.getWeight(u,v) < d[v]: return False return d,pi def copy(graph): g = Graph() for u in graph.vertices(): g.addVertex(u) for u in graph.vertices(): for v in graph.getVertices(u): w = graph.getWeight(u,v) g.addEdge(u,v,w) return g def johnson(graph): g = copy(graph) s = 'THIS IS A STRING' g.addVertex(s) for k in g.vertices(): if k != s: g.addEdge(s,k,0) h = bellman_ford(g,s) if not h: return False h = h[0] g = copy(graph) for u in g.vertices(): for v in g.getVertices(u): w = g.getWeight(u,v) g.setWeight(u,v,w + h[u] - h[v]) results = [] for u in g.vertices(): d,pi = dijkstra(g,u) for v in d.keys(): d[v] = d[v] - h[u] + h[v] results.append((u,d,pi)) return results |