0%

匈牙利算法和KM算法


本来是在写最小费用最大流的,然后找到Going Home这道题。这道题的确是可以使用费用流来解,但是复杂度有点儿高, 然后想了一下,其实可以用匹配来解决。那么既然都用匹配了,这里就将算法进行一下记录,毕竟匹配算法还是比较常用的。

在这里,主要还是针对偶图来说,下面一步一步记录。


偶图

偶图也称为二部图,是指具有二分类(X,Y)的图,它的点集可以分解为两个(非空)子集X和Y,使得每条边的一个端点在X中, 另一个端点在Y中。

完全偶图就是X中每一个点与Y中每一个点都有连边。

下图就是一个偶图,但不是完全偶图:

偶图


匹配

给定图G=(V,E)。设M是E的一个不包含环的子集,它的任意两条边在G中均不相邻,则称M为G的匹配。 M中一条边的两个端点称为在M下是配对的,同时这两个点也是饱和的。若在M下,G的每一个顶点都是饱和的, 则称M为G的完美匹配。若G没有另外的匹配M’,使得|M|<|M’|,则称为G的最大匹配完美匹配也是最大匹配, 反之不成立。

对于偶图来说,一条边连接了X中的一个点和Y中的一个点,通常在这里就是为了找到最大匹配或者是完美匹配, 使得X中(或Y中)的点能够都达到饱和。

偶图都有最大匹配,但是只有在X的点数等于Y的点数的情况下,才可能有完美匹配。X与Y点数相等的完全偶图显然一定有完美匹配。

对于一个X与Y点数相等的完全偶图,如果每一条边都有一个权重,那么在其中寻找一个具有最大权和的完美匹配,就称这个匹配为最优匹配

最大匹配

上图是一个最大匹配,注意到它不是完美匹配(V8不饱和)。


匈牙利算法

匈牙利算法可以用于在偶图中寻找最大匹配

它的思路的确很简单,直接先举例:

匈牙利算法

要在上图中寻找最大匹配。

匈牙利算法

按照顺序来一个点一个点的进行匹配。

  • x1直接匹配到y1。
  • x2先查看y1,发现y1已经匹配了,所以又查看y2,它还没有匹配,所以x2匹配上y2。
  • x3先查看y2,发现y2已经匹配了,所以又查看y3,它还没有匹配,所以x3匹配上y3。

匈牙利算法

但是对于x4,很尴尬,它会发现y2,y3都已经匹配了。

这时就是体现匈牙利算法思想的时候了,那就是去寻求一条可扩路

可扩路就是起点与终点均为非饱和点的交错路交错路就是处于匹配中的边和未处于匹配中的边相互交错组成的路。

匈牙利算法

如上图所示,黄色边表示不在当前匹配中的边,红色边表示在匹配中的边,红黄交错,所以这里就是交错路。

对于上面一条路,它的起点和终点都是非饱和点,所有就是一条可扩路。

那么将黄色红色相互交换,就能得到下面那条路,相当于匹配边与未匹配边相互交换,这样这条路上所有的点就都饱和了,岂不美哉?

那么回到x4,我们希望以x4这个未饱和点为起点,去找到一条可扩路,然后匹配边与非匹配边一交换,x4就能饱和了。

匈牙利算法

这里使用深度优先搜索的方式来找可扩路,先找到y2,然后y2找到x2,然后x2再找到y4。于是就找到了可扩路(x4,y2,x2,y4)。

PS:事实上在这里深度优先搜索时,要进过几次搜索失败才能到找到上面这条可扩路:

  • x4,y2,x2,y1,x1,y2。(重复y2,搜索失败)
  • x4,y2,x2,y1,x1,y3,x3,y2(重复y2,搜索失败)
  • x4,y2,x2,y4。(搜索成功)

匈牙利算法

进行扩充,也就是将匹配(x2,y2)换成了(x2,y4)、(x4,y2)。

匈牙利算法

最后将(x5,y5)匹配上,这里就找到了最大匹配。

