1. 程式人生 > >帶權最短路 Dijkstra、SPFA、Bellman-Ford、ASP、Floyd-Warshall 演算法分析

帶權最短路 Dijkstra、SPFA、Bellman-Ford、ASP、Floyd-Warshall 演算法分析

圖論中,用來求最短路的方法有很多,適用範圍和時間複雜度也各不相同。

本文主要介紹的演算法的程式碼主要來源如下:

  1. Dijkstra: Algorithms(《演算法概論》)Sanjoy Dasgupta, Christos Papadimitriou, Umesh Vazirani;《演算法競賽入門經典—訓練指南》劉汝佳、陳峰。
  2. SPFA (Shortest Path Faster Algorithm): Data Structure and Algorithms Analysis in C, 2nd ed.(《資料結構與演算法分析》)Mark Allen Weiss.
  3. Bellman-Ford: Algorithms
    (《演算法概論》)Sanjoy Dasgupta, Christos Papadimitriou, Umesh Vazirani.
  4. ASP (Acyclic Shortest Paths): Introduction to Algorithms - A Creative Approach(《演算法引論—一種創造性方法》)Udi Manber.
  5. Floyd-Warshall: Data Structure and Algorithms Analysis in C, 2nd ed.(《資料結構與演算法分析》)Mark Allen Weiss.

它們的使用限制和執行時間如下:

Dijkstra: 不含負權。執行時間依賴於優先佇列的實現,如 O((∣V∣+∣E∣)log∣V∣)O((∣V∣+∣E∣)log∣V∣)


SPFA: 無限制。執行時間O(k⋅∣E∣) (k≪∣V∣)O(k⋅∣E∣) (k≪∣V∣)
Bellman-Ford:無限制。執行時間O(∣V∣⋅∣E∣)O(∣V∣⋅∣E∣)
ASP: 無圈。執行時間O(∣V∣+∣E∣)O(∣V∣+∣E∣)
Floyd-Warshall: 無限制。執行時間O(∣V∣3)

其中 1~4 均為單源最短路徑 (Single Source Shortest Paths) 演算法; 5 為全源最短路徑 (All Pairs Shortest Paths) 演算法。順便說一句,為什麼沒有點對點的最短路徑?如果我們只需要一個起點和一個終點,不是比計算一個起點任意終點更節省時間麼?答案還真不是,目前還沒有發現比算從源點到所有點更快的演算法。

圖的表示

本文中,前四個演算法的圖都採用鄰接表表示法,如下:

[cpp] view plain copy print?
  1. struct Edge  
  2. {  
  3.     int from;  
  4.     int to;  
  5.     int weight;  
  6.     Edge(int f, int t, int w):from(f), to(t), weight(w) {}  
  7. };  
  8. int num_nodes;  
  9. int num_edges;  
  10. vector<Edge> edges;  
  11. vector<int> G[max_nodes]; // 每個節點出發的邊編號
  12. int p[max_nodes]; // 當前節點單源最短路中的上一條邊
  13. int d[max_nodes]; // 單源最短路徑長
struct Edge
{
    int from;
    int to;
    int weight;

    Edge(int f, int t, int w):from(f), to(t), weight(w) {}
};

int num_nodes;
int num_edges;
vector<Edge> edges;
vector<int> G[max_nodes]; // 每個節點出發的邊編號
int p[max_nodes]; // 當前節點單源最短路中的上一條邊
int d[max_nodes]; // 單源最短路徑長

Dijkstra 方法

Dijkstra 方法依據其優先佇列的實現不同,可以寫成幾種時間複雜度不同的演算法。它是圖論-最短路中最經典、常見的演算法。關於這個方法,網上有許多分析,但是我最喜歡的還是《演算法概論》中的講解。為了理解 Dijkstra 方法,首先回顧一下無權最短路的演算法。無權最短路演算法基於 BFS,每次從源點向外擴充套件一層,並且給擴充套件到的頂點標明距離,這個距離就是最短路的長。我們完全可以仿照這個思路,把帶權圖最短路問題規約到無權圖最短路問題——只要把長度大於 1 的邊填充進一些「虛頂點」即可。如下圖所示。

Dijkstra

