在微软苏州已经入职了一个月了,生活开始渐渐稳定下来。和在北京的各种不堪相比,我要感谢这次调职,得以让我从 easy 模式开始社畜生活。同时搁置好久的算法又可以得以继续研究下去了。

今天讲的算法是最大流最小割算法。实际上这个是我在《算法导论》上看来看去都看不明白的一个章节:术语太多,而且代码很少(我好菜啊)。幸运的是,hihoCoder 上面的一系列专题为我垫了不少砖头,让我得以掌握这一块知识。

对于一个最大流问题,其定义一般都是这样的:有一个有向图 G,源点为 s,汇点为 t,每条边都有其对应的流最大容量 c,源点 s 只有发出的流量,汇点 t 只有汇入的流量,求能从源点到汇点的最大可能流量。

那么,我先不想说证明过程,我直接给出结论:我们用一种反复查找的方法,试图在当前图 G 上找到一条从 s 到 t 的路径,其边上允许的流量非零。每次找到这一条路径之后,也就可以确定这条路径上的流量瓶颈,将路径上的边的可行流量减去这一流量瓶颈,在新图上进行下一次查找。我们一直持续这样的查找,直到无法从 s 到 t 走出一条流量瓶颈大于 0 的路径。这个算法叫做 Ford-Fulkerson 算法。这个查找方式可以使用 BFS 进行。

嘛,基本上就是这样。但是由于图中可能会存在反向边,对于我上面讲的这个情况,对于下面这种情况是不成立的:

An tricky example of flow network

如果我们用 BFS 找到了 A-B-D-F 这条可行路径,其瓶颈是 3,减去该瓶颈之后无法找到更新的路径。然而其最大流是 5:由 B 分出的流量有一份到 D,另外两份到 E,然后 A-C-D 有一个流量为 2 的路径。这时候,我们需要一个“残余网络”的概念帮助我们解决这个问题:每条边维护一条对应的反向边,其大小为当前在正向上面已经消耗掉的流量,而正向边的容量为当前正向剩余的可行流量。换句话说,正向边和反向边的存在,帮助我们维护了每条边还有多少剩余流量可用。如果有一条流量容量为 1 的边,其正向流量为 1000000,反向流量为 999999,这也是可行的!我们在残余网络上找到的一条可行路径,叫做增广路径。于是我上面描述的算法可以如下表达:

1
2
3
4
While ( findAugmentPath() ) // 判断是否有增广路
maxFlow = maxFlow + delta // 最大流增加
modifyGraph() // 对增广路进行修改
End While

与最小割的关系

我们通常说 max flow, minimun cut。那么什么被定义为 cut?在一个带有 s 和 t 的网络流里,有一种划分把点划分到两个不相交的 S 和 T 集合中,$s \in S, t \in T$. 净流 f(S, T) 被定义为穿过割 (S, T) 的流量之和(当一个割经过反向边的时候,反向边上的流量应为负值)。割的容量 C(S, T) 被定义为这条割上所有的容量之和(不包括反向边)。也就是说,$f(S, T) \le C(S, T)$。 可以证明,当前网络的流量总是等于任意一个割的净流。但是,任意的割不会有相同的容量。其中最小的那个割的容量对应的割,我们称之为最小割。

现在有一个结论:对于任意网络流图来说,最大流一定等于最小割的容量。这个结论是最大流最小割定理的直接推论,这个定理由三个等价的表达组成:

  1. f 是网络的最大流;
  2. 该网络的残余网络不包含任何增广路径。
  3. 流网络的某个切割 (S, T) 的容量等于 f。

因此,求出了最大流,也就求出了最小割的最大容量。使用上文所提及的 Ford-Fulkerson 算法,能够快速算出割、网络流量。

Production Code