步骤总结:

  1. 初始化顶点编号,令n=1。
  2. 若n大于顶点数N,结束。否则步骤3。
  3. 深度优先搜索寻找顶点n的可扩路。如果找到,进行扩充,如果未找到,则这个顶点没法饱和。n = n + 1,返回步骤2.

KM算法

匈牙利算法用于在偶图中寻找最大匹配

那么如果要在完全偶图中寻找最优匹配,就需要用到这里的KM算法

KM算法的思想实际上并不复杂,举个简单的例如描述一下它的思想:

假如有M个工人去做M份工作,这些工作每个工人都会做,但是效率不一样,所以希望分配工作能够使得总效率最大。 先假设你不会KM算法,那么你就按照正常人的思路去想:

  1. 先让每个人都挑他们效率最高的那份工作,如果都不冲突,那显然分配已经最优了。否则有冲突,那么继续想办法。
  2. 对于冲突了的几个人,让他们中最叼的那个哥们换另一份工作试试,看看能不能解决冲突,这样损失的总效率最小。 比如A、B、C之间冲突了,假如一共5份工作(a,b,c,d,e), 他们的效率分别是A(100,100,99,99,99)B(8,7,4,3,2)C(5,5,1,1,1), 在第一轮里面,按照最大效率的选择,A想选ab中的一个,B想选a,C想选ab中的一个,显然他们三个冲突了。 假如这时d没人选,你如果拍脑袋想,“B你去做a,A你去做b,C你去做d,因为C你最垃圾,A做b的效率可是100”, 那么三人的总效率加起来就是8 + 100 + 1 = 109,但是很遗憾,不是最优。 这时你的分配应该是“B你去做a,C你去做b,A你去做d,因为C你就只能做这玩意,换一个别的你做不了,A换成d也差不多效率”, 这是三人的总效率加起来就是8 + 5 + 99 = 112,这时就是最好的分配方案。

当然这个例子比较简单,在复杂情况下冲突会发生多次。

那么按照完全偶图(X,Y)来说,思路就是:

  1. 控制目前图上能选边的数量,第一轮每个X都只有权重最大边能选。
  2. 使用匈牙利算法来分配,如果分配成功,这就是最优匹配,否则转步骤3。
  3. 对于冲突的部分,想办法添加边进来解决冲突,添加的边是冲突的几个点中权重损失最小的那一条,然后转步骤2。

上面的都是语言上的描述,很不准确,下面搬一下书上的数学语言描述。


可行顶点标号

定义:若在顶点集$X \bigcup Y$上的实值函数L适合下述条件:对所有的$x \in X$及$y \in Y$均有

则把这个函数称为该偶图的一个可行顶点标号。

可行顶点标号也就是每条边的权重都没有它的两端点标号之和大,不管偶图是什么样子,都存在一个可行顶点标号:

其实这就是上面的“每个X都只有权重最大边能选”。


相等子图

定义:若L是可行顶点标号,则$E_L$表示可行性顶点标号定义中,使得等式成立的那些边,即:

则称具有边集$E_L$的G的生成子图为对应于可行顶点标号L的相等子图,记为$G_L$。


相关定理

设L是G=(V,E)的可行顶点标号。若$G_L$包含完美匹配$M^*$,则$M^*$是G的最优匹配。

证明:

假设$G_L$包含完美匹配$M^*$。由于$G_L$是G的生成子图,所以$M^*$也就是G的完美匹配。于是

这是因为每个$e \in M^*$都属于这个相等子图,且$M^*$的边的端点覆盖V的每个顶点恰好一次。 另一方面,若M是G的任一完美匹配,则有

从上面两个式子可以看出,$W(M^*) \geqslant W(M)$。于是$M^*$是最优匹配。


Kuhn-Munkres

KM算法由Kuhn(1955)和Munkers(1957)提出,所以取名叫KM

它的思想就是:

  1. 给定一个可行顶点标号。
  2. 在相等子图上使用匈牙利算法,如果完美匹配,则也是最优匹配(上面已经证明)。否则进入下一步。
  3. 修改可行顶点标号,使得相等子图中冲突的几个点的边更多,再回到步骤2。