這個辦法雖然可行,但是顯然效率很低。不過,Dijkstra 方法EC,EB,ED分別出發,經過一系列「虛節點」,依次到達D,B,C 。為了不在虛節點處浪費時間,出發之前,我們設定三個鬧鐘,時間分別為4,3,2提醒我們預計在這些時刻會有重要的事情發生(經過實際節點)。更一般地說,假設現在我們處理到了某個頂點u,和u相鄰接的頂點為v1,v2,…,vn,它們和uu的距離為d1,d2,…,dn。我們為v1,v2,…,vn各設定一個鬧鐘。如果還沒有設定鬧鐘,那麼設定為d ;如果設定的時間比d晚,那麼重新設定為d(此時我們沿著這條路比之前的某一條路會提前趕到)。每次鬧鐘響起,都說明可能經過了實際節點,我們都會更新這些資訊,直到不存在任何鬧鐘。綜上所述,也就是隨著 BFS 的進行,我們一旦發現更近的路徑,就立即更新路徑長,直到處理完最後(最遠)的一個頂點。由此可見,由於上述「虛頂點」並非我們關心的實際頂點,因此 Dijkstra 方法的處理方式為:直接跳過了它們。

還需要解決的一個問題,就是鬧鐘的管理。鬧鐘一定是從早到晚按順序響起的,然而我們設鬧鐘的順序卻不一定按照時間升序,因此需要一個優先佇列來管理。Dijkstra 方法實現的效率嚴重依賴於優先佇列的實現。一個使用標準庫容器介面卡 priority_queue 的演算法版本如下:

[cpp] view plain copy print?
  1. typedef pair<intint> HeapNode;  
  2. void Dijkstra(int s)  
  3. {  
  4.     priority_queue< HeapNode, vector<HeapNode>, greater<HeapNode> > Q;  
  5.     for (int i=0; i<num_nodes; ++i)  
  6.         d[i] = __inf;  
  7.     d[s] = 0;  
  8.     Q.push(make_pair(0, s));  
  9.     while (!Q.empty()) {  
  10.         pair<intint> N = Q.top();  
  11.         Q.pop();  
  12.         int u = N.second;  
  13.         if (N.first != d[u]) continue;  
  14.         for (int i=0; i<G[u].size(); ++i) {  
  15.             Edge &e = edges[G[u][i]];  
  16.             if (d[e.to] > d[u] + e.weight) {  
  17.                 d[e.to] = d[u] + e.weight;  
  18.                 p[e.to] = G[u][i];  
  19.                 Q.push(make_pair(d[e.to], e.to));  
  20.             }  
  21.         }  
  22.     }  
  23. }  
typedef pair<int, int> HeapNode;
void Dijkstra(int s)
{
    priority_queue< HeapNode, vector<HeapNode>, greater<HeapNode> > Q;
    for (int i=0; i<num_nodes; ++i)
        d[i] = __inf;

    d[s] = 0;
    Q.push(make_pair(0, s));
    while (!Q.empty()) {
        pair<int, int> N = Q.top();
        Q.pop();
        int u = N.second;
        if (N.first != d[u]) continue;
        for (int i=0; i<G[u].size(); ++i) {
            Edge &e = edges[G[u][i]];
            if (d[e.to] > d[u] + e.weight) {
                d[e.to] = d[u] + e.weight;
                p[e.to] = G[u][i];
                Q.push(make_pair(d[e.to], e.to));
            }
        }
    }
}

Bellman-Ford:一個簡單的想法

Dijkstra 方法的本質是進行一系列如下的更新操作

d(v)=min{d(v),d(u)+l(u,v)}

然而,如果邊權含有負值,那麼 Dijkstra 方法將不再適用。原因解釋如下

假設最終的最短路徑為:

su1u2u3ukt

不難看出,如果按照 (s,u1),(u1,u2),,(uk,t) 的順序執行上述更新操作,最終t的最短路徑一定是正確的。而且,只要保證上述更新操作全部按順序執行即可,並不要求上述更新操作是連續進行的。Dijkstra 演算法所執行的更新序列是經過選擇的,而選擇基於這一假設:st的最短路一定不會經過s距離大於l(s,t)的點。對於正權圖這一假設是顯然的,對於負權圖這一假設是錯誤的。

因此,為了求出負權圖的最短路徑,我們需要保證一個合理的更新序列。但是,我們並不知道最終的最短路徑!因此一個簡單的想法就是:更新所有的邊,每條邊都更新∣V∣−1次。由於多餘的更新操作總是無害的,因此演算法(幾乎)可以正確執行。等等,為什麼是∣V∣−1次?這是由於,任何含有∣V∣個頂點的圖兩個點之間的最短路徑最多含有∣V∣−1條邊這意味著最短路不會包含環。理由是,如果是負環,最短路不存在;如果是正環,去掉後變短;如果是零環,去掉後不變。