使用 BFS 进行每次增广路径查找的 Ford-Fulkerson 实现叫做 Edmonds-Karp 算法。由于 BFS 时间复杂度(约)为 O(E),流量递增操作的总次数为 O(VE),该算法的时间复杂度为 $O(VE^2)$。这并不是最快的算法。现在有一种算法,通过记录每个点到汇点 t 的最短距离维持搜索顺序,使用 DFS 的方法进行增广路径的查找,能大大降低时间复杂度。这一算法的代码如下:

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
// 本代码对应于 http://hihocoder.com/problemset/problem/1369
#include <cstdio>
#include <cstring>
#include <queue>
const int inf = 0x3f3f3f3f;
const int maxn = 505, maxm = 40005;
int N, M;
typedef struct {
int u, v, c, prev;
}edge_t;
int head[maxn];
edge_t edges[maxm];
int source, sink;
int dep[maxn], depcnt[maxn], cur[maxn], stack[maxn];
void addedge(int u, int v, int c) {
static int id = 0;
edges[id].u = u; edges[id].v = v; edges[id].c = c; edges[id].prev = head[u]; head[u] = id++;
edges[id].u = v; edges[id].v = u; edges[id].c = 0; edges[id].prev = head[v]; head[v] = id++;
}
void bfs() {
memset(dep, 0xff, sizeof(dep));
memset(depcnt, 0, sizeof(depcnt));
dep[sink] = 0;
depcnt[dep[sink]] = 1;
std::queue<int> Q;
Q.push(sink);
int u, v;
while (!Q.empty()) {
u = Q.front();
Q.pop();
for (int e = head[u]; ~e; e = edges[e].prev) if (edges[e].c == 0 && dep[edges[e].v] == -1) {
v = edges[e].v;
Q.push(v);
dep[v] = dep[u] + 1;
depcnt[dep[v]]++;
}
}
}
int SAP() {
bfs();
memcpy(cur, head, sizeof(head));
int u = source, maxflow = 0, top = 0, i;
while (dep[source]<N) {
if (u == sink) {
//find the bottleneck
int delta = inf, p;
for (int t = top - 1; t >= 0; t--) {
int e = stack[t];
if (edges[e].c<delta) {
delta = edges[e].c;
p = t;
}
}
for (int t = top - 1; t >= 0; t--) {
int e = stack[t];
edges[e].c -= delta;
edges[e ^ 1].c += delta;
}
maxflow += delta;
top = p;
u = edges[stack[p]].u;
}
for (i = cur[u]; ~i; i = edges[i].prev) {
if (edges[i].c>0 && dep[edges[i].v] == dep[u] - 1) {
cur[u] = i;
u = edges[i].v;
stack[top++] = i;
break;
}
}
if (i == -1) {
depcnt[dep[u]]--;
if (depcnt[dep[u]] == 0) break;
int mindep = N;
for (int e = head[u]; ~e; e = edges[e].prev) {
if (edges[e].c>0 && mindep > dep[edges[e].v]) {
mindep = dep[edges[e].v];
cur[u] = e;
}
}
dep[u] = mindep + 1;
depcnt[dep[u]]++;
if (u != source) {
u = edges[stack[--top]].u;
}
}
}
return maxflow;
}
int main() {
freopen("testcase", "r", stdin);
scanf("%d%d", &N, &M);
source = 1, sink = N;
memset(head, 0xff, sizeof(head));
for (int i = 1; i <= M; i++) {
int u, v, c;
scanf("%d%d%d", &u, &v, &c);
addedge(u, v, c);
}
printf("%d\n", SAP());
return 0;
}

嘛,这个优化后的算法叫什么我也不知道,我只是从别人的解法里学习过来的。但是这个代码做了至少这么几个工作:

  1. 在进行 DFS 寻找增广路径之前用 BFS 从汇点搜索建立 dep 和 depcnt,即到汇点的最短距离,这样为在 DFS 的时候总是走最短路径提供了依据。
  2. 正向边和反向边成对存储,索引为偶数的总是正向边;正向边和反向边可以通过索引的异或切换;
  3. stack 里面存储通过 DFS 寻找的增广路径。
  4. 在每一次 DFS 寻找过后,回溯到瓶颈路径的起点重新搜索,节省时间;
  5. 维护一个 cur,其在物理意义上与 head 相同,也是为了直接从当前的可行边搜索,节省时间;
  6. 但是如果从某点找不到可行边,需要更新该点的 dep 与 depcnt。这也就意味着 cur 也要更新。

就是这样!但是有的文献指出,第一步建立 dep 与 depcnt 可以不进行。这主要的实现参考 http://hihocoder.com/contest/hiho117/solution/898175。

主要参考文献:

  1. hihocoder hiho一下 115、116、117、118、119 周
  2. 《算法导论》(第三版)第 26 章
  3. Max Flow-SAP-Improved Shortest Augmenting

留言