它的步骤如下(主要需要注意的就是它修改可行顶点标号的方式):

从任一可行顶点标号L开始,然后决定$G_L$,并且在$G_L$中选取任一匹配M。

  1. 若X是M饱和的,则M是完美匹配,且由定理15知M是最优匹配,在这种情况下,算法终止。 否则,令u是一个M非饱和点,置$S=\{u\}$,$T= \varnothing$。

  2. 若$N_{G_L}(S) \supset T$,则转到步骤3。否则$N_{G_L}(S)=T$。计算

    且由

    给出可行顶点标号$\hat{L}$(注意$\alpha_L > 0$且$N_{G_{\hat{L}}} \supseteq T$)。 以$\hat{L}$代替L,以$G_{\hat{L}}$代替$G_L$。

  3. 在$N_{G_L}(S) \setminus T$中选择一个顶点y。和上节中树的生长程序一样,考察y是否M饱和。 若y是M饱和的,且$yz \in M$,则用$S \cup \{z\}$代替S,用$T \cup \{y\}$代替T,再转到步骤2. 否则,设P是$G_L$中的M可扩(u,y)路,用${M}’=M \triangle E(P)$代替M,并转到步骤1。

上面的$N_{G_L}(S)$表示在$G_L$这个图中,S这个点集的所有相邻的点。
$N_{G_L}(S) \setminus T$表示$N_{G_L}(S)$这个集合减去T这个集合。
$M \triangle E(P)$在这里表示$M - M \cap E(p) + (E(p) - M \cap E(p))$,其实也就是在可扩路上进行匹配扩展,边集交换。


例子实际操作

只看上面的步骤其实还不太会实际操作,这里给一个实际的例子。

KM算法

上图就是一个具体的例子,左边就是权重矩阵,行表示x点,列表示y点, x顶点的标号已经初始化(按照上面描述的可行顶点标号中的方法初始化)。 图的右边就是现在的相等子图(边的权重等于两端点的标号之和)。

那么下面的步骤就是在相等子图上使用匈牙利算法

KM算法

很简单的,前三个x点就能找到匹配。

KM算法

但是到x4时,发现找不到匹配,即使经过搜索也找不到(肉眼就能看出来x1,x3,x4三个点,一起抢y2,y3两个点,显然不能完美匹配)。

那么修改顶点标号。

按照上面的KM算法步骤2中的公式$\alpha_L = \min_{x \in S \ y \notin T} \{L(x) + L(y) - w(xy)\}$, 这里的S就是{x1,x3,x4}这三个点, T就时{y2,y3}这两个点,可以算出$\alpha = 1$。

同样按照KM算法步骤2中的公式,用算出来的$\alpha = 1$,修改顶点标号,x1,x3,x4三个点标号减去$\alpha$,y2,y3两个点标号加上$\alpha$,

KM算法

图中左边的顶点标号已经修改,右边相等子图发生了变化,蓝色边就是新加入的边。

另外需要注意的是(x2,y2)这条边不在相等子图里面了,这个细节表明了相等子图中的边可不是越来越多!增加边的同时也会有减少的边! 所以上面算法说了,可以从任一可行标点开始!当然,通常写代码时从X标号值最大开始,思考起来会比较连贯。

KM算法

然后很简单的就能得到完美匹配,也就是这里的最优匹配。


KM算法复杂度

注意到每一次找可扩路失败会修改顶点标号,在修改顶点标号后,都会至少有一个新的y端点加入, 也就是对于一个未饱和点,最多修改$O(n)$次顶点标号就能找到可扩路。

  1. KM算法需要找到可扩路的数量为$O(n)$。
  2. 每次尝试寻找可扩路的复杂度$O(n^2)$。
  3. 尝试找到可扩路,最多尝试$O(n)$次,也就是最多需要修改$O(n)$次顶点标号。
  4. 每一次修改顶点标号需要$O(n^2)$的复杂度。