演算法實現中唯一一個需要注意的問題就是負值圈 (negative-cost cycle)。負值圈指的是,權值總和為負的圈。如果存在這種圈,我們可以在裡面滯留任意長而不斷減小最短路徑長,因此這種情況下最短路徑可能是不存在的,可能使程式陷入無限迴圈。好在,本文介紹的幾種演算法都可以判斷負值圈是否存在。對於 Bellman-Ford 演算法來說,判斷負值圈存在的方法是:在V1次迴圈之後再執行一次迴圈,如果還有更新操作發生,則說明存在負值圈。

Bellman-Ford 演算法的程式碼如下:

[cpp] view plain copy print?
  1. bool Bellman_Ford(int s)  
  2. {  
  3.     for (int i=0; i<num_nodes; ++i)  
  4.         d[i] = __inf;  
  5.     d[s] = 0;  
  6.     for (int i=0; i<num_nodes; ++i) {  
  7.         bool changed = false;  
  8.         for (int e=0; e<num_edges; ++e) {  
  9.             if (d[edges[e].to] > d[edges[e].from] + edges[e].weight   
  10.                && d[edges[e].from] != __inf) {  
  11.                 d[edges[e].to] = d[edges[e].from] + edges[e].weight;  
  12.                 p[edges[e].to] = e;  
  13.                 changed = true;  
  14.             }  
  15.         }  
  16.         if (!changed) returntrue;  
  17.         if (i == num_nodes && changed) returnfalse;  
  18.     }  
  19.     returnfalse// 程式應該永遠不會執行到這裡
  20. }  
bool Bellman_Ford(int s)
{
    for (int i=0; i<num_nodes; ++i)
        d[i] = __inf;

    d[s] = 0;
    for (int i=0; i<num_nodes; ++i) {
        bool changed = false;
        for (int e=0; e<num_edges; ++e) {
            if (d[edges[e].to] > d[edges[e].from] + edges[e].weight 
               && d[edges[e].from] != __inf) {
                d[edges[e].to] = d[edges[e].from] + edges[e].weight;
                p[edges[e].to] = e;
                changed = true;
            }
        }
        if (!changed) return true;
        if (i == num_nodes && changed) return false;
    }
    return false; // 程式應該永遠不會執行到這裡
}

註記:

  1. 如果某次迴圈沒有更新操作發生,以後也不會有了。我們可以就此結束程式,避免無效的計算。
  2. 上述程式中第 11 行的判斷:如果去掉這個判斷,只要圖中存在負值圈函式就會返回 false。否則,僅在給定源點可以達到負值圈時才返回 false

SPFA:一個改進的想法

看來,Bellman-Ford 演算法多少有些浪費。這裡介紹的 SPFA 可以算作是 Bellman-Ford 的佇列改進版。在 Dijkstra 方法中,隨著 BFS 的進行,最短路一點一點地「生長」。然而如果存在負權,我們的演算法必須允許「繞回去」的情況發生。換句話說,我們需要在某些時候撤銷已經形成的最短路。同時,我們還要改變 Bellman-Ford 演算法盲目更新的特點,只更新有用的節點。SPFA 中,一開始,我們把源點 s放入佇列中,然後每次迴圈讓一個頂點u出隊,找出所有與u鄰接的頂點v,更新最短路,並當v不在佇列裡時將它入隊。迴圈直到佇列為空(沒有需要更新的頂點)。
可以顯示出 SPFA 和 Bellman-Ford 演算法相比的一個重大改進的最經典的一個例子,就是圖為一條鏈的情形。
容易看出,如果存在負值圈,這個演算法將無限迴圈下去。因此需要一個機制來確保演算法得以中止。由於最短路最長只含有∣V∣−1條邊,因此如果任何一個頂點已經出隊∣V∣+1次,演算法停止執行

SPFA 的程式碼如下:

[cpp] view plain copy print?
  1. bool in_queue[max_nodes];  
  2. int cnt[max_nodes];  
  3. bool SPFA(int s)  
  4. {  
  5.     int u;  
  6.     queue<int> Q;  
  7.     memset(in_queue, 0, sizeof(in_queue));  
  8.     memset(cnt, 0, sizeof(cnt));  
  9.     for (int i=0; i<num_nodes; ++i)  
  10.         d[i] = __inf;  
  11.     d[s] = 0;  
  12.     in_queue[s] = true;  
  13.     Q.push(s);  
  14.     while (!Q.empty()) {  
  15.         in_queue[u=Q.front()] = false;  
  16.         Q.pop();  
  17.         for (int i=0; i<G[u].size(); ++i) {  
  18.             Edge &e = edges[G[u][i]];  
  19.             if (d[e.to] > d[u] + e.weight) {  
  20.                 d[e.to] = d[u] + e.weight;  
  21.                 p[e.to] = G[u][i];  
  22.                 if (!in_queue[e.to]) {  
  23.                     Q.push(e.to);  
  24.                     in_queue[e.to] = true;  
  25.                     if (++cnt[e.to] > num_nodes)  
  26.                         returnfalse;  
  27.                 }  
  28.             }  
  29.         }  
  30.     }  
  31.     returntrue;  
  32. }  
