最大流

本頁面主要介紹最大流問題相關的算法知識。

概述

網絡流基本概念參見網絡流簡介。

令 \(G=(V,E)\) 是一個有源匯點的網絡，我們希望在 \(G\) 上指定合適的流 \(f\)，以最大化整個網絡的流量 \(|f|\)（即 \(\sum_{x \in V} f(s, x) - \sum_{x \in V} f(x, s)\)），這一問題被稱作最大流問題（Maximum flow problem）。

Ford–Fulkerson 增廣

Ford–Fulkerson 增廣是計算最大流的一類算法的總稱。該方法運用貪心的思想，通過尋找增廣路來更新並求解最大流。

概述

給定網絡 \(G\) 及 \(G\) 上的流 \(f\)，我們做如下定義。

對於邊 \((u, v)\)，我們將其容量與流量之差稱為剩餘容量 \(c_f(u,v)\)（Residual Capacity），即 \(c_f(u,v)=c(u,v)-f(u,v)\)。

我們將 \(G\) 中所有結點和剩餘容量大於 \(0\) 的邊構成的子圖稱為殘量網絡 \(G_f\)（Residual Network），即 \(G_f=(V,E_f)\)，其中 \(E_f=\left\{(u,v) \mid c_f(u,v)>0\right\}\)。

Warning

正如我們馬上要提到的，流量可能是負值，因此，\(E_f\) 的邊有可能並不在 \(E\) 中。引入增廣的概念後，下文將具體解釋這一點。

我們將 \(G_f\) 上一條從源點 \(s\) 到匯點 \(t\) 的路徑稱為增廣路（Augmenting Path）。對於一條增廣路，我們給每一條邊 \((u, v)\) 都加上等量的流量，以令整個網絡的流量增加，這一過程被稱為增廣（Augment）。由此，最大流的求解可以被視為若干次增廣分別得到的流的疊加。

此外，在 Ford–Fulkerson 增廣的過程中，對於每條邊 \((u, v)\)，我們都新建一條反向邊 \((v, u)\)。我們約定 \(f(u, v) = -f(v, u)\)，這一性質可以通過在每次增廣時引入退流操作來保證，即 \(f(u, v)\) 增加時 \(f(v, u)\) 應當減少同等的量。

Tip

在最大流算法的代碼實現中，我們往往需要支持快速訪問反向邊的操作。在鄰接矩陣中，這一操作是 trivial 的（\(g_{u, v} \leftrightarrow g_{v, u}\)）。但主流的實現是更加優秀的鏈式前向星。其中，一個常用的技巧是，我們令邊從偶數（通常為 \(0\)）開始編號，並在加邊時總是緊接着加入其反向邊使得它們的編號相鄰。由此，我們可以令編號為 \(i\) 的邊和編號為 \(i \oplus 1\) 的邊始終保持互為反向邊的關係。

初次接觸這一方法的讀者可能察覺到一個違反直覺的情形——反向邊的流量 \(f(v, u)\) 可能是一個負值。實際上我們可以注意到，在 Ford–Fulkerson 增廣的過程中，真正有意義的是剩餘容量 \(c_f\)，而 \(f(v, u)\) 的絕對值是無關緊要的，我們可以將反向邊流量的減少視為反向邊剩餘容量 \(c_f(v, u)\) 的增加——這也與退流的意義相吻合——反向邊剩餘容量的增加意味着我們接下來可能通過走反向邊來和原先正向的增廣抵消，代表一種「反悔」的操作。

以下案例有可能幫助你理解這一過程。假設 \(G\) 是一個單位容量的網絡，我們考慮以下過程：

\(G\) 上有多條增廣路，其中，我們選擇進行一次先後經過 \(u, v\) 的增廣（如左圖所示），流量增加 \(1\)。
我們注意到，如果進行中圖上的增廣，這個局部的最大流量不是 \(1\) 而是 \(2\)。但由於指向 \(u\) 的邊和從 \(v\) 出發的邊在第一次增廣中耗盡了容量，此時我們無法進行中圖上的增廣。這意味着我們當前的流是不夠優的，但局部可能已經沒有其他（只經過原圖中的邊而不經過反向邊的）增廣路了。
現在引入退流操作。第一次增廣後，退流意味着 \(c_f(v, u)\) 增加了 \(1\) 剩餘容量，即相當於新增 \((v, u)\) 這條邊，因此我們可以再進行一次先後經過 \(p, v, u, q\) 的增廣（如右圖橙色路徑所示）。無向邊 \((u, v)\) 上的流量在兩次增廣中抵消，我們驚奇地發現兩次增廣疊加得到的結果實際上和中圖是等價的。

以上案例告訴我們，退流操作帶來的「抵消」效果使得我們無需擔心我們按照「錯誤」的順序選擇了增廣路。

容易發現，只要 \(G_f\) 上存在增廣路，那麼對其增廣就可以令總流量增加；否則説明總流量已經達到最大可能值，求解過程完成。這就是 Ford–Fulkerson 增廣的過程。

最大流最小割定理

我們大致瞭解了 Ford–Fulkerson 增廣的思想，可是如何證明這一方法的正確性呢？為什麼增廣結束後的流 \(f\) 是一個最大流？

實際上，Ford–Fulkerson 增廣的正確性和最大流最小割定理（The Maxflow-Mincut Theorem）等價。這一定理指出，對於任意網絡 \(G = (V, E)\)，其上的最大流 \(f\) 和最小割 \(\{S, T\}\) 總是滿足 \(|f| = ||S, T||\)。

為了證明最大流最小割定理，我們先從一個引理出發：對於網絡 \(G = (V, E)\)，任取一個流 \(f\) 和一個割 \(\{S, T\}\)，總是有 \(|f| \leq ||S, T||\)，其中等號成立當且僅當 \(\{(u, v) | u \in S, v \in T\}\) 的所有邊均滿流，且 \(\{(u, v) | u \in T, v \in S\}\) 的所有邊均空流。

證明

\[ \begin{aligned} |f| & = f(s) \\ & = \sum_{u \in S} f(u) \\ & = \sum_{u \in S} \left( \sum_{v \in V} f(u, v) - \sum_{v \in V} f(v, u) \right) \\ & = \sum_{u \in S} \left( \sum_{v \in T} f(u, v) + \sum_{v \in S} f(u, v) - \sum_{v \in T} f(v, u) - \sum_{v \in S} f(v, u) \right) \\ & = \sum_{u \in S} \left( \sum_{v \in T} f(u, v) - \sum_{v \in T} f(v, u) \right) + \sum_{u \in S} \sum_{v \in S} f(u, v) - \sum_{u \in S} \sum_{v \in S} f(v, u) \\ & = \sum_{u \in S} \left( \sum_{v \in T} f(u, v) - \sum_{v \in T} f(v, u) \right) \\ & \leq \sum_{u \in S} \sum_{v \in T} f(u, v) \\ & \leq \sum_{u \in S} \sum_{v \in T} c(u, v) \\ & = ||S, T|| \\ \end{aligned} \]

為了取等，第一個不等號需要 \(\{(u, v) \mid u \in T, v \in S\}\) 的所有邊均空流，第二個不等號需要 \(\{(v, u) \mid u \in S, v \in T\}\) 的所有邊均滿流。原引理得證。

那麼，對於任意網絡，以上取等條件是否總是能被滿足呢？如果答案是肯定的，則最大流最小割定理得證。以下我們嘗試證明。

證明

假設某一輪增廣後，我們得到流 \(f\) 使得 \(G_f\) 上不存在增廣路，即 \(G_f\) 上不存在 \(s\) 到 \(t\) 的路徑。此時我們記從 \(s\) 出發可以到達的結點組成的點集為 \(S\)，並記 \(T = V \setminus S\)。

顯然，\(\{S, T\}\) 是 \(G_f\) 的一個割，且 \(||S, T|| = \sum_{u \in S} \sum_{v \in T} c_f(u, v) = 0\)。由於剩餘容量是非負的，這也意味着對於任意 \(u \in S, v \in T, (u, v) \in E_f\)，均有 \(c_f(u, v) = 0\)。以下我們將這些邊分為存在於原圖中的邊和反向邊兩種情況討論：