从两个方面看,

  1. 可扩路的数量 * 尝试寻找可扩路次数 * 尝试寻找可扩路的复杂度为$O(n^4)$。
  2. 可扩路的数量 * 最多需要修改顶点标号的次数 * 每一次修改顶点标号的复杂度为$O(n^4)$。

所以KM算法的复杂度为上面两个部分相加,所有最后还是$O(n^4)$。

PS:网上很多博客都说可以把修改顶点标号的复杂度降到$O(n^3)$(的确可以),然后KM算法的复杂度就能降到$O(n^3)$??? 但是他们用的都是DFS,对此我表示深刻的怀疑,感觉应该是把复杂度降到$O(n^4) + O(n^3)$,但是实际上应该还是$O(n^4)$的复杂度, 不过降低了常数因子,应该是快了一倍,但是并不能有数量级上的变化。
只有使用BFS的写法才能达到$O(n^3)$,后面进行解释。


KM算法个人理解

可以从上面例子看出KM算法的精髓,那就是在相等子图找不到完美匹配时,修改顶点标号, 让冲突的那几个点在相等子图里面能够有更多的边,并至少多连接一个y顶点。

这里之所以$\alpha$取min,是为了顶点标号之和损失最小。因为:

  1. 最优匹配的值就等于所有点顶点标号之和。
  2. 在最优匹配里面每一个顶点出现一次仅一次,相等子图的每条边权重等于它两端点标号和。

这里对于X标号值从大到小的算法步骤而言,每一次发生冲突,都是因为几个x点它们共同可以选择的y点的数量比它们小(例如上面例子中三个x点竞争两个y点), 也就是KM算法步骤2中的S集合的规模一定大于T集合,$|S|>|T|$。所以在进行标号修改时,由于S中的点减去标号值,T中的点加上标号值, 所以所有点的标号之和就一定会比之前小

  • 如果发生冲突,S集合的规模一定大于T集合,$|S|>|T|$。
  • 标号修改量$\alpha$一定是个大于零的值。
  • 在进行标号修改时,S中的每个点减去标号值,T中的每个点加上标号值。
  • 所有点的标号之和一定减小。

上面这段是为了表达KM算法有一种这个意味:“最开始期望一个最大的最优匹配,如果不行再一点一点降低期望。”


代码

是时候上代码了。


DFS版本

DFS版本(复杂度$O(n^4)$):

虽然我没有参考一个网络上的模板,但是写出来也大同小异…

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
public class KM_DFS {

private int[][] table = null; // 权重矩阵(方阵)
private int[] xl = null; // X标号值
private int[] yl = null; // Y标号值
private int[] xMatch = null; // X点对应的匹配点
private int[] yMatch = null; // Y点对应的匹配点
private int n = 0; // 矩阵维度
private int a = 0; // 标号修改量

public int solve(int[][] table) { // 入口,输入权重矩阵
this.table = table;
init();

for ( int x = 0; x < n; x++ ) { // 为每一个x寻找匹配
boolean[] S = new boolean[n], T = new boolean[n]; // S集合,T集合
a = Integer.MAX_VALUE;
while ( !dfs(S, T, x) ) { // 找到可扩路结束,否则修改标号值
LModified(S, T);
Arrays.fill(S, false);
Arrays.fill(T, false);
a = Integer.MAX_VALUE;
}
}

int value = 0;
for ( int x = 0; x < n; x++ ) {
value += table[x][xMatch[x]];
}
return value;
}

private boolean dfs(boolean[] S, boolean[] T, int x) { // 深度优先搜索
S[x] = true;
for ( int y = 0; y < n; y++ ) {
if ( T[y] ) {
continue;
}
int tmp = xl[x] + yl[y] - table[x][y];
if ( tmp == 0 ) { // 在相等子树中
T[y] = true;
if ( yMatch[y] == -1 || dfs(S, T, yMatch[y]) ) { // 1. y顶点没有匹配,那么进行匹配
xMatch[x] = y; // 2. dfs寻找可扩路成功,那么这条x,y就会因为可扩路的扩展而交换到匹配中
yMatch[y] = x;
return true;
}
} else { // 不在相等子树中
a = Math.min(tmp, a);
}
}
return false;
}

private void init() {
this.n = table.length;
this.xl = new int[n];
this.yl = new int[n];
Arrays.fill(xl, Integer.MIN_VALUE);
for ( int x = 0; x < n; x++ ) {
for ( int y = 0; y < n; y++ ) {
if ( table[x][y] > xl[x] ) {
xl[x] = table[x][y];
}
}
}
this.xMatch = new int[n];
this.yMatch = new int[n];
Arrays.fill(xMatch, -1);
Arrays.fill(yMatch, -1);
}

private void LModified(boolean[] S, boolean[] T) { // 修改标号值
for ( int i = 0; i < n; i++ ) {
if ( S[i] ) {
xl[i] -= a;
}
if ( T[i] ) {
yl[i] += a;
}
}
}

}