bool in_queue[max_nodes];
int cnt[max_nodes];
bool SPFA(int s)
{
    int u;
    queue<int> Q;
    memset(in_queue, 0, sizeof(in_queue));
    memset(cnt, 0, sizeof(cnt));
    for (int i=0; i<num_nodes; ++i)
        d[i] = __inf;

    d[s] = 0;
    in_queue[s] = true;
    Q.push(s);
    while (!Q.empty()) {
        in_queue[u=Q.front()] = false;
        Q.pop();
        for (int i=0; i<G[u].size(); ++i) {
            Edge &e = edges[G[u][i]];
            if (d[e.to] > d[u] + e.weight) {
                d[e.to] = d[u] + e.weight;
                p[e.to] = G[u][i];
                if (!in_queue[e.to]) {
                    Q.push(e.to);
                    in_queue[e.to] = true;
                    if (++cnt[e.to] > num_nodes)
                        return false;
                }
            }
        }
    }
    return true;
}

我們已經給出 SPFA 的執行時間為O(kE)(kV)O。實際上,可以期望k<2。但是可以構造出使 SPFA 演算法變得很慢的針對性資料。

Acyclic Shortest Path

如果圖是無圈的(acyclic)(或稱為有向無環圖,DAG),那麼情況就變的容易了。我們可以構造一個拓撲升序序列,由拓撲排序的性質,無圈圖的任意路徑中,頂點都是沿著「拓撲升序序列」排列的,因此我們只需要按照這個序列執行更新操作即可。顯然,這樣可以線上性時間內解決問題。

實現上,拓撲排序和更新可以一趟完成。這種演算法的程式碼如下:

[cpp] view plain copy print?
  1. int indegree[max_nodes];  
  2. void asp(int s)  
  3. {  
  4.     queue<int> Q;  
  5.     for (int i=0; i<num_nodes; ++i) {  
  6.         d[i] = __inf;  
  7.         indegree[i] = 0;  
  8.     }  
  9.     for (int i=0; i<num_edges; ++i)  
  10.         ++indegree[edges[i].to];  
  11.     for (int i=0; i<num_nodes; ++i)  
  12.         if (indegree[i] == 0) Q.push(i);  
  13.     d[s] = 0;  
  14.     while (!Q.empty()) {  
  15.         int w = Q.front();  
  16.         Q.pop();  
  17.         for (int i=0; i<G[w].size(); ++i) {  
  18.             if (d[edges[G[w][i]].to] > d[w] + edges[G[w][i]].weight && d[w] != __inf) {   
  19.                 d[edges[G[w][i]].to] = d[w] + edges[G[w][i]].weight;  
  20.                 p[edges[G[w][i]].to] = G[w][i];  
  21.             }  
  22.             if (–indegree[edges[G[w][i]].to] == 0)  
  23.                 Q.push(edges[G[w][i]].to);  
  24.         }  
  25.     }  
  26. }  
  27. Floyd-Warshall  
int indegree[max_nodes];
void asp(int s)
{
    queue<int> Q;
    for (int i=0; i<num_nodes; ++i) {
        d[i] = __inf;
        indegree[i] = 0;
    }
    for (int i=0; i<num_edges; ++i)
        ++indegree[edges[i].to];
    for (int i=0; i<num_nodes; ++i)
        if (indegree[i] == 0) Q.push(i);

    d[s] = 0;
    while (!Q.empty()) {
        int w = Q.front();
        Q.pop();
        for (int i=0; i<G[w].size(); ++i) {
            if (d[edges[G[w][i]].to] > d[w] + edges[G[w][i]].weight && d[w] != __inf) { 
                d[edges[G[w][i]].to] = d[w] + edges[G[w][i]].weight;
                p[edges[G[w][i]].to] = G[w][i];
            }
            if (--indegree[edges[G[w][i]].to] == 0)
                Q.push(edges[G[w][i]].to);
        }
    }
}
Floyd-Warshall

Floyd-Warshall