\((u, v) \in E\)：此時，\(c_f(u, v) = c(u, v) - f(u, v) = 0\)，因此有 \(c(u, v) = f(u, v)\)，即 \(\{(u, v) \mid u \in S, v \in T\}\) 的所有邊均滿流；
\((v, u) \in E\)：此時，\(c_f(u, v) = c(u, v) - f(u, v) = 0 - f(u, v) = f(v, u) = 0\)，即 \(\{(v, u) \mid u \in S, v \in T\}\) 的所有邊均空流。

因此，增廣停止後，上述流 \(f\) 滿足取等條件。根據引理指出的大小關係，自然地，\(f\) 是 \(G\) 的一個最大流，\(\{S, T\}\) 是 \(G\) 的一個最小割。

容易看出，Kőnig 定理是最大流最小割定理的特殊情形。實際上，它們都和線性規劃中的對偶有關。

時間複雜度分析

在整數流量的網絡 \(G = (V, E)\) 上，平凡地，我們假設每次增廣的流量都是整數，則 Ford–Fulkerson 增廣的時間複雜度的一個上界是 \(O(|E||f|)\)，其中 \(f\) 是 \(G\) 上的最大流。這是因為單輪增廣的時間複雜度是 \(O(|E|)\)，而增廣會導致總流量增加，故增廣輪數不可能超過 \(|f|\)。

對於 Ford–Fulkerson 增廣的不同實現，時間複雜度也各不相同。其中較主流的實現有 Edmonds–Karp, Dinic, SAP, ISAP 等算法，我們將在下文中分別介紹。

Edmonds–Karp 算法

算法思想

如何在 \(G_f\) 中尋找增廣路呢？當我們考慮 Ford–Fulkerson 增廣的具體實現時，最自然的方案就是使用 BFS。此時，Ford–Fulkerson 增廣表現為 Edmonds–Karp 算法。其具體流程如下：

如果在 \(G_f\) 上我們可以從 \(s\) 出發 BFS 到 \(t\)，則我們找到了新的增廣路。
對於增廣路 \(p\)，我們計算出 \(p\) 經過的邊的剩餘容量的最小值 \(\Delta = \min_{(u, v) \in p} c_f(u, v)\)。我們給 \(p\) 上的每條邊都加上 \(\Delta\) 流量，並給它們的反向邊都退掉 \(\Delta\) 流量，令最大流增加了 \(\Delta\)。
因為我們修改了流量，所以我們得到新的 \(G_f\)，我們在新的 \(G_f\) 上重複上述過程，直至增廣路不存在，則流量不再增加。

以上算法即 Edmonds–Karp 算法。

時間複雜度分析

接下來讓我們嘗試分析 Edmonds–Karp 算法的時間複雜度。

顯然，單輪 BFS 增廣的時間複雜度是 \(O(|E|)\)。

增廣總輪數的上界是 \(O(|V||E|)\)。這一論斷在網絡資料中常被偽證（或被含糊其辭略過）。以下我們嘗試給出一個較正式的證明¹。

增廣總輪數的上界的證明

首先，我們引入一個引理——最短路非遞減引理。具體地，我們記 \(d_f(u)\) 為 \(G_f\) 上結點 \(u\) 到源點 \(s\) 的距離（即最短路長度，下同）。對於某一輪增廣，我們用 \(f\) 和 \(f'\) 分別表示增廣前的流和增廣後的流，我們斷言，對於任意結點 \(u\)，增廣總是使得 \(d_{f'}(u) \geq d_f(u)\)。我們將在稍後證明這一引理。

不妨稱增廣路上剩餘容量最小的邊是飽和邊（存在多條邊同時最小則取任一）。如果一條有向邊 \((u, v)\) 被選為飽和邊，增廣會清空其剩餘容量導致飽和邊的消失，並且退流導致反向邊的新增（如果原先反向邊不存在），即 \((u, v) \not \in E_{f'}\) 且 \((v, u) \in E_{f'}\)。以上分析使我們知道，對於無向邊 \((u, v)\)，其被增廣的兩種方向總是交替出現。

在 \(G_f\) 上沿 \((u, v)\) 增廣時，\(d_f(u) + 1 = d_f(v)\)，此後殘量網絡變為 \(G_{f'}\)。在 \(G_{f'}\) 上沿 \((v, u)\) 增廣時，\(d_{f'}(v) + 1 = d_{f'}(u)\)。根據最短路非遞減引理又有 \(d_{f'}(v) \geq d_f(v)\)，我們連接所有式子，得到 \(d_{f'}(u) \geq d_{f}(u) + 2\)。換言之，如果有向邊 \((u, v)\) 被選為飽和邊，那麼與其上一次被選為飽和邊時相比，\(u\) 到 \(s\) 的距離至少增加 \(2\)。

\(s\) 到任意結點的距離不可能超過 \(|V|\)，結合上述性質，我們發現每條邊被選為飽和邊的次數是 \(O(|V|)\) 的，與邊數相乘後得到增廣總輪數的上界 \(O(|V||E|)\)。