BFS版本

版本1(复杂度$O(n^4)$):

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
public class KM_BFS_1 {
private int[][] table = null; // 权重矩阵(方阵)
private int[] xl = null; // X标号值
private int[] yl = null; // Y标号值
private int[] xMatch = null; // X点对应的匹配点
private int[] yMatch = null; // Y点对应的匹配点
private int n = 0; // 矩阵维度

public int solve(int[][] table) { // 入口,输入权重矩阵
this.table = table;
init();

for ( int x = 0; x < n; x++ ) {
bfs(x);
}

int value = 0;
for ( int x = 0; x < n; x++ ) {
value += table[x][xMatch[x]];
}
return value;
}

private void bfs(int startX) { // 为一个x点寻找匹配
boolean find = false;
int endY = -1;
int[] yPre = new int[n]; // 标识搜索路径上y点的前一个点
boolean[] S = new boolean[n], T = new boolean[n]; // S集合,T集合
Arrays.fill(yPre, -1);

int a = Integer.MAX_VALUE;
int[] queue = new int[n]; // 队列
int qs = 0, qe = 0; // 队列开始结束索引
queue[qe++] = startX;
while (true) { // 循环直到找到匹配
while (qs < qe && !find) { // 队列不为空
int x = queue[qs++];
S[x] = true;
for (int y = 0; y < n; y++) {
int tmp = xl[x] + yl[y] - table[x][y];
if ( tmp == 0 ) { // 相等子树中的边
if (T[y]) {
continue;
}
T[y] = true;
yPre[y] = x;
if (yMatch[y] == -1) {
endY = y;
find = true;
break;
} else {
queue[qe++] = yMatch[y];
}
} else { // 不在相等子树中的边,记录一下最小差值
a = Math.min(a, tmp);
}
}
}
if ( find ) {
break;
}
qs = qe = 0;
for ( int i = 0; i < n; i++ ) { // 根据a修改标号值
if ( S[i] ) {
xl[i] -= a;
queue[qe++] = i; // 把所有在S中的点加回到队列中
}
if ( T[i] ) {
yl[i] += a;
}
}
a = Integer.MAX_VALUE;
}

while ( endY != -1 ) { // 找到可扩路最后的y点后,回溯并扩充
int preX = yPre[endY], preY = xMatch[preX];
xMatch[preX] = endY;
yMatch[endY] = preX;
endY = preY;
}
}

private void init() {
this.n = table.length;
this.xl = new int[n];
this.yl = new int[n];
Arrays.fill(xl, Integer.MIN_VALUE);
for ( int x = 0; x < n; x++ ) {
for ( int y = 0; y < n; y++ ) {
if ( table[x][y] > xl[x] ) {
xl[x] = table[x][y];
}
}
}
this.xMatch = new int[n];
this.yMatch = new int[n];
Arrays.fill(xMatch, -1);
Arrays.fill(yMatch, -1);
}
}