與前面四種不同,Floyd-Warshall 演算法是所謂的「全源最短路徑」,也就是任意兩點間的最短路徑。它並不是對單源最短路徑V次迭代的一種漸進改進,但是對非常稠密的圖卻可能更快,因為它的迴圈更加緊湊。而且,這個演算法支援負的權值。

Floyd-Warshall 演算法如下:

[cpp] view plain copy print?
  1. int dist[max_nodes][max_nodes]; // 記錄路徑長
  2. int path[max_nodes][max_nodes]; // 記錄實際路徑
  3. bool Floyd_Warshall()  
  4. {  
  5.     for (int i=0; i<num_nodes; ++i)  
  6.         for (int j=0; j<num_nodes; ++j)  
  7.             path[i][j] = j;  
  8.     for (int k=0; k<num_nodes; ++k) {  
  9.         for (int i=0; i<num_nodes; ++i) {  
  10.             for (int j=0; j<num_nodes; ++j) {  
  11.                 if (dist[i][j] > dist[i][k] + dist[k][j]  
  12.                  && dist[i][k] != __inf && dist[k][j] != __inf) {  
  13.                     dist[i][j] = dist[i][k] + dist[k][j];  
  14.                     path[i][j] = path[i][k];  
  15.                     if (i == j && dist[i][j] < 0)  
  16.                         returnfalse;  
  17.                 }  
  18.             }  
  19.         }  
  20.     }  
  21.     returntrue;  
  22. }  
int dist[max_nodes][max_nodes]; // 記錄路徑長
int path[max_nodes][max_nodes]; // 記錄實際路徑
bool Floyd_Warshall()
{
    for (int i=0; i<num_nodes; ++i)
        for (int j=0; j<num_nodes; ++j)
            path[i][j] = j;

    for (int k=0; k<num_nodes; ++k) {
        for (int i=0; i<num_nodes; ++i) {
            for (int j=0; j<num_nodes; ++j) {
                if (dist[i][j] > dist[i][k] + dist[k][j]
                 && dist[i][k] != __inf && dist[k][j] != __inf) {
                    dist[i][j] = dist[i][k] + dist[k][j];
                    path[i][j] = path[i][k];
                    if (i == j && dist[i][j] < 0)
                        return false;
                }
            }
        }
    }
    return true;
}

其中 dist 陣列應初始化為鄰接矩陣。需要提醒的是, dist[i][i] 實際上表示「從頂點 i 繞一圈再回來的最短路徑」,因此圖存在負環當且僅當 dist[i][i]<0。初始化時, dist[i][i]可以初始化為0,也可以初始化為

顯示實際路徑

顯示實際路徑實際上非常簡單。利用前四個演算法提供的 int *p,還原實際路徑的一個方法如下:

[cpp] view plain copy print?
  1. void printpath(int from, int to, bool firstcall = true)  
  2. {  
  3.     if (from == to) {  
  4.         printf(”%d”, from);  
  5.         return;  
  6.     }  
  7.     if (p[to] == -1) return;  
  8.     if (firstcall) printf(“%d ->”, from);  
  9.     int v = edges[p[to]].from;  
  10.     if (v == from) {  
  11.         if (firstcall) printf(“ %d”, to);  
  12.         return;  
  13.     }  
  14.     printpath(from, v, false);  
  15.     printf(” %d ->”, v+1);  
  16.     if (firstcall) printf(“ %d”, to);  
  17. <span style=”font-size:14px;”>}</span>  
void printpath(int from, int to, bool firstcall = true)
{
    if (from == to) {
        printf("%d", from);
        return;
    }
    if (p[to] == -1) return;
    if (firstcall) printf("%d ->", from);
    int v = edges[p[to]].from;
    if (v == from) {
        if (firstcall) printf(" %d", to);
        return;
    }
    printpath(from, v, false);
    printf(" %d ->", v+1);
    if (firstcall) printf(" %d", to);
<span style="font-size:14px;">}</span>
利用 Floyd-Warshall 演算法提供的 int **path,還原實際路徑的一個方法如下:
[cpp] view plain copy print?
  1. void showpath(int from, int to)  
  2. {  
  3.     int c = from;  
  4.     printf(”%d -> %d:(%2d)  %d”, from, to, dist[from][to], from);  
  5.     while (c != to) {  
  6.         printf(” -> %d”, path[c][to]);  
  7.         c = path[c][to];  
  8.     }  
  9.     printf(”\n”);  
  10. }  
void showpath(int from, int to)
{
    int c = from;
    printf("%d -> %d:(%2d)  %d", from, to, dist[from][to], from);
    while (c != to) {
        printf(" -> %d", path[c][to]);
        c = path[c][to];
    }
    printf("\n");
}