[轉]從入門到精通: 最小費用流的“zkw算法”
>>>> 原文地址：最小費用流的“zkw算法” <<<<
1. 網絡流的一些基本概念
很多同學建立過網絡流模型做題目, 也學過了各種算法, 但是對於基本的概念反而說不清楚. 雖然不同的模型在具體叫法上可能不相同, 但是不同叫法對應的思想是一致的. 下面的討論力求規範, 個別地方可能需要對通常的叫法加以澄清.
求解可行流: 給定一個網絡流圖, 初始時每個節點不一定平衡 (每個節點可以有盈余或不足), 每條邊的流量可以有上下界, 每條邊的當前流量可以不滿足上下界約束. 可行流求解中沒有源和匯的概念, 算法的目的是尋找一個可以使所有節點都能平衡, 所有邊都能滿足流量約束的方案, 同時可能附加有最小費用的條件 (最小費用可行流).
求解最大流: 給定一個網絡流圖, 其中有兩個特殊的節點稱為源和匯. 除源和匯之外, 給定的每個節點一定平衡. 源可以產生無限大的流量, 匯可以吸收無限大的流量. 標準的最大流模型, 初始解一定是可行的 (例如, 所有邊流量均為零), 因此邊上不能有下界. 算法的目的是尋找一個從源到匯流量最大的方案, 同時不破壞可行約束, 並可能附加有最小費用的條件 (最小費用最大流).
擴展的最大流: 在有上下界或有節點盈余的網絡流圖中求解最大流. 實際上包括兩部分, 先是消除下界, 消除盈余, 可能還需要消除不滿足最優條件的流量 (最小費用流), 找到一個可行流, 再進一步得到最大流. 因此這裏我們的轉化似乎是從最大流轉化為可行流再變回最大流, 但其實質是將一個過程 (擴展的最大流) 變為了兩個過程 (可行流 + 最大流).
以上概念同時適用於最小費用流.
2. 最小費用流的各種轉化
不少同學認為消圈算法比增廣路算法使用面廣, 可以有負邊, 負圈, 每個點都能有盈余虧空等等. 實際上增廣路算法也一樣可以有負邊, 負圈, 上下界等等, 且一般來說速度快於消圈算法.
以下討論中 
為盈余量, 
為費用, 
為容量.
1.最小費用(可行)流 
最小費用最大流
建立超級源 
和超級匯 
, 對頂點 
, 若 
添加邊 
, 若 
添加邊 
, 之後求從 到 的最小費用最大流, 如果流量等於 
, 就存在可行流, 殘量網絡已在原圖上求出.
2.最小費用最大流 最小費用(可行)流
連邊 
, 所有點 有 
3.最小費用(可行)流中負權邊的消除
直接將負權邊滿流.
4.最小費用最大流中負權邊的消除
先連邊
, 使用 (3.) 中的方法消除負權邊, 使用 (1.) 中的方法求出最小費用 (可行) 流, 之後距離標號不變, 再求最小費用最大流; 註意此時增廣費用不能機械使用源點的標號——應該是源點匯點標號之差.
3. 費用流中的負邊和負圈
在費用流的求解中難免會遇到負邊和負圈的問題. 對於消圈算法, 負圈就是算法本身的一部分; 但對於增廣路算法, 負圈是我們所不願意看到的. 鑒於競賽中使用消圈算法將帶來嚴重的效率問題, 我們必須探索使用增廣路算法處理費用流負邊和負圈的方法.
對於單純的負邊, 如果這些負邊沒有組成負圈, 可以使用重賦權技術來處理, 即通過對每個節點合適的頂標, 使得 Reduced Cost 非負. 這個頂標通常可以使用到匯點的距離, 修改之後的邊權變為 
. 根據流量平衡條件, 容易根據這個新的費用算出原費用, 同時可以證明這樣的重賦權不改變最優解. 這樣做的代價是一次 SPFA 操作, 時間耗費是相對較小的.
如果存在負圈, SPFA 算法將不會終止. 很多同學使用人為的限制條件使其終止, 這是錯誤的. 舉一個例子: 源點和匯點均為孤立點, 圖中另外存在一個負費用, 正容量的圈. 不管怎樣初始標號, 僅憑增廣無法消除這個負圈, 從而得不到最優解 (根據最小費用最大流的定義, 要在所有最大流 (流量都為 
) 中尋找一個費用最小的方案). 也許你會說這個例子太過極端, 但是也很容易構造出其他一些例子, 說明不處理負圈, 或者簡單地對算法打補丁並不能解決本質問題.
解決方法就是將所有負費用的邊強制滿流, 稱為“退流”操作. 這樣做之後我們破壞了平衡條件, 但滿足了最優條件. 之後我們可以先求可行流, 再求最大流, 其間一直維持最優條件, 並逐步解決平衡條件, 最大流條件等問題. 這也就是 (2.) 中提到的方法. 這個方法可以消除所有負權邊 (在 Reduced Cost 意義下), 同時正確處理所有負圈問題. 那麽, 這個方法的速度如何呢?
費用流相關的圖大致可以分為兩類, 一類圖側重於費用, 即總的流量不大, 如 Two Shortest (只有 
), KaKa‘s Matrix Travels 等. 較為通用的強制滿流方法在這類圖上表現不佳, 原因是原本的流量較小而負費用較多, 經過轉化的圖在可行流求解階段流量與原流量相比很大, 嚴重影響速度. 另一方面, 這類圖往往有深刻的最短路背景, 從而不會出現負圈, 可以使用 SPFA 重賦權. 另一類圖側重於流量, 即邊的費用相差不大, 如 employee, 最優圖像的求解. 這種圖中流量不小而負費用相對有限 (也可以稍大, 關鍵是前者), 經過轉化增加的流量可以很快完成計算, 因此強制滿流方法就沒有問題.
另外, 不同的建圖方式有可能造成不同的模型. 在 employee 這道題中, 如果建圖時從每個點向後連邊, 從前向後運行最小費用最大流, 這時的算法就只有負邊, 沒有負圈. 而如果從每個點向前連邊, 在整個圖上求最小費用平衡流, 就會有負圈出現. 但是這裏並沒有本質區別: 將第二種模型中的所有向前連的邊強制滿流即可得到第一種模型, 如果將其他的負邊也強制滿流還能得到一個完全沒有負邊的模型. 大家可能也猜到了, 這幾種模型在速度上相去不遠, 畢竟根據我們剛才的討論, 增加的流量在算法的執行中不占主要地位.
4. 最小費用流的“zkw算法”
最小費用流在 OI 競賽中應當算是比較偏門的內容, 但是 NOI2008 中 employee 的突然出現確實讓許多人包括 zkw 自己措手不及. 可憐的 zkw 當時想出了最小費用流模型, 可是他從來沒有實現過, 所以不敢寫, 此題 0 分. zkw 現在對費用流的心得是: 雖然理論上難, 但是寫一個能 AC 題的費用流還算簡單. 先貼一個我寫的 costflow 程序: 只有不到 70 行, 費用流比最大流還好寫～
1 #include <cstdio> 2 #include <cstring> 3 using namespace std; 4 const int maxint=~0U>>1; 5 6 int n,m,pi1,cost=0; 7 bool v[550]; 8 struct etype 9 { 10 int t,c,u; 11 etype *next,*pair; 12 etype(){} 13 etype(int t_,int c_,int u_,etype* next_): 14 t(t_),c(c_),u(u_),next(next_){} 15 void* operator new(unsigned,void* p){return p;} 16 } *e[550]; 17 18 int aug(int no,int m) 19 { 20 if(no==n)return cost+=pi1*m,m; 21 v[no]=true; 22 int l=m; 23 for(etype *i=e[no];i;i=i->next) 24 if(i->u && !i->c && !v[i->t]) 25 { 26 int d=aug(i->t,l<i->u?l:i->u); 27 i->u-=d,i->pair->u+=d,l-=d; 28 if(!l)return m; 29 } 30 return m-l; 31 } 32 33 bool modlabel() 34 { 35 int d=maxint; 36 for(int i=1;i<=n;++i)if(v[i]) 37 for(etype *j=e[i];j;j=j->next) 38 if(j->u && !v[j->t] && j->c<d)d=j->c; 39 if(d==maxint)return false; 40 for(int i=1;i<=n;++i)if(v[i]) 41 for(etype *j=e[i];j;j=j->next) 42 j->c-=d,j->pair->c+=d; 43 pi1 += d; 44 return true; 45 } 46 47 int main() 48 { 49 freopen("costflow.in","r",stdin); 50 freopen("costflow.out","w",stdout); 51 scanf("%d %d",&n,&m); 52 etype *Pe=new etype[m+m]; 53 while(m--) 54 { 55 int s,t,c,u; 56 scanf("%d%d%d%d",&s,&t,&u,&c); 57 e[s]=new(Pe++)etype(t, c,u,e[s]); 58 e[t]=new(Pe++)etype(s,-c,0,e[t]); 59 e[s]->pair=e[t]; 60 e[t]->pair=e[s]; 61 } 62 do do memset(v,0,sizeof(v)); 63 while(aug(1,maxint)); 64 while(modlabel()); 65 printf("%d\n",cost); 66 return 0; 67 }
這裏使用的是連續最短路算法. 最短路算法? 為什麽程序裏沒有 SPFA? Dijkstra? 且慢, 先讓我們回顧一下圖論中最短路算法中的距離標號. 定義 
為點 的距離標號, 任何一個最短路算法保證, 算法結束時對任意指向頂點 、從頂點 
出發的邊滿足 
(條件1), 且對於每個 存在一個 使得等號成立 (條件2). 換句話說, 任何一個滿足以上兩個條件的算法都可以叫做最短路, 而不僅僅是 SPFA、Dijkstra, 算法結束後, 恰在最短路上的邊滿足 
.
在最小費用流的計算中, 我們每次沿 的路徑增廣後都不會破壞條件 1, 但是可能破壞了條件 2. 不滿足條件 2 的後果是什麽呢? 使我們找不到每條邊都滿足 
新的增廣路. 只好每次增廣後使用 Dijkstra, SPFA 等等算法重新計算新的滿足條件 2 的距離標號. 這無疑是一種浪費. KM 算法中我們可以修改不斷修改可行頂標, 不斷擴大可行子圖, 這裏也同樣, 我們可以在始終滿足條件 1 的距離標號上不斷修改, 直到可以繼續增廣 (滿足條件 2).
回顧一下 KM 算法修改頂標的方法. 根據最後一次尋找交錯路不成功的 DFS, 找到 
, 左邊的點增加 
, 右邊的點減少 . 這裏也一樣, 根據最後一次尋找增廣路不成功的 DFS, 找到 $ d = \min_{i \in V, j \notin V, u_{ij} > 0} \left\{ c_{ij} - D_i + D_j } \right\}$ , 所有訪問過的點距離標號增加 . 可以證明, 這樣不會破壞性質 1, 而且至少有一條新的邊進入了 的子圖.
算法的步驟就是初始標號設為 , 不斷增廣, 如果不能增廣, 修改標號繼續增廣, 直到徹底不能增廣: 源點的標號已經被加到了 
. 註意: 在程序中所有的 cost 均表示的是 reduced cost, 即 . 另外, 這個算法不能直接用於有任何負權邊的圖. 更不能用於負權圈的情況. 有關這兩種情況的處理, 參見 (2.) 和 (3.) 中的說明.
這樣我們得到了一個簡單的算法, 只需要增廣, 改標號, 各自只有 7 行, 不需要 BFS, 隊列, SPFA, 編程復雜度很低. 由於實際的增廣都是沿最短路進行的, 所以理論時間復雜度與使用 SPFA 等等方法的連續最短路算法一致, 但節省了 SPFA 或者 Dijkstra 的運算時間. 實測發現這種算法常數很小, 速度較快, employee 這道題所有數據加在一起耗時都在 2s 之內.
5. “zkw” 費用流算法在哪些圖上慢
實踐中, 上面的這個算法非常奇怪. 在某一些圖上, 算法速度非常快, 另一些圖上卻比純 SPFA 增廣的算法慢. 不少同學經過實測總結的結果是稠密圖上比較快, 稀疏圖上比較慢, 但也不盡然. 這裏我從理論上分析一下, 究竟這個算法用於哪些圖可以得到理想的效果.
先分析算法的增廣流程. 和 SPFA 直接算法相比, 由於同屬於沿最短路增廣的算法, 實際進行的增流操作並沒有太多的區別, 每次的增流路徑也大同小異. 因此不考慮多路增廣時, 增廣次數應該基本相同. 運行時間上主要的差異應當在於如何尋找增流路徑的部分.
那麽 zkw 算法的優勢在於哪裏呢? 與 SPFA 相比, KM 的重標號方式明顯在速度上占優, 每次只是一個對邊的掃描操作而已. 而 SPFA 需要維護較為復雜的標號和隊列操作, 同時為了修正標號, 需要不止一次地訪問某些節點, 速度會慢不少. 另外, 在 zkw 算法中, 增廣是多路進行的, 同時可能在一次重標號後進行多次增廣. 這個特點可以在許多路徑都費用相同的時候派上用場, 進一步減少了重標號的時間耗費.
下面想一想 zkw 算法的劣勢, 也就是 KM 重標號方式存在的問題. KM 重標號的主要問題就是, 不保證經過一次重標號之後能夠存在增廣路. 最差情況下, 一次只能在零權網絡中增加一條邊而已. 這時算法就會反復重標號, 反復嘗試增廣而次次不能增廣, 陷入弄巧成拙的境地.
接下來要說什麽, 大家可能已經猜到了. 對於最終流量較大, 而費用取值範圍不大的圖, 或者是增廣路徑比較短的圖 (如二分圖), zkw 算法都會比較快. 原因是充分發揮優勢. 比如流多說明可以同一費用反復增廣, 費用窄說明不用改太多距離標號就會有新增廣路, 增廣路徑短可以顯著改善最壞情況, 因為即使每次就只增加一條邊也可以很快湊成最短路. 如果恰恰相反, 流量不大, 費用不小, 增廣路還較長, 就不適合 zkw 算法了.
不適合怎麽辦? 繼續向下看.
6. 最小費用流的原始對偶 (Primal-Dual) 算法
最小費用流的直接 SPFA 算法和前面說的 KM 重標號算法, 各自都有一些情況非常慢. 這裏要寫的就是一個所謂的“新算法” (其實非常經典), 融合兩者優勢.
先從 zkw 算法的本質開始. 細心的同學已經註意到了, 算法的主要過程就是反復交替進行最短路和最大流的計算. 具體地說, 每一次增廣過程就是在零權網絡中找一個最大流增廣, 而重標號過程就是找到一組新的可行的距離標號. 既然如此, 最大流和最短路就完全可以用其他方法來替換. 這種交替進行最短路和最大流計算用來求得最小費用最大流的算法, 其實就是非常經典的原始對偶算法.
費用流的算法大致分為兩種, 一種是經典的解法, 如消圈, 增廣路, 原始對偶等等, 特點是步步為營, 維持可行性或者最優性其中之一, 再不斷對另一方面作出改進. 另一種就比較現代一些, 典型的例子是松弛算法和網絡單純形, 由於放松了對求解過程中解的限制條件, 使得其速度遠遠超過經典解法, 同時也增加了編程難度和理解障礙. 下面要說的原始對偶算法, 速度自然不可能比松弛和網絡單純形快, 但應該是經典解法中的佼佼者了.
先考慮最短路算法的選擇. 在流大, 費用小, 距離短的圖上, KM 重標號就可以取得不錯的效果, 但是 SPFA 在另外的圖上就很有優勢.
再考慮最大流算法的選擇. 最樸素的 SPFA 暴力增廣中根本就不維護距離標號, 這裏隱含地使用了一個單路增廣. 在 zkw 費用流的算法裏, 這裏使用的是多路增廣. 是不是有必要用更高級的網絡流算法呢? 經過實際測試, 除非在特別稠密的圖上, 否則效果一般. 原因是可增廣的流量並不是很大, 如果去維護 SAP 的頂標反而增加了不少開銷, 預流推進亦然.
總結一下, 折中的選擇就是使用 Small Label First 優化 的 SPFA 來維護 zkw 算法中的距離標號, 保留多路增廣. 註意, 除了第一次之外, SPFA 算法都工作在一個所有邊均為正權的圖上 (Reduced cost 意義下), 這是和不維護頂標直接 SPFA 的不同之處, 也是原來 zkw 算法的精神所在. 也正因為如此, 這裏的 SPFA 甚至可以用 Dijkstra 替換.
這個特殊的原始對偶算法在稠密二分費用小的圖上不敵原來的 zkw 算法, 但遠遠勝過暴力 SPFA. 在另外的圖上, 對兩者都是穩勝. 比暴力 SPFA 快原因是, 多路增廣, 同時使用了 Reduced Cost 縮小了費用範圍, 從而利於 SPFA 算法的工作 (需要的松弛次數減少), 而且使用 Reduced Cost 後不再有負邊, 使 SLF 的優化落到了實處 (回憶: SLF 優化只有當所有邊均為正的時候才能發揮出最佳效果), 甚至允許用 Dijkstra 來完成後面的工作. 比 zkw 算法快的原因是, 在流小費用大距離長的圖上, 一次性把距離標號改對往往比反復調整更有效率.
最後是代碼, 這次的例子是 POJ 3680 zkw Accepted 216K 360MS. 另外還有一個 POJ 3422, 是原來的 zkw 算法不擅長的一個例子, 這種方法只用 47MS, 代碼大同小異, 就不放在這裏了.
1 #include <cstdio> 2 #include <cstring> 3 #include <deque> 4 #include <algorithm> 5 using namespace std; 6 7 const int V=440, E=V*2, maxint=0x3F3F3F3F; 8 9 struct etype 10 { 11 int t, c, u; 12 etype *next, *pair; 13 etype() {} 14 etype(int T, int C, int U, etype* N): t(T), c(C), u(U), next(N) {} 15 void* operator new(unsigned, void* p){return p;} 16 } *e[V], Te[E+E], *Pe; 17 18 int S, T, n, piS, cost; 19 bool v[V]; 20 21 void addedge(int s, int t, int c, int u) 22 { 23 e[s] = new(Pe++) etype(t, +c, u, e[s]); 24 e[t] = new(Pe++) etype(s, -c, 0, e[t]); 25 e[s]->pair = e[t]; 26 e[t]->pair = e[s]; 27 } 28 29 int aug(int no, int m) 30 { 31 if (no == T) return cost += piS * m, m; 32 v[no] = true; 33 int l = m; 34 for (etype *i = e[no]; i; i = i->next) 35 if (i->u && !i->c && !v[i->t]) 36 { 37 int d = aug(i->t, l < i->u ? l : i->u); 38 i->u -= d, i->pair->u += d, l -= d; 39 if (!l) return m; 40 } 41 return m - l; 42 } 43 44 bool modlabel() 45 { 46 static int d[V]; memset(d, 0x3F, sizeof(d)); d[T] = 0; 47 static deque<int> Q; Q.push_back(T); 48 while(Q.size()) 49 { 50 int dt, no = Q.front(); Q.pop_front(); 51 for(etype *i = e[no]; i; i = i->next) 52 if(i->pair->u && (dt = d[no] - i->c) < d[i->t]) 53 (d[i->t] = dt) <= d[Q.size() ? Q.front() : 0] 54 ? Q.push_front(i->t) : Q.push_back(i->t); 55 } 56 for(int i = 0; i < n; ++i) 57 for(etype *j = e[i]; j; j = j->next) 58 j->c += d[j->t] - d[i]; 59 piS += d[S]; 60 return d[S] < maxint; 61 } 62 63 int ab[V], *pab[V], w[V]; 64 65 struct lt 66 { 67 bool operator()(int* p1,int* p2) {return *p1 < *p2;} 68 }; 69 70 int main() 71 { 72 int t; 73 scanf("%d",&t); 74 while(t--) 75 { 76 memset(e,0,sizeof(e)); 77 Pe = Te; 78 static int m, k; 79 scanf("%d %d", &m, &k); 80 int abz = 0; 81 for(int i = 0; i < m; ++i) 82 { 83 scanf("%d", pab[abz] = &ab[abz]), abz++; 84 scanf("%d", pab[abz] = &ab[abz]), abz++; 85 scanf("%d", &w[i]); 86 } 87 sort(&pab[0], &pab[abz], lt()); 88 int c=0xDEADBEEF; n=0; 89 for(int i = 0; i < abz; ++i) 90 { 91 if(c != *pab[i]) c = *pab[i], ++n; 92 *pab[i] = n; 93 } 94 ++n, S = 0, T = n++; 95 for(int i = 0; i < T; ++i) addedge(i, i+1, 0, k); 96 for(int i = 0; i < m; ++i) addedge(ab[i+i], ab[i+i+1], -w[i], 1); 97 piS = cost = 0; 98 while(modlabel()) 99 do memset(v, 0, sizeof(v)); 100 while(aug(S, maxint)); 101 printf("%d\n", -cost); 102 } 103 return 0; 104 }
[轉]從入門到精通: 最小費用流的“zkw算法”