版本2(复杂度$O(n^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
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
public class KM_BFS_2 {
private int[][] table = null; // 权重矩阵(方阵)
private int[] xl = null; // X标号值
private int[] yl = null; // Y标号值
private int[] xMatch = null; // X点对应的匹配点
private int[] yMatch = null; // Y点对应的匹配点
private int n = 0; // 矩阵维度

public int solve(int[][] table) { // 入口,输入权重矩阵
this.table = table;
init();

for ( int x = 0; x < n; x++ ) {
bfs(x);
}

int value = 0;
for ( int x = 0; x < n; x++ ) {
value += table[x][xMatch[x]];
}
return value;
}

private void bfs(int startX) { // 为一个x点寻找匹配

boolean find = false;
int endY = -1;
int[] yPre = new int[n]; // 标识搜索路径上y点的前一个点
boolean[] S = new boolean[n], T = new boolean[n]; // S集合,T集合
int[] slackY = new int[n]; // Y点的松弛变量
Arrays.fill(yPre, -1);
Arrays.fill(slackY, Integer.MAX_VALUE);

int[] queue = new int[n]; // 队列
int qs = 0, qe = 0; // 队列开始结束索引
queue[qe++] = startX;
while (!find) { // 循环直到找到匹配
while (qs < qe && !find) { // 队列不为空
int x = queue[qs++];
S[x] = true;
for ( int y = 0; y < n; y++ ) {
if ( T[y] ) {
continue;
}
int tmp = xl[x] + yl[y] - table[x][y];
if ( tmp == 0 ) { // 相等子树中的边
T[y] = true;
yPre[y] = x;
if ( yMatch[y] == -1 ) {
endY = y;
find = true;
break;
} else {
queue[qe++] = yMatch[y];
}
} else if ( slackY[y] > tmp ) { // 不在相等子树中的边,看是否能够更新松弛变量
slackY[y] = tmp;
yPre[y] = x;
}
}
}
if ( find ) {
break;
}
int a = Integer.MAX_VALUE;
for ( int y = 0; y < n; y++ ) { // 找到最小的松弛值
if ( !T[y] ) {
a = Math.min(a, slackY[y]);
}
}
for ( int i = 0; i < n; i++ ) { // 根据a修改标号值
if ( S[i] ) {
xl[i] -= a;
}
if ( T[i] ) {
yl[i] += a;
}
}

qs = qe = 0;
for ( int y = 0; y < n; y++ ) { // 重要!!!控制修改标号之后需要检查的x点
if ( !T[y] && slackY[y] == a ) { // 查看那些y点新加入到T集合,注意,这些y点的前向x点都记录在了yPre里面,所以这些x点不用再次入队
T[y] = true;
if ( yMatch[y] == -1 ) { // 新加入的y点没有匹配,那么就找到可扩路了
endY = y;
find = true;
break;
} else { // 新加入的y点已经有匹配了,将它匹配的x加到队列
queue[qe++] = yMatch[y];
}
}
slackY[y] -= a; // 所有松弛值减去a。(对于T集合中的松弛值已经没用了,对于不在T集合里面的y点,
} // 它们的松弛值是通过S集合中的x点求出的,S集合中的x点的标号值在上面都减去了a,所以这里松弛值也要减去a)
}

while ( endY != -1 ) { // 找到可扩路最后的y点后,回溯并扩充
int preX = yPre[endY], preY = xMatch[preX];
xMatch[preX] = endY;
yMatch[endY] = preX;
endY = preY;
}
}

private void init() {
this.n = table.length;
this.xl = new int[n];
this.yl = new int[n];
Arrays.fill(xl, Integer.MIN_VALUE);
for ( int x = 0; x < n; x++ ) {
for ( int y = 0; y < n; y++ ) {
if ( table[x][y] > xl[x] ) {
xl[x] = table[x][y];
}
}
}
this.xMatch = new int[n];
this.yMatch = new int[n];
Arrays.fill(xMatch, -1);
Arrays.fill(yMatch, -1);
}
}

代码解释

不把一个算法写成代码,真的不能算看懂了这个算法。

写这几份代码花费了一天的时间,大部分时间都用在了BFS代码上,巧妙的地方太多了,先说一下BFS算法的思想:

BFS算法思想:BFS_1代码的思想其实没什么特殊的,只是和DFS的遍历方式不一样,所有减小了一些边的重复遍历次数。 但是对于BFS_2这份代码来说,第一个就是松弛变量,它使得每个x节点不用重复入队,第二个是数组实现队列,省掉了很多的容器操作。

松弛变量(slack数组)这个东西,在DFS里面没有很大的作用,但是在BFS里面就变得很重要,它是BFS_2算法的核心。

松弛变量:

  1. 每一次修改标号值,都在相等子图里面会增加y点。
  2. 每一次增加的y点,就是松弛变量值最小的那些y点。
  3. 加入的y点在可扩路中的前向点是在计算最小松弛变量$L(x) + L(y) - w(xy)$中的那个x点。
  4. 修改标号值时使用的最小松弛变量值是从那些不在T集合中的y点的松弛变量值中选出来的。
  5. 所有y点的松弛变量值都是由在S集合中的那些x点算出来的。
  6. 未被加入T集合的y点,在修改标号值时,由于S集合的所有x点的标号值减去了$\alpha$,所以这些y点的松弛变量值也要减去$\alpha$。

BFS中的queue:本来应该是使用集合Queue来操作,但是那样增加了许多复杂度,实际上每个x最多入队一次,所以这里可以使用数组来进行操作。


BenchMark

对上面三份代码进行测试,查看他们实际上的性能。

例子数 矩阵最大维度 权重范围 DFS用时(秒) BFS_1用时(秒) BFS_2用时(秒)
1000 300 10 2.013 0.34 0.389
1000 300 100 1.418 0.437 0.779
1000 300 1000 2.734 0.923 1.202
1000 300 10000 10.455 3.413 1.435
1000 300 100000 29.347 11.564 1.542
50 1000 10 3.651 0.231 0.263
50 1000 100 2.213 0.296 0.365
50 1000 1000 1.357 0.445 0.716
50 1000 10000 3.677 1.102 0.893
50 1000 100000 12.689 5.531 1.178

BFS全面优于DFS,BFS_2受到权重范围的影响较小。

上面的结果和之前的算法复杂度很接近,果然DFS是$O(n^4)$,BFS_1其实也是$O(n^4)$,BFS_2是真正的$O(n^3)$,当然在权重范围不大时,BFS_1是效果最好的。


例题

下面使用一些例题来参考。


1. Going Home

这道题来自杭州电子科技大学的OJ,题号1533,Going Home

Going Home

题目的意思就是,有几个人要回到房子里面去,每一个房子只能住一个人,现在人数和房子数一样,但是要怎么分配房子才能使的大家回房子的总路程最小。

我知道的解法有两种,第一种是将这看作图,然后使用最小费用最大流,这个方法可以AC,但是感觉太过于复杂了。

第二种方法就是把房子和人看成两个集合,相互连线,权重就是人到房子的距离,这样就转化为了一个寻找最优匹配的问题,不过这里是找路径总长最小, 那么为了沿用模板,把所有权重取一个负号就行。

代码就是用的上面的模板,这里就不贴了。

用时:

Use Time


2. 奔小康赚大钱

这道题来自HDU-2255

奔小康赚大钱

这道题直接就是一个求最优匹配的问题。

这道题对java不友好!!!,我上面的KM算法都试了一下,java就是过不了,艹了。看了提交统计,就没有java版本通过。

改成c++之后就通过了???

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
#include<iostream>
#include<cstdio>
#include<cstring>

using namespace std;

const int N = 310;
const int INF = 0x3f3f3f3f;

int table[N][N]; // 权重矩阵(方阵)
int xl[N], yl[N]; // 标号值
int xMatch[N], yMatch[N]; // 点对应的匹配点
bool S[N], T[N];
int yPre[N], slackY[N];
int queue[N]; // 队列
int n; // 矩阵维度


void bfs(int startX) { // 为一个x点寻找匹配
memset(yPre, -1, sizeof(yPre));
memset(S, 0, sizeof(S));
memset(T, 0, sizeof(T));
memset(T, 0, sizeof(T));
memset(slackY, INF, sizeof(slackY));

bool find = false;
int endY = -1;
int qs = 0, qe = 0; // 队列开始结束索引
queue[qe++] = startX;
while (!find) { // 循环直到找到匹配
while (qs < qe && !find) { // 队列不为空
int x = queue[qs++];
S[x] = true;
for (int y = 0; y < n; y++) {
if (T[y]) {
continue;
}
int tmp = xl[x] + yl[y] - table[x][y];
if (tmp == 0) { // 相等子树中的边
T[y] = true;
yPre[y] = x;
if (yMatch[y] == -1) {
endY = y;
find = true;
break;
}
else {
queue[qe++] = yMatch[y];
}
}
else if (slackY[y] > tmp) { // 不在相等子树中的边,看是否能够更新松弛变量
slackY[y] = tmp;
yPre[y] = x;
}
}
}
if (find) {
break;
}
int a = INF;
for (int y = 0; y < n; y++) { // 找到最小的松弛值
if (!T[y] && slackY[y] < a) {
a = slackY[y];
}
}
for (int i = 0; i < n; i++) { // 根据a修改标号值
if (S[i]) {
xl[i] -= a;
}
if (T[i]) {
yl[i] += a;
}
}

qs = qe = 0;
for (int y = 0; y < n; y++) {
if (!T[y] && slackY[y] == a) { // 查看那些y点新加入到T集合,注意,这些y点的前向x点都记录在了yPre里面,所以这些x点不用再次入队
T[y] = true;
if (yMatch[y] == -1) { // 新加入的y点没有匹配,那么就找到可扩路了
endY = y;
find = true;
break;
}
else { // 新加入的y点已经有匹配了,将它匹配的x加到队列
queue[qe++] = yMatch[y];
}
}
slackY[y] -= a; // 所有松弛值减去a。(对于T集合中的松弛值已经没用了,对于不在T集合里面的y点,
} // 它们的松弛值是通过S集合中的x点求出的,S集合中的x点的标号值在上面都减去了a,所以这里松弛值也要减去a)
}

while (endY != -1) { // 找到可扩路最后的y点后,回溯并扩充
int preX = yPre[endY], preY = xMatch[preX];
xMatch[preX] = endY;
yMatch[endY] = preX;
endY = preY;
}
}

int solve() {
memset(xMatch, -1, sizeof(xMatch));
memset(yMatch, -1, sizeof(yMatch));
memset(yl, 0, sizeof(yl));
for (int x = 0; x < n; x++) {
for (int y = 0; y < n; y++) {
if (table[x][y] > xl[x]) {
xl[x] = table[x][y];
}
}
}

for (int x = 0; x < n; x++) {
bfs(x);
}

int value = 0;
for (int x = 0; x < n; x++) {
value += table[x][xMatch[x]];
}
return value;
}

int main() {
while (~scanf("%d", &n)) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
scanf("%d", &table[i][j]);
int ans = solve();
printf("%d\n", ans);
}
return 0;
}

用时:

奔小康赚大钱


3. Cyclic Tour

这道题来自HDU-1853

Cyclic Tour

每一个点需要经过一次仅一次,所有的点都是单向边,所以:

  1. 每个点出入度都至少要大于1,否则不能形成环路。
  2. 每个点的出入度都大于1,那么一定有环路。
  3. 在最后的环游方案中,每个点被入度一次仅一次,出度一次仅一次。

那么把所有点既放到X集合里面,又放到Y集合里面,然后这个偶图的最优匹配就是这个题目的解。

代码从上面那份C++代码改一下就完成了。

Cyclic Tour


总结

  1. KM算法步骤比较简单,思想没那么简单,想写出一个好的代码更不简单。
  2. HDU对java不友好。

参考资料:

《图论及其应用》-高等教育出版社-张先迪、李正良