接下來我們證明最短路非遞減引理，即 \(d_{f'}(u) \geq d_f(u)\)。這一證明並不難，但可能稍顯繞口，讀者可以停下來認真思考片刻。

最短路非遞減引理的證明

考慮反證。對於某一輪增廣，我們假設存在若干結點，它們在該輪增廣後到 \(s\) 的距離較增廣前減小。我們記 \(v\) 為其中到 \(s\) 的距離最小的一者（即 \(v = \arg \min_{x \in V, d_{f'}(x) < d_f(x)} d_{f'}(x)\)）。注意，根據反證假設，此時 \(d_{f'}(v) < d_f(v)\) 是已知條件。

在 \(G_{f'}\) 中 \(s\) 到 \(v\) 的最短路上，我們記 \(u\) 是 \(v\) 的上一個結點，即 \(d_{f'}(u) + 1 = d_{f'}(v)\)。

為了不讓 \(u\) 破壞 \(v\) 的「距離最小」這一性質，\(u\) 必須滿足 \(d_{f'}(u) \geq d_f(u)\)。

對於上式，我們令不等號兩側同加，得 \(d_{f'}(v) \geq d_f(u) + 1\)。根據反證假設進行放縮，我們得到 \(d_f(v) > d_f(u) + 1\)。

以下我們嘗試討論 \((u, v)\) 上的增廣方向。

假設有向邊 \((u, v) \in E_f\)。根據 BFS「廣度優先」的性質，我們有 \(d_f(u) + 1 \geq d_f(v)\)。該式與放縮結果衝突，導出矛盾。
假設有向邊 \((u, v) \not \in E_f\)。根據 \(u\) 的定義我們已知 \((u, v) \in E_{f'}\)，因此這條邊的存在必須是當前輪次的增廣經過了 \((v, u)\) 並退流產生反向邊的結果，也即 \(d_f(v) + 1 = d_f(u)\)。該式與放縮結果衝突，導出矛盾。

由於 \((u, v)\) 沿任何方向增廣都會導出矛盾，我們知道反證假設不成立，最短路非遞減引理得證。

將單輪 BFS 增廣的複雜度與增廣輪數的上界相乘，我們得到 Edmonds–Karp 算法的時間複雜度是 \(O(|V||E|^2)\)。

代碼實現

Edmonds–Karp 算法的可能實現如下。

參考代碼

#define maxn 250
#define INF 0x3f3f3f3f

struct Edge {
  int from, to, cap, flow;

  Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};

struct EK {
  int n, m;             // n：點數，m：邊數
  vector<Edge> edges;   // edges：所有邊的集合
  vector<int> G[maxn];  // G：點 x -> x 的所有邊在 edges 中的下標
  int a[maxn], p[maxn];  // a：點 x -> BFS 過程中最近接近點 x 的邊給它的最大流
                         // p：點 x -> BFS 過程中最近接近點 x 的邊

  void init(int n) {
    for (int i = 0; i < n; i++) G[i].clear();
    edges.clear();
  }

  void AddEdge(int from, int to, int cap) {
    edges.push_back(Edge(from, to, cap, 0));
    edges.push_back(Edge(to, from, 0, 0));
    m = edges.size();
    G[from].push_back(m - 2);
    G[to].push_back(m - 1);
  }

  int Maxflow(int s, int t) {
    int flow = 0;
    for (;;) {
      memset(a, 0, sizeof(a));
      queue<int> Q;
      Q.push(s);
      a[s] = INF;
      while (!Q.empty()) {
        int x = Q.front();
        Q.pop();
        for (int i = 0; i < G[x].size(); i++) {  // 遍歷以 x 作為起點的邊
          Edge& e = edges[G[x][i]];
          if (!a[e.to] && e.cap > e.flow) {
            p[e.to] = G[x][i];  // G[x][i] 是最近接近點 e.to 的邊
            a[e.to] =
                min(a[x], e.cap - e.flow);  // 最近接近點 e.to 的邊賦給它的流
            Q.push(e.to);
          }
        }
        if (a[t]) break;  // 如果匯點接受到了流，就退出 BFS
      }
      if (!a[t])
        break;  // 如果匯點沒有接受到流，説明源點和匯點不在同一個連通分量上
      for (int u = t; u != s;
           u = edges[p[u]].from) {  // 通過 u 追尋 BFS 過程中 s -> t 的路徑
        edges[p[u]].flow += a[t];      // 增加路徑上邊的 flow 值
        edges[p[u] ^ 1].flow -= a[t];  // 減小反向路徑的 flow 值
      }
      flow += a[t];
    }
    return flow;
  }
};

Dinic 算法

算法思想

考慮在增廣前先對 \(G_f\) 做 BFS 分層，即根據結點 \(u\) 到源點 \(s\) 的距離 \(d(u)\) 把結點分成若干層。令經過 \(u\) 的流量只能流向下一層的結點 \(v\)，即刪除 \(u\) 向層數標號相等或更小的結點的出邊，我們稱 \(G_f\) 剩下的部分為層次圖（Level Graph）。形式化地，我們稱 \(G_L = (V, E_L)\) 是 \(G_f = (V, E_f)\) 的層次圖，其中 \(E_L = \left\{ (u, v) \mid (u, v) \in E_f, d(u) + 1 = d(v) \right\}\)。

如果我們在層次圖 \(G_L\) 上找到一個最大的增廣流 \(f_b\)，使得僅在 \(G_L\) 上是不可能找出更大的增廣流的，則我們稱 \(f_b\) 是 \(G_L\) 的阻塞流（Blocking Flow）。

Warning

儘管在上文中我們僅在單條增廣路上定義了增廣/增廣流，廣義地，「增廣」一詞不僅可以用於單條路徑上的增廣流，也可以用於若干增廣流的並——後者才是我們定義阻塞流時使用的意義。

定義層次圖和阻塞流後，Dinic 算法的流程如下。

在 \(G_f\) 上 BFS 出層次圖 \(G_L\)。
在 \(G_L\) 上 DFS 出阻塞流 \(f_b\)。
將 \(f_b\) 併到原先的流 \(f\) 中，即 \(f \leftarrow f + f_b\)。
重複以上過程直到不存在從 \(s\) 到 \(t\) 的路徑。

此時的 \(f\) 即為最大流。

在分析這一算法的複雜度之前，我們需要特別説明「在 \(G_L\) 上 DFS 出阻塞流 \(f_b\)」的過程。儘管 BFS 層次圖對於本頁面的讀者應當是 trivial 的，但 DFS 阻塞流的過程則稍需技巧——我們需要引入當前弧優化。

注意到在 \(G_L\) 上 DFS 的過程中，如果結點 \(u\) 同時具有大量入邊和出邊，並且 \(u\) 每次接受來自入邊的流量時都遍歷出邊表來決定將流量傳遞給哪條出邊，則 \(u\) 這個局部的時間複雜度最壞可達 \(O(|E|^2)\)。為避免這一缺陷，如果某一時刻我們已經知道邊 \((u, v)\) 已經增廣到極限（邊 \((u, v)\) 已無剩餘容量或 \(v\) 的後側已增廣至阻塞），則 \(u\) 的流量沒有必要再嘗試流向出邊 \((u, v)\)。據此，對於每個結點 \(u\)，我們維護 \(u\) 的出邊表中第一條還有必要嘗試的出邊。習慣上，我們稱維護的這個指針為當前弧，稱這個做法為當前弧優化。

多路增廣

多路增廣是 Dinic 算法的一個常數優化——如果我們在層次圖上找到了一條從 \(s\) 到 \(t\) 的增廣路 \(p\)，則接下來我們未必需要重新從 \(s\) 出發找下一條增廣路，而可能從 \(p\) 上最後一個仍有剩餘容量的位置出發尋找一條岔路進行增廣。考慮到其與回溯形式的一致性，這一優化在 DFS 的代碼實現中也是自然的。

常見誤區

可能是由於大量網絡資料的錯誤表述引發以訛傳訛的情形，相當數量的選手喜歡將當前弧優化和多路增廣並列稱為 Dinic 算法的兩種優化。實際上，當前弧優化是用於保證 Dinic 時間複雜度正確性的一部分，而多路增廣只是一個不影響複雜度的常數優化。

時間複雜度分析

應用當前弧優化後，對 Dinic 算法的時間複雜度分析如下。

首先，我們嘗試證明單輪增廣中 DFS 求阻塞流的時間複雜度是 \(O(|V||E|)\)。

單輪增廣的時間複雜度的證明

考慮阻塞流 \(f_b\) 中的每條增廣路，它們都是在 \(G_L\) 上每次沿當前弧跳轉而得到的結果，其中每條增廣路經歷的跳轉次數不可能多於 \(|V|\)。

每找到一條增廣路就有一條飽和邊消失（剩餘容量清零）。考慮阻塞流 \(f_b\) 中的每條增廣路，我們將被它們清零的飽和邊形成的邊集記作 \(E_1\)。考慮到 \(G_L\) 分層的性質，飽和邊消失後其反向邊不可能在同一輪增廣內被其他增廣路經過，因此，\(E_1\) 是 \(E_L\) 的子集。

此外，對於沿當前弧跳轉但由於某個位置阻塞所以沒有成功得到增廣路的情形，我們將這些不完整的路徑上的最後一條邊形成的邊集記作 \(E_2\)。\(E_2\) 的成員不飽和，所以 \(E_1\) 與 \(E_2\) 不交，且 \(E_1 \cup E_2\) 仍是 \(E_L\) 的子集。

由於 \(E_1 \cup E_2\) 的每個成員都沒有花費超過 \(|V|\) 次跳轉（且在使用多路增廣優化後一些跳轉將被重複計數），因此，綜上所述，DFS 過程中的總跳轉次數不可能多於 \(|V||E_L|\)。

常見偽證一則

對於每個結點，我們維護下一條可以增廣的邊，而當前弧最多變化 \(|E|\) 次，從而單輪增廣的最壞時間複雜度為 \(O(|V||E|)\)。

Bug

「當前弧最多變化 \(|E|\) 次」並不能推得「每個結點最多訪問其出邊 \(|E|\) 次」。這是因為，訪問當前弧並不一定耗盡上面的剩餘容量，結點 \(u\) 可能多次訪問同一條當前弧。

注意到層次圖的層數顯然不可能超過 \(|V|\)，如果我們可以證明層次圖的層數在增廣過程中嚴格單增，則 Dinic 算法的增廣輪數是 \(O(|V|)\) 的。接下來我們嘗試證明這一結論²。

層次圖層數單調性的證明

我們需要引入預流推進類算法（另一類最大流算法）中的一個概念——高度標號。為了更方便地結合高度標號表述我們的證明，在證明過程中，我們令 \(d_f(u)\) 為 \(G_f\) 上結點 \(u\) 到匯點 \(t\) 的距離，從匯點而非源點出發進行分層（這並沒有本質上的區別）。對於某一輪增廣，我們用 \(f\) 和 \(f'\) 分別表示增廣前的流和增廣後的流。在該輪增廣中求解並加入阻塞流後，記層次圖由 \(G_L = (V, E_L)\) 變為 \(G'_{L} = (V, E'_L)\)。

我們給高度標號一個不嚴格的臨時定義——在網絡 \(G = (V, E)\) 上，令 \(h\) 是點集 \(V\) 到整數集 \(N\) 上的函數，\(h\) 是 \(G\) 上合法的高度標號當且僅當 \(h(u) \leq h(v) + 1\) 對於 \((u, v) \in E\) 恆成立。

考察所有 \(E_{f'}\) 的成員 \((u, v)\)，我們發現 \((u, v) \in E_{f'}\) 的原因是以下二者之一。

\((u, v) \in E_f\)，且剩餘容量在該輪增廣過程中未耗盡——根據最短路的定義，此時我們有 \(d_f(u) \leq d_f(v) + 1\)；
\((u, v) \not \in E_f\)，但在該輪增廣過程中阻塞流經過 \((v, u)\) 並退流產生反向邊——根據層次圖和阻塞流的定義，此時我們有 \(d_f(u) + 1 = d_f(v)\)。

以上觀察讓我們得出一個結論——\(d_f\) 在 \(G_{f'}\) 上是一個合法的高度標號。當然，在 \(G_{f'}\) 的子圖 \(G'_L\) 上也是。

現在，對於一條 \(G'_L\) 上的增廣路 \(p = (s, \dots, u, v, \dots, t)\)，按照 \(p\) 上結點的反序（從 \(t\) 到 \(s\) 的順序）考慮從空路徑開始每次添加一個結點的過程。假設結點 \(v\) 已加入，結點 \(u\) 正在加入，我們發現，加入結點 \(u\) 後，根據層次圖的定義，\(d_{f'}(u)\) 的值較 \(d_{f'}(v)\) 增加 \(1\)；與此同時，由於 \(d_f\) 是 \(G'_L\) 上的高度標號，\(d_f(u)\) 的值既可能較 \(d_f(v)\) 增加 \(1\)，也可能保持不變或減少。因此，在整條路徑被添加完成後，我們得到 \(d_{f'}(s) \geq d_f(s)\)，其中取等的充要條件是 \(d_f(u) = d_f(v) + 1\) 對於 \((u, v) \in p\) 恆成立。如果該不等式不能取等，則有 \(d_{f'}(s) > d_f(s)\)——即我們想要的結論「層次圖的層數在增廣過程中嚴格單增」。以下我們嘗試證明該不等式不能取等。

考慮反證，我們假設 \(d_{f'}(s) = d_f(s)\) 成立，並嘗試導出矛盾。現在我們斷言，在 \(G'_L\) 上，\(p\) 至少包含一條邊 \((u, v)\) 滿足 \((u, v)\) 在 \(G_L\) 上不存在。如果沒有這樣的邊，考慮到 \(d_f(s) = d_{f'}(s)\)，結合層次圖和阻塞流的定義，\(G_L\) 上的增廣應尚未完成。為了不產生以上矛盾，我們的斷言只好是正確的。

令 \((u, v)\) 是滿足斷言條件的那條邊，其滿足斷言的原因只能是以下二者之一。

\((u, v) \in E_f\) 但 \(d_f(u) \leq d_f(v) + 1\) 未取等，故根據層次圖的定義可知 \((u, v) \not \in E_L\)，並在增廣後新一輪重分層中被加入到 \(E'_L\) 中；
\((u, v) \not \in E_f\)，這意味着 \((u, v)\) 這條邊的產生是當前輪次增廣中阻塞流經過 \((v, u)\) 並退流產生反向邊的結果，也即 \(d_f(u) = d_f(v) - 1\)。

由於我們無論以何種方式滿足斷言均得到 \(d_f(u) \neq d_f(v) + 1\)，也即 \(d_{f'}(s) \geq d_f(s)\) 取等的充要條件無法被滿足，這與反證假設 \(d_{f'}(s) = d_f(s)\) 衝突，原命題得證。

常見偽證另一則

考慮反證。假設層次圖的層數在一輪增廣結束後較原先相等，則層次圖上應仍存在至少一條從 \(s\) 到 \(t\) 的增廣路滿足相鄰兩點間的層數差為 \(1\)。這條增廣路未被增廣説明該輪增廣尚未結束。為了不產生上述矛盾，原命題成立。

Bug

「一輪增廣結束後新的層次圖上 \(s\)-\(t\) 最短路較原先相等」並不能推得「舊的層次圖上該輪增廣尚未結束」。這是因為，沒有理由表明兩張層次圖的邊集相同，新的層次圖上的 \(s\)-\(t\) 最短路有可能經過舊的層次圖上不存在的邊。

將單輪增廣的時間複雜度 \(O(|V||E|)\) 與增廣輪數 \(O(|V|)\) 相乘，Dinic 算法的時間複雜度是 \(O(|V|^2|E|)\)。

如果需要令 Dinic 算法的實際運行時間接近其理論上界，我們需要構造有特殊性質的網絡作為輸入。由於在算法競賽實踐中，對於網絡流知識相關的考察常側重於將原問題建模為網絡流問題的技巧。此時，我們的建模通常不包含令 Dinic 算法執行緩慢的特殊性質；恰恰相反，Dinic 算法在大部分圖上效率非常優秀。因此，網絡流問題的數據範圍通常較大，「將 \(|V|, |E|\) 的值代入 \(|V|^2|E|\) 以估計運行時間」這一方式並不適用。實際上，進行準確的估計需要選手對 Dinic 算法的實際效率有一定的經驗，讀者可以多加練習。

特殊情形下的時間複雜度分析

在一些性質良好的圖上，Dinic 算法有更好的時間複雜度。

對於網絡 \(G = (V, E)\)，如果其所有邊容量均為 \(1\)，即 \(c(u, v) \in \{0, 1\}\) 對於 \((u, v) \in E\) 恆成立，則我們稱 \(G\) 是單位容量（Unit Capacity）的。

在單位容量的網絡中，Dinic 算法的單輪增廣的時間複雜度為 \(O(|E|)\)。

證明

這是因為，每次增廣都會導致增廣路上的所有邊均飽和並消失，故單輪增廣中每條邊只能被增廣一次。

在單位容量的網絡中，Dinic 算法的增廣輪數是 \(O(|E|^{\frac{1}{2}})\) 的。

證明

以源點 \(s\) 為中心分層，記 \(d_f(u)\) 為 \(G_f\) 上結點 \(u\) 到源點 \(s\) 的距離。另外，我們定義將點集 \(\left\{u \mid u \in V, d_f(u) = k \right\}\) 定義為編號為 \(k\) 的層次 \(D_k\)，並記 \(S_k = \cup_{i \leq k} D_i\)。

假設我們已經進行了 \(|E|^{\frac{1}{2}}\) 輪增廣。根據鴿巢原理，至少存在一個 \(k\) 滿足邊集 \(\left\{ (u, v) \mid u \in D_k, v \in D_{k+1}, (u, v) \in E_f \right\}\) 的大小不超過 \(\frac {|E|} {|E|^{\frac{1}{2}}} \approx |E|^{\frac{1}{2}}\)。顯然，\(\{S_k, V - S_k\}\) 是 \(G_f\) 上的 \(s\)-\(t\) 割，且其割容量不超過 \(|E|^{\frac{1}{2}}\)。根據最大流最小割定理，\(G_f\) 上的最大流不超過 \(|E|^{\frac{1}{2}}\)，也即 \(G_f\) 上最多還能執行 \(|E|^{\frac{1}{2}}\) 輪增廣。因此，總增廣輪數是 \(O(|E|^{\frac{1}{2}})\) 的。

在單位容量的網絡中，Dinic 算法的增廣輪數是 \(O(|V|^{\frac{2}{3}})\) 的。

證明

假設我們已經進行了 \(2 |V|^{\frac{2}{3}}\) 輪增廣。由於至多有半數的（\(|V|^{\frac{2}{3}}\) 個）層次包含多於 \(|V|^{\frac{1}{3}}\) 個點，故無論我們如何分配所有層次的大小，至少存在一個 \(k\) 滿足相鄰兩個層次同時包含不多於 \(|V|^{\frac{1}{3}}\) 個點，即 \(|D_k| \leq |V|^{\frac{1}{3}}\) 且 \(|D_{k+1}| \leq |V|^{\frac{1}{3}}\)。

為最大化 \(D_k\) 和 \(D_{k+1}\) 之間的邊數，我們假定這是一個完全二分圖，此時邊集 \(\left\{ (u, v) \mid u \in D_k, v \in D_{k+1}, (u, v) \in E_f \right\}\) 的大小不超過 \(|V|^{\frac{2}{3}}\)。顯然，\(\{S_k, V - S_k\}\) 是 \(G_f\) 上的 \(s\)-\(t\) 割，且其割容量不超過 \(|V|^{\frac{2}{3}}\)。根據最大流最小割定理，\(G_f\) 上的最大流不超過 \(|V|^{\frac{2}{3}}\)，也即 \(G_f\) 上最多還能執行 \(|V|^{\frac{2}{3}}\) 輪增廣。因此，總增廣輪數是 \(O(|V|^{\frac{2}{3}})\) 的。

在單位容量的網絡中，如果除源匯點外每個結點 \(u\) 都滿足 \(\mathit{deg}_{\mathit{in}}(u) = 1\) 或 \(\mathit{deg}_{\mathit{out}}(u) = 1\)，則 Dinic 算法的增廣輪數是 \(O(|V|^{\frac{1}{2}})\) 的。其中，\(\mathit{deg}_{\mathit{in}}(u)\) 和 \(\mathit{deg}_{\mathit{out}}(u)\) 分別代表結點 \(u\) 的入度和出度。

證明

我們引入以下引理——對於這一形式的網絡，其上的任意流總是可以分解成若干條單位流量的、點不交 的增廣路。

假設我們已經進行了 \(|V|^{\frac{1}{2}}\) 輪增廣。根據層次圖的定義，此時任意新的增廣路的長度至少為 \(|V|^{\frac{1}{2}}\)。

考慮 \(G_f\) 上的最大流的增廣路分解，我們得到的增廣路的數量不能多於 \(\frac {|V|} {|V|^{\frac{1}{2}}} \approx |V|^{\frac{1}{2}}\)，這意味着 \(G_f\) 上最多還能執行 \(|V|^{\frac{1}{2}}\) 輪增廣。因此，總增廣輪數是 \(O(|V|^{\frac{1}{2}})\) 的。

綜上，我們得出一些推論。

在單位容量的網絡上，Dinic 算法的總時間複雜度是 \(O(|E| \min(|E|^\frac{1}{2}, |V|^{\frac{2}{3}}))\)。
在單位容量的網絡上，如果除源匯點外每個結點 \(u\) 都滿足 \(\mathit{deg}_{\mathit{in}}(u) = 1\) 或 \(\mathit{deg}_{\mathit{out}}(u) = 1\)，Dinic 算法的總時間複雜度是 \(O(|E||V|^{\frac{1}{2}})\)。對於二分圖最大匹配問題，我們常使用 Hopcroft–Karp 算法解決，而這一算法實際上是 Dinic 算法在滿足上述度數限制的單位容量網絡上的特例。

代碼實現

參考代碼

struct MF {
  struct edge {
    int v, nxt, cap, flow;
  } e[N];

  int fir[N], cnt = 0;

  int n, S, T;
  ll maxflow = 0;
  int dep[N], cur[N];

  void init() {
    memset(fir, -1, sizeof fir);
    cnt = 0;
  }

  void addedge(int u, int v, int w) {
    e[cnt] = {v, fir[u], w, 0};
    fir[u] = cnt++;
    e[cnt] = {u, fir[v], 0, 0};
    fir[v] = cnt++;
  }

  bool bfs() {
    queue<int> q;
    memset(dep, 0, sizeof(int) * (n + 1));

    dep[S] = 1;
    q.push(S);
    while (q.size()) {
      int u = q.front();
      q.pop();
      for (int i = fir[u]; ~i; i = e[i].nxt) {
        int v = e[i].v;
        if ((!dep[v]) && (e[i].cap > e[i].flow)) {
          dep[v] = dep[u] + 1;
          q.push(v);
        }
      }
    }
    return dep[T];
  }

  int dfs(int u, int flow) {
    if ((u == T) || (!flow)) return flow;

    int ret = 0;
    for (int& i = cur[u]; ~i; i = e[i].nxt) {
      int v = e[i].v, d;
      if ((dep[v] == dep[u] + 1) &&
          (d = dfs(v, min(flow - ret, e[i].cap - e[i].flow)))) {
        ret += d;
        e[i].flow += d;
        e[i ^ 1].flow -= d;
        if (ret == flow) return ret;
      }
    }
    return ret;
  }

  void dinic() {
    while (bfs()) {
      memcpy(cur, fir, sizeof(int) * (n + 1));
      maxflow += dfs(S, INF);
    }
  }
} mf;

MPM 算法

MPM(Malhotra, Pramodh-Kumar and Maheshwari) 算法得到最大流的方式有兩種：使用基於堆的優先隊列，時間複雜度為 \(O(n^3\log n)\)；常用 BFS 解法，時間複雜度為 \(O(n^3)\)。注意，本章節只專注於分析更優也更簡潔的 \(O(n^3)\) 算法。

MPM 算法的整體結構和 Dinic 算法類似，也是分階段運行的。在每個階段，在 \(G\) 的殘量網絡的分層網絡中找到增廣路。它與 Dinic 算法的主要區別在於尋找增廣路的方式不同：MPM 算法中尋找增廣路的部分的只花了 \(O(n^2)\), 時間複雜度要優於 Dinic 算法。

MPM 算法需要考慮頂點而不是邊的容量。在分層網絡 \(L\) 中，如果定義點 \(v\) 的容量 \(p(v)\) 為其傳入殘量和傳出殘量的最小值，則有：

\[ \begin{aligned} p_{in}(v) &= \sum\limits_{(u,v) \in L} (c(u, v) - f(u, v)) \\ p_{out}(v) &= \sum\limits_{(v,u) \in L} (c(v, u) - f(v, u)) \\ p(v) &= \min (p_{in}(v), p_{out}(v)) \end{aligned} \]

我們稱節點 \(r\) 是參考節點當且僅當 \(p(r) = \min {p(v)}\)。對於一個參考節點 \(r\)，我們一定可以讓經過 \(r\) 的流量增加 \(p(r)\) 以使其容量變為 \(0\)。這是因為 \(L\) 是有向無環圖且 \(L\) 中節點容量至少為 \(p(r)\)，所以我們一定能找到一條從 \(s\) 經過 \(r\) 到達 \(t\) 的有向路徑。那麼我們讓這條路上的邊流量都增加 \(p(r)\) 即可。這條路即為這一階段的增廣路。尋找增廣路可以用 BFS。增廣完之後所有滿流邊都可以從 \(L\) 中刪除，因為它們不會在此階段後被使用。同樣，所有與 \(s\) 和 \(t\) 不同且沒有出邊或入邊的節點都可以刪除。

時間複雜度分析

MPM 算法的每個階段都需要 \(O(V^2)\)，因為最多有 \(V\) 次迭代（因為至少刪除了所選的參考節點），並且在每次迭代中，我們刪除除最多 \(V\) 之外經過的所有邊。求和，我們得到 \(O(V^2+E)=O(V^2)\)。由於階段總數少於 \(V\)，因此 MPM 算法的總運行時間為 \(O(V^3)\)。

階段總數小於 V 的證明

MPM 算法在少於 \(V\) 個階段內結束。為了證明這一點，我們必須首先證明兩個引理。

引理 1：每次迭代後，從 \(s\) 到每個點的距離不會減少，也就是説，\(level_{i+1}[v] \ge level_{i}[v]\)。

證明：固定一個階段 \(i\) 和點 \(v\)。考慮 \(G_{i}^R\) 中從 \(s\) 到 \(v\) 的任意最短路徑 \(P\)。\(P\) 的長度等於 \(level_{i}[v]\)。注意 \(G_{i}^R\) 只能包含 \(G_{i}^R\) 的後向邊和前向邊。如果 \(P\) 沒有 \(G_{i}^R\) 的後邊，那麼 \(level_{i+1}[v] \ge level_{i}[v]\)。因為 \(P\) 也是 \(G_{i}^R\) 中的一條路徑。現在，假設 \(P\) 至少有一個後向邊且第一個這樣的邊是 \((u,w)\)，那麼 \(level_{i+1}[u] \ge level_{i}[u]\)（因為第一種情況）。邊 \((u,w)\) 不屬於 \(G_{i}^R\)，因此 \((u,w)\) 受到前一次迭代的增廣路的影響。這意味着 \(level_{i}[u] = level_{i}[w]+1\)。此外，\(level_{i+1}[w] = level_{i+1}[u]+1\)。從這兩個方程和 \(level_{i+1}[u] \ge level_{i}[u]\) 我們得到 \(level_{i+1}[w] \ge level_{i}[w]+2\)。路徑的剩餘部分也可以使用相同思想。

引理 2：\(level_{i+1}[t] > level_{i}[t]\)。

證明：從引理一我們得出，\(level_{i+1}[t] \ge level_{i}[t]\)。假設 \(level_{i+1}[t] = level_{i}[t]\)，注意 \(G_{i}^R\) 只能包含 \(G_{i}^R\) 的後向邊和前向邊。這意味着 \(G_{i}^R\) 中有一條最短路徑未被增廣路阻塞。這就形成了矛盾。

實現

參考代碼

struct MPM {
  struct FlowEdge {
    int v, u;
    long long cap, flow;

    FlowEdge() {}

    FlowEdge(int _v, int _u, long long _cap, long long _flow)
        : v(_v), u(_u), cap(_cap), flow(_flow) {}

    FlowEdge(int _v, int _u, long long _cap)
        : v(_v), u(_u), cap(_cap), flow(0ll) {}
  };

  const long long flow_inf = 1e18;
  vector<FlowEdge> edges;
  vector<char> alive;
  vector<long long> pin, pout;
  vector<list<int> > in, out;
  vector<vector<int> > adj;
  vector<long long> ex;
  int n, m = 0;
  int s, t;
  vector<int> level;
  vector<int> q;
  int qh, qt;

  void resize(int _n) {
    n = _n;
    ex.resize(n);
    q.resize(n);
    pin.resize(n);
    pout.resize(n);
    adj.resize(n);
    level.resize(n);
    in.resize(n);
    out.resize(n);
  }

  MPM() {}

  MPM(int _n, int _s, int _t) {
    resize(_n);
    s = _s;
    t = _t;
  }

  void add_edge(int v, int u, long long cap) {
    edges.push_back(FlowEdge(v, u, cap));
    edges.push_back(FlowEdge(u, v, 0));
    adj[v].push_back(m);
    adj[u].push_back(m + 1);
    m += 2;
  }

  bool bfs() {
    while (qh < qt) {
      int v = q[qh++];
      for (int id : adj[v]) {
        if (edges[id].cap - edges[id].flow < 1) continue;
        if (level[edges[id].u] != -1) continue;
        level[edges[id].u] = level[v] + 1;
        q[qt++] = edges[id].u;
      }
    }
    return level[t] != -1;
  }

  long long pot(int v) { return min(pin[v], pout[v]); }

  void remove_node(int v) {
    for (int i : in[v]) {
      int u = edges[i].v;
      auto it = find(out[u].begin(), out[u].end(), i);
      out[u].erase(it);
      pout[u] -= edges[i].cap - edges[i].flow;
    }
    for (int i : out[v]) {
      int u = edges[i].u;
      auto it = find(in[u].begin(), in[u].end(), i);
      in[u].erase(it);
      pin[u] -= edges[i].cap - edges[i].flow;
    }
  }

  void push(int from, int to, long long f, bool forw) {
    qh = qt = 0;
    ex.assign(n, 0);
    ex[from] = f;
    q[qt++] = from;
    while (qh < qt) {
      int v = q[qh++];
      if (v == to) break;
      long long must = ex[v];
      auto it = forw ? out[v].begin() : in[v].begin();
      while (true) {
        int u = forw ? edges[*it].u : edges[*it].v;
        long long pushed = min(must, edges[*it].cap - edges[*it].flow);
        if (pushed == 0) break;
        if (forw) {
          pout[v] -= pushed;
          pin[u] -= pushed;
        } else {
          pin[v] -= pushed;
          pout[u] -= pushed;
        }
        if (ex[u] == 0) q[qt++] = u;
        ex[u] += pushed;
        edges[*it].flow += pushed;
        edges[(*it) ^ 1].flow -= pushed;
        must -= pushed;
        if (edges[*it].cap - edges[*it].flow == 0) {
          auto jt = it;
          ++jt;
          if (forw) {
            in[u].erase(find(in[u].begin(), in[u].end(), *it));
            out[v].erase(it);
          } else {
            out[u].erase(find(out[u].begin(), out[u].end(), *it));
            in[v].erase(it);
          }
          it = jt;
        } else
          break;
        if (!must) break;
      }
    }
  }

  long long flow() {
    long long ans = 0;
    while (true) {
      pin.assign(n, 0);
      pout.assign(n, 0);
      level.assign(n, -1);
      alive.assign(n, true);
      level[s] = 0;
      qh = 0;
      qt = 1;
      q[0] = s;
      if (!bfs()) break;
      for (int i = 0; i < n; i++) {
        out[i].clear();
        in[i].clear();
      }
      for (int i = 0; i < m; i++) {
        if (edges[i].cap - edges[i].flow == 0) continue;
        int v = edges[i].v, u = edges[i].u;
        if (level[v] + 1 == level[u] && (level[u] < level[t] || u == t)) {
          in[u].push_back(i);
          out[v].push_back(i);
          pin[u] += edges[i].cap - edges[i].flow;
          pout[v] += edges[i].cap - edges[i].flow;
        }
      }
      pin[s] = pout[t] = flow_inf;
      while (true) {
        int v = -1;
        for (int i = 0; i < n; i++) {
          if (!alive[i]) continue;
          if (v == -1 || pot(i) < pot(v)) v = i;
        }
        if (v == -1) break;
        if (pot(v) == 0) {
          alive[v] = false;
          remove_node(v);
          continue;
        }
        long long f = pot(v);
        ans += f;
        push(v, s, f, false);
        push(v, t, f, true);
        alive[v] = false;
        remove_node(v);
      }
    }
    return ans;
  }
};

ISAP

在 Dinic 算法中，我們每次求完增廣路後都要跑 BFS 來分層，有沒有更高效的方法呢？

答案就是下面要介紹的 ISAP 算法。

過程

和 Dinic 算法一樣，我們還是先跑 BFS 對圖上的點進行分層，不過與 Dinic 略有不同的是，我們選擇在反圖上，從 \(t\) 點向 \(s\) 點進行 BFS。

執行完分層過程後，我們通過 DFS 來找增廣路。

增廣的過程和 Dinic 類似，我們只選擇比當前點層數少 \(1\) 的點來增廣。

與 Dinic 不同的是，我們並不會重跑 BFS 來對圖上的點重新分層，而是在增廣的過程中就完成重分層過程。

具體來説，設 \(i\) 號點的層為 \(d_i\)，當我們結束在 \(i\) 號點的增廣過程後，我們遍歷殘量網絡上 \(i\) 的所有出邊，找到層最小的出點 \(j\)，隨後令 \(d_i \gets d_j+1\)。特別地，若殘量網絡上 \(i\) 無出邊，則 \(d_i \gets n\)。

容易發現，當 \(d_s \geq n\) 時，圖上不存在增廣路，此時即可終止算法。

和 Dinic 類似，ISAP 中也存在 當前弧優化。

而 ISAP 還存在另外一個優化，我們記錄層數為 \(i\) 的點的數量 \(num_i\)，每當將一個點的層數從 \(x\) 更新到 \(y\) 時，同時更新 \(num\) 數組的值，若在更新後 \(num_x=0\)，則意味着圖上出現了斷層，無法再找到增廣路，此時可以直接終止算法（實現時直接將 \(d_s\) 標為 \(n\)），該優化被稱為 GAP 優化。

實現

參考代碼

struct Edge {
  int from, to, cap, flow;

  Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};

bool operator<(const Edge& a, const Edge& b) {
  return a.from < b.from || (a.from == b.from && a.to < b.to);
}

struct ISAP {
  int n, m, s, t;
  vector<Edge> edges;
  vector<int> G[maxn];
  bool vis[maxn];
  int d[maxn];
  int cur[maxn];
  int p[maxn];
  int num[maxn];

  void AddEdge(int from, int to, int cap) {
    edges.push_back(Edge(from, to, cap, 0));
    edges.push_back(Edge(to, from, 0, 0));
    m = edges.size();
    G[from].push_back(m - 2);
    G[to].push_back(m - 1);
  }

  bool BFS() {
    memset(vis, 0, sizeof(vis));
    queue<int> Q;
    Q.push(t);
    vis[t] = 1;
    d[t] = 0;
    while (!Q.empty()) {
      int x = Q.front();
      Q.pop();
      for (int i = 0; i < G[x].size(); i++) {
        Edge& e = edges[G[x][i] ^ 1];
        if (!vis[e.from] && e.cap > e.flow) {
          vis[e.from] = 1;
          d[e.from] = d[x] + 1;
          Q.push(e.from);
        }
      }
    }
    return vis[s];
  }

  void init(int n) {
    this->n = n;
    for (int i = 0; i < n; i++) G[i].clear();
    edges.clear();
  }

  int Augment() {
    int x = t, a = INF;
    while (x != s) {
      Edge& e = edges[p[x]];
      a = min(a, e.cap - e.flow);
      x = edges[p[x]].from;
    }
    x = t;
    while (x != s) {
      edges[p[x]].flow += a;
      edges[p[x] ^ 1].flow -= a;
      x = edges[p[x]].from;
    }
    return a;
  }

  int Maxflow(int s, int t) {
    this->s = s;
    this->t = t;
    int flow = 0;
    BFS();
    memset(num, 0, sizeof(num));
    for (int i = 0; i < n; i++) num[d[i]]++;
    int x = s;
    memset(cur, 0, sizeof(cur));
    while (d[s] < n) {
      if (x == t) {
        flow += Augment();
        x = s;
      }
      int ok = 0;
      for (int i = cur[x]; i < G[x].size(); i++) {
        Edge& e = edges[G[x][i]];
        if (e.cap > e.flow && d[x] == d[e.to] + 1) {
          ok = 1;
          p[e.to] = G[x][i];
          cur[x] = i;
          x = e.to;
          break;
        }
      }
      if (!ok) {
        int m = n - 1;
        for (int i = 0; i < G[x].size(); i++) {
          Edge& e = edges[G[x][i]];
          if (e.cap > e.flow) m = min(m, d[e.to]);
        }
        if (--num[d[x]] == 0) break;
        num[d[x] = m + 1]++;
        cur[x] = 0;
        if (x != s) x = edges[p[x]].from;
      }
    }
    return flow;
  }
};

Push-Relabel 預流推進算法

該方法在求解過程中忽略流守恆性，並每次對一個結點更新信息，以求解最大流。

通用的預流推進算法

首先我們介紹預流推進算法的主要思想，以及一個可行的暴力實現算法。

預流推進算法通過對單個結點的更新操作，直到沒有結點需要更新來求解最大流。

算法過程維護的流函數不一定保持流守恆性，對於一個結點，我們允許進入結點的流超過流出結點的流，超過的部分被稱為結點 \(u(u\in V-\{s,t\})\) 的 超額流 \(e(u)\)：

\[ e(u)=\sum_{(x,u)\in E}f(x,u)-\sum_{(u,y)\in E}f(u,y) \]

若 \(e(u)>0\)，稱結點 \(u\) 溢出⁶，注意當我們提到溢出結點時，並不包括 \(s\) 和 \(t\)。

預流推進算法維護每個結點的高度 \(h(u)\)，並且規定溢出的結點 \(u\) 如果要推送超額流，只能向高度小於 \(u\) 的結點推送；如果 \(u\) 沒有相鄰的高度小於 \(u\) 的結點，就修改 \(u\) 的高度（重貼標籤）。

高度函數⁷

準確地説，預流推進維護以下的一個映射 \(h:V\to \mathbf{N}\)：

\(h(s)=|V|,h(t)=0\)
\(\forall (u,v)\in E_f,h(u)\leq h(v)+1\)

稱 \(h\) 是殘量網絡 \(G_f=(V_f,E_f)\) 的高度函數。

引理 1：設 \(G_f\) 上的高度函數為 \(h\)，對於任意兩個結點 \(u,v\in V\)，如果 \(h(u)>h(v)+1\)，則 \((u,v)\) 不是 \(G_f\) 中的邊。

算法只會在 \(h(u)=h(v)+1\) 的邊執行推送。

推送（Push）

適用條件：結點 \(u\) 溢出，且存在結點 \(v((u,v)\in E_f,c(u,v)-f(u,v)>0,h(u)=h(v)+1)\)，則 push 操作適用於 \((u,v)\)。

於是，我們儘可能將超額流從 \(u\) 推送到 \(v\)，推送過程中我們只關心超額流和 \(c(u,v)-f(u,v)\) 的最小值，不關心 \(v\) 是否溢出。

如果 \((u,v)\) 在推送完之後滿流，將其從殘量網絡中刪除。

重貼標籤（Relabel）

適用條件：如果結點 \(u\) 溢出，且 \(\forall (u,v)\in E_f,h(u)\leq h(v)\)，則 relabel 操作適用於 \(u\)。

則將 \(h(u)\) 更新為 \(\min_{(u,v)\in E_f}h(v)+1\) 即可。

初始化

\[ \forall (u,v)\in E,~~f(u,v)=\begin{cases} c(u,v),&u=s\\ 0,&u\neq s \end{cases} \]

\[ \forall u\in V,~~h(u)=\begin{cases} |V|,&u=s\\ 0,&u\neq s \end{cases} \]

\[ e(u)=\sum_{(x,u)\in E}f(x,u)-\sum_{(u,y)\in E}f(u,y) \]

上述將 \((s,v)\in E\) 充滿流，並將 \(h(s)\) 抬高，使得 \((s,v)\notin E_f\)，因為 \(h(s)>h(v)\)，而且 \((s,v)\) 畢竟滿流，沒必要留在殘量網絡中；上述還將 \(e(s)\) 初始化為 \(\sum_{(s,v)\in E}f(s,v)\) 的相反數。

過程

我們每次掃描整個圖，只要存在結點 \(u\) 滿足 push 或 relabel 操作的條件，就執行對應的操作。

如圖，每個結點中間表示編號，左下表示高度值 \(h(u)\)，右下表示超額流 \(e(u)\)，結點顏色的深度也表示結點的高度；邊權表示 \(c(u,v)-f(u,v)\)，綠色的邊表示滿足 \(h(u)=h(v)+1\) 的邊 \((u,v)\)（即殘量網絡的邊 \(E_f\)）：

整個算法我們大致瀏覽一下過程，這裏筆者使用的是一個暴力算法，即暴力掃描是否有溢出的結點，有就更新

最後的結果

可以發現，最後的超額流一部分回到了 \(s\)，且除了源點匯點，其他結點都沒有溢出；這時的流函數 \(f\) 滿足流守恆性，為最大流，流量即為 \(e(t)\)。

但是實際上論文³指出只處理高度小於 \(n\) 的溢出節點也能獲得正確的最大流值，不過這樣一來算法結束的時候預流還不滿足流函數性質，不能知道每條邊上真實的流量。

實現

核心代碼

const int N = 1e4 + 4, M = 1e5 + 5, INF = 0x3f3f3f3f;
int n, m, s, t, maxflow, tot;
int ht[N], ex[N];

void init() {  // 初始化
  for (int i = h[s]; i; i = e[i].nex) {
    const int &v = e[i].t;
    ex[v] = e[i].v, ex[s] -= ex[v], e[i ^ 1].v = e[i].v, e[i].v = 0;
  }
  ht[s] = n;
}

bool push(int ed) {
  const int &u = e[ed ^ 1].t, &v = e[ed].t;
  int flow = min(ex[u], e[ed].v);
  ex[u] -= flow, ex[v] += flow, e[ed].v -= flow, e[ed ^ 1].v += flow;
  return ex[u];  // 如果 u 仍溢出，返回 1
}

void relabel(int u) {
  ht[u] = INF;
  for (int i = h[u]; i; i = e[i].nex)
    if (e[i].v) ht[u] = min(ht[u], ht[e[i].t]);
  ++ht[u];
}

HLPP 算法

最高標號預流推進算法（Highest Label Preflow Push）在上述通用的預流推送算法中，在每次選擇結點時，都優先選擇高度最高的溢出結點，其算法算法複雜度為 \(O(n^2\sqrt m)\)。

過程

具體地説，HLPP 算法過程如下：

初始化（基於預流推進算法）；
選擇溢出結點中高度最高的結點 \(u\)，並對它所有可以推送的邊進行推送；
如果 \(u\) 仍溢出，對它重貼標籤，回到步驟 2；
如果沒有溢出的結點，算法結束。

一篇對最大流算法實際表現進行測試的論文⁴表明，實際上基於預流的算法，有相當一部分時間都花在了重貼標籤這一步上。以下介紹兩種來自論文⁵的能顯著減少重貼標籤次數的優化。

BFS 優化

HLPP 的上界為 \(O(n^2\sqrt m)\)，但在使用時卡得比較緊；我們可以在初始化高度的時候進行優化。具體來説，我們初始化 \(h(u)\) 為 \(u\) 到 \(t\) 的最短距離；特別地，\(h(s)=n\)。

在 BFS 的同時我們順便檢查圖的連通性，排除無解的情況。

GAP 優化

HLPP 推送的條件是 \(h(u)=h(v)+1\)，而如果在算法的某一時刻，存在某個 \(k\)，使得 \(h(u)=k\) 的結點個數為 \(0\)，那麼對於 \(h(u)>k\) 的結點就永遠無法推送超額流到 \(t\)，因此只能送回 \(s\)，那麼我們就在這時直接讓他們的高度變成至少 \(n+1\)，以儘快推送回 \(s\)，減少重貼標籤的操作。

以下的實現採取論文⁴中的實現方法，使用 \(N*2-1\) 個桶 B，其中 B[i] 中記錄所有當前高度為 \(i\) 的溢出節點。加入了以上提到的兩種優化，並且只處理了高度小於 \(n\) 的溢出節點。

值得注意的是論文⁴中使用的桶是基於鏈表的棧，而 STL 中的 stack 默認的容器是 deque。經過簡單的測試發現 vector，deque，list 在本題的實際運行過程中效率區別不大。

實現

LuoguP4722【模板】最大流加強版/預流推進

#include <cstdio>
#include <cstring>
#include <queue>
#include <stack>
using namespace std;
const int N = 1200, M = 120000, INF = 0x3f3f3f3f;
int n, m, s, t;

struct qxx {
  int nex, t;
  long long v;
};

qxx e[M * 2 + 1];
int h[N + 1], cnt = 1;

void add_path(int f, int t, long long v) {
  e[++cnt] = qxx{h[f], t, v}, h[f] = cnt;
}

void add_flow(int f, int t, long long v) {
  add_path(f, t, v);
  add_path(t, f, 0);
}

int ht[N + 1];        // 高度;
long long ex[N + 1];  // 超額流;
int gap[N];           // gap 優化. gap[i] 為高度為 i 的節點的數量
stack<int> B[N];      // 桶 B[i] 中記錄所有 ht[v]==i 的v
int level = 0;        // 溢出節點的最高高度

int push(int u) {      // 儘可能通過能夠推送的邊推送超額流
  bool init = u == s;  // 是否在初始化
  for (int i = h[u]; i; i = e[i].nex) {
    const int &v = e[i].t;
    const long long &w = e[i].v;
    // 初始化時不考慮高度差為1
    if (!w || (init == false && ht[u] != ht[v] + 1) || t[v] == INF) continue;
    long long k = init ? w : min(w, ex[u]);
    // 取到剩餘容量和超額流的最小值，初始化時可以使源的溢出量為負數。
    if (v != s && v != t && !ex[v]) B[ht[v]].push(v), level = max(level, ht[v]);
    ex[u] -= k, ex[v] += k, e[i].v -= k, e[i ^ 1].v += k;  // push
    if (!ex[u]) return 0;  // 如果已經推送完就返回
  }
  return 1;
}

void relabel(int u) {  // 重貼標籤（高度）
  ht[u] = INF;
  for (int i = h[u]; i; i = e[i].nex)
    if (e[i].v) ht[u] = min(ht[u], ht[e[i].t]);
  if (++ht[u] < n) {  // 只處理高度小於 n 的節點
    B[ht[u]].push(u);
    level = max(level, ht[u]);
    ++gap[ht[u]];  // 新的高度，更新 gap
  }
}

bool bfs_init() {
  memset(ht, 0x3f, sizeof(ht));
  queue<int> q;
  q.push(t), ht[t] = 0;
  while (q.size()) {  // 反向 BFS, 遇到沒有訪問過的結點就入隊
    int u = q.front();
    q.pop();
    for (int i = h[u]; i; i = e[i].nex) {
      const int &v = e[i].t;
      if (e[i ^ 1].v && ht[v] > ht[u] + 1) ht[v] = ht[u] + 1, q.push(v);
    }
  }
  return ht[s] != INF;  // 如果圖不連通，返回 0
}

// 選出當前高度最大的節點之一, 如果已經沒有溢出節點返回 0
int select() {
  while (B[level].size() == 0 && level > -1) level--;
  return level == -1 ? 0 : B[level].top();
}

long long hlpp() {            // 返回最大流
  if (!bfs_init()) return 0;  // 圖不連通
  memset(gap, 0, sizeof(gap));
  for (int i = 1; i <= n; i++)
    if (ht[i] != INF) gap[ht[i]]++;  // 初始化 gap
  ht[s] = n;
  push(s);  // 初始化預流
  int u;
  while ((u = select())) {
    B[level].pop();
    if (push(u)) {  // 仍然溢出
      if (!--gap[ht[u]])
        for (int i = 1; i <= n; i++)
          if (i != s && ht[i] > ht[u] && ht[i] < n + 1)
            ht[i] = n + 1;  // 這裏重貼成 n+1 的節點都不是溢出節點
      relabel(u);
    }
  }
  return ex[t];
}

int main() {
  scanf("%d%d%d%d", &n, &m, &s, &t);
  for (int i = 1, u, v, w; i <= m; i++) {
    scanf("%d%d%d", &u, &v, &w);
    add_flow(u, v, w);
  }
  printf("%lld", hlpp());
  return 0;
}

感受一下運行過程

HLPP

其中 pic13 到 pic14 執行了 Relabel(4)，並進行了 GAP 優化。

腳註

http://pisces.ck.tp.edu.tw/~peng/index.php?action=showfile&file=f6cdf7ef750d7dc79c7d599b942acbaaee86a2e3e ↩
https://people.orie.cornell.edu/dpw/orie633/LectureNotes/lecture9.pdf ↩
Cherkassky B V, Goldberg A V. On implementing push-relabel method for the maximum flow problem[C]//International Conference on Integer Programming and Combinatorial Optimization. Springer, Berlin, Heidelberg, 1995: 157-171. ↩
Ahuja R K, Kodialam M, Mishra A K, et al. Computational investigations of maximum flow algorithms[J]. European Journal of Operational Research, 1997, 97(3): 509-542. ↩↩↩
Derigs U, Meier W. Implementing Goldberg's max-flow-algorithm—A computational investigation[J]. Zeitschrift für Operations Research, 1989, 33(6): 383-403. ↩
英語文獻中通常稱為「active」。 ↩
在英語文獻中，一個結點的高度通常被稱為「distance label」。此處使用的「高度」這個術語源自算法導論中的相關章節。你可以在機械工業出版社算法導論（原書第 3 版）的 P432 腳註中找到這麼做的理由。 ↩

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用

最大流

概述

Ford–Fulkerson 增廣

概述

最大流最小割定理

時間複雜度分析

Edmonds–Karp 算法

算法思想

時間複雜度分析

代碼實現

Dinic 算法

算法思想

時間複雜度分析

特殊情形下的時間複雜度分析

代碼實現

MPM 算法

時間複雜度分析

階段總數小於 V 的證明

實現

ISAP

過程

實現

Push-Relabel 預流推進算法

通用的預流推進算法

高度函數7

推送（Push）

重貼標籤（Relabel）

初始化

過程

實現

HLPP 算法

過程

BFS 優化

GAP 優化

實現

腳註

高度函數⁷