算法的乐趣——7.4　二部图与二分匹配

之前讨论稳定匹配问题的时候，我们把完美匹配定义为每个男人和女人都属于匹配中的某个对，并不是很直观，现在我们准备用图的术语更一般地表达完美匹配的概念。首先介绍一下二部图，二部图G=(V,E)是这样的一个图，它的顶点集合V可以划分为X和Y两个集合，它的边集合E中的每条边都有一个端点在X集合，另一个端点在Y集合。图7-2就是一个二部图。 现在给出针对二部图的匹配的定义，给定一个二部图G=(V,E)的子图M，如果M的边集中任意两条边都不依附于同一个顶点，则称M是一个匹配。简单地说，图7-2中x2、x3、x4等点都有多条边与之连接，也就是说有多个边依附于这些点，因此图7-2所示的二部图不是一个匹配。现在考虑删除一些边，最终得到如图7-3所示的一个G的子图。该子图中没有任何边同时连接X或Y中的同一个顶点，因此这是一个匹配。 如果G的一系列子图M0,M1,…,Mn都是匹配，那么包含边数最多的那个匹配就是图G的最大匹配。如果一个最大匹配中所有的点都有边与之相连，没有未覆盖点，则这个最大匹配就是完美匹配。未覆盖点的定义是：图G的一个顶点Vi，如果Vi不与任何一条属于匹配M的边相连，则成Vi是一个未覆盖点。图7-3就是一个完美匹配。 图7-2　简单的二部图 图7-3　一个完美匹配的二部图 根据以上定义，如果G的一个匹配M是最大匹配，并且没有未覆盖点，则这个匹配就是完美匹配。可见，图G的匹配和完美匹配正好就是之前介绍的“稳定婚姻问题”中的匹配和完美匹配。用图论的方法寻找完美匹配，需要首先找到最大匹配，当二部图中两个顶点集合中的顶点个数相等时，这个最大匹配同时也是完美匹配。求二部图的最大匹配可以使用最大流（maximal flow）或匈牙利算法（Hungarian algorithm），接下来我们就来介绍匈牙利算法。 7.4.1　最大匹配与匈牙利算法 寻找二部图最大匹配的匈牙利数学家埃德蒙德斯（Edmonds）在1965年提出的一个简化的最大流算法。该算法根据二部图匹配这个问题的特点将最大流算法进行了简化，提高了效率。普通的最大流算法一般都是基于带权网络模型的，二部图匹配问题不需要区分图中的源点和汇点，也不关心边的方向，因此不需要复杂的网络图模型，这就是匈牙利算法简化的原因。正是因为这个原因，匈牙利算法成为一种很简单的二分匹配算法，其基本流程是： 将图G最大匹配初始化为空 while(从Xi点开始在图G中找到新的增广路径) { 将增广路径假如到最大匹配中; } 输出图G的最大匹配; 根据匈牙利算法的流程看，寻找图G中的增广路径（Augment Path）是匈牙利算法的关键。先来看看什么是增广路径，二部图中的增广路径具有以下性质：  路径中边的条数是奇数；  路径的起点在二部图的左半边，终点在二部图的右半边；  路径上的点一个在左半边，一个在右半边，交替出现，整条路径上没有重复的点；  只有路径的起点和终点都是未覆盖的点，路径上其他的点都已经配对；  对路径上的边按照顺序编号，所有奇数编号的边都不在已知的匹配中，所有偶数编号的边都在已知的匹配中；  对增广路径进行“取反”操作，新的匹配数就比已知匹配数增加一个，也就是说，可以得到一个更大的匹配。 所谓的增广路径取反操作，就是把增广路径上奇数编号的边加入到已知匹配中，并把增广路径上偶数编号的边从已知匹配中删除。每做一次“取反”操作，得到的匹配就比原匹配多一个。匈牙利算法的思路就是不停地寻找增广路径，增加匹配的个数，当不能再找到增广路径时，算法就结束了，得到的匹配就是最大匹配。 增广路径的起点总是在二部图的左边，因此寻找增广路径的算法总是从一侧的顶点开始，逐个顶点搜索。从Xi顶点开始搜索增广路径的流程如下： while(从Xi的邻接表中找到下一个关联顶点Yj) { if(顶点Yj不在增广路径上) { 将Yj加入增广路径; if(Yj是未覆盖点或者从与Yj相关连的顶点（Xk）能找到增广路径) { 将Yj的关联顶点修改为Xi； 从顶点Xi开始有增广路径，返回true； } } 从顶点Xi开始没有增广路径，返回false； } 在这个算法流程中，“从与Yj相关连的顶点（Xk）能找到增广路径”这一步体现的是一个递归过程。因为如果之前的搜索已经将Yj加入到增广路径中，说明Yj在X集合中一定有一个关联点，我们假设Yj在X集合中的这个关联点是Xk，所以要从Xk开始继续寻找增广路径。当从Xk开始的递归搜索完成后，通过“将Yj的关联顶点修改为Xi”这一步操作，将其与Xi连在一起，形成一条更长的增广路径。 到现在为止，匈牙利算法的流程已经很清楚了，现在我们来给出实现代码。首先定义求最大匹配的数据结构，这个数据结构要能表示二部图的边的关系，还要能体现最终的增广路径结果，我们给出如下定义： typedef struct tagMaxMatch { int edge[UNIT_COUNT][UNIT_COUNT]; bool on_path[UNIT_COUNT]; int path[UNIT_COUNT]; int max_match; }GRAPH_MATCH; edge是顶点与边的关系表，用来表示二部图，on_path用来表示顶点Yj是否已经在当前搜索过程中形成的增广路径上了，path是当前找到的增广路径，max_match是当前增广路径中边的条数，当算法结束时，如果max_match不等于顶点个数，说明有顶点不在最大增广路径上，也就是说，找不到能覆盖所有点的增广路径，此二部图没有最大匹配。从Xi寻找增广路径的算法实现如下： bool FindAugmentPath(GRAPH_MATCH *match, int xi) { for(int yj = 0; yj < UNIT_COUNT; yj++) { if((match->edge[xi][yj] == 1) && !match->on_path[yj]) { match->on_path[yj] = true; if( (match->path[yj] == -1) || FindAugmentPath(match, match->path[yj]) ) { match->path[yj] = xi; return true; } } } return false; } 算法实现基本上是按照之前的算法流程实现的，不需要做特别说明，唯一需要注意的是path中存放增广路径的方式。读者可能已经注意到了，存放的方式是以Y集合中的顶点为索引存放，其值是对应的关联顶点在X集合中的索引。搜索是按照X集合中的顶点索引进行的，增广路径以Y集合中的顶点为索引存储，关系是反的。输出结果的时候，需要结合Y集合中的顶点索引输出，如果需要以X集合的顺序输出结果，需要反向转换，转换的方法非常简单： int path[UNIT_COUNT] = { 0 }; for(int i = 0; i < match->max_match; i++) { path[match->path[i]] = i; } 转换后path中就是以X集合的顺序存放的结果。 结合之前给出的匈牙利算法基本流程，最后给出匈牙利算法的入口函数实现： bool Hungary_Match(GRAPH_MATCH *match) { for(int xi = 0; xi < UNIT_COUNT; xi++) { if(FindAugmentPath(match, xi)) { match->max_match++; } ClearOnPathSign(match); } return (match->max_match == UNIT_COUNT); } 每完成一个顶点的搜索，需要重置Y集合中相关顶点的on_path标志，ClearOnPathSign()函数就负责干这个事情。 我们用图7-2中的二部图数据初始化GRAPH_MATCH中的顶点关系表edge，然后调用Hungary_Match()函数得到一组匹配： X1<--->Y3 X2<--->Y1 X3<--->Y4 X4<--->Y2 X5<--->Y5 结果与图7-3一致，因为这个最大匹配没有未覆盖点，所以是完美匹配。 匈牙利算法的实现以顶点集合V为基础，每次X集合中选一个顶点Xi做增广路径的起点搜索增广路径。搜索增广路径需要遍历边集E内的所有边，遍历方法可以采用深度优先遍历（DFS），也可以采用广度优先遍历（BFS），无论什么方法，其时间复杂度都是O(E)。匈牙利算法每个顶点Vi只能选择一次，因此算法的整体时间复杂度是O(V*E)，总的来说，是一个相当高效的算法。除了匈牙利算法，求二部图的最大匹配还可以使用Hopcroft-Karp算法。Hopcroft-Karp算法是由Hopcroft和Karp在1972年提出的一种算法，也是最大流算法的一种改进算法。算法的基本思想就是在每次搜索增广路径的时候不是只找一条增广路径，而是同时找几条互不相交的增广路径，形成最大增广路径集合，然后沿着集合中的几条增广路径同时扩大增广路径长度。通过进一步的分析，Hopcroft-Karp算法的时间复杂度可以达到O(sqrt(V)*E)，也非常高效。 7.4.2　带权匹配与Kuhn-Munkres算法 上一节我们介绍了二部图的最大匹配算法，用匈牙利算法寻找最大匹配，不要求每个个体给出的偏爱列表包含全部异性成员。比如在舞伴问题中，如果Albert非常讨厌Marcy，那么Albert的偏爱列表无论如何也不会包含Marcy。在这一节，我们让舞伴问题再复杂一点，于是为舞伴问题引入带权优先表的概念，为每一个配对指定一个权重，表明我们更希望哪一对成为舞伴。通过控制每一对舞伴关系的权重，使得最后的完美匹配结果中有尽量多的舞伴是我们所期望的配对关系。 这个问题变得有点像最优解问题了，一提到与图有关的最优解问题，你会想到穷举法。穷举所有的完美匹配，然后计算每个完美匹配中各边的权重之和，取权重之和最大的一个作为最后的结果。这是一种解决方案，但是穷举法虽然是万能方法，但是不到万不得已最好不要用穷举法。仔细思考一下，其实这个问题已经演化成了求解二部图的带权匹配问题，所谓二部图的带权匹配其实就是求出一个匹配集合，使得集合中各边的权值之和最大或最小。对于本问题，给每一个配对（图中的边）指定一个权重之后问题就变成了求二部图的带权最大匹配问题。 通过之前对Gale-Shapley算法和匈牙利算法的介绍，我们已经了解了完美匹配、稳定匹配和最大匹配这些概念，那么带权匹配和之前的这些概念是什么关系呢？答案是没有半毛钱关系，至少和完美匹配与最大匹配之间不存在包含或等于关系。二部图的最大权或最小权匹配，只是要求得到的一个匹配中各边的权值之和最大或最小，并不要求这个匹配是完美匹配或最大匹配。如果这个权值最大（或最小）的匹配同时又是完美匹配，则这样的结果就被称为最佳匹配。本节我们要介绍的Kuhn-Munkres算法是求最大权或最小权匹配的算法，如果期望Kuhn-Munkres算法得到的结果同时是一个完美匹配（最佳匹配），那么要求算法运行的数据必须存在完美匹配（比如两个顶点集合的顶点个数必须相等之类的条件）。很多同学会忽略这一点，以为Kuhn-Munkres算法可以在任何情况下得到带权的最大匹配，这个理解是错误的。 Kuhn-Munkres算法，也称KM算法，是Kuhn和Munkres二人在1955～1957年各自独立提出的一种算法，是一种求解最大最小权匹配问题的经典算法。最初的Kuhn-Munkres算法以矩阵为基础结构，但是Edmonds在1965年发布了匈牙利算法之后，Kuhn-Munkres算法也基于匈牙利算法进行了改进。当给定的二部图存在完美匹配的情况下，Kuhn-Munkres算法通过给每个顶点设置一个标号（叫作顶标）的方式把求最大权匹配的问题转化为求完美匹配的问题的，最终得到一个最大权完美匹配。那么这个转换是如何实现的呢？这就需要分析一下Kuhn-Munkres算法的原理了。 我们假设二部图中X顶点集合中每个顶点Xi的顶标是A[i]，Y顶点集合中每个顶点Yi的顶标是B[i]，顶点Xi和Yi之间的边的权重是weight[i][j]，则Kuhn-Munkres算法的原理就是基于以下定理： 若由二部图中所有满足A[i]+B[j] = weight[i,j]的边(Xi,Yj)构成的子图（称作相等子图）有完美匹配，那么这个完美匹配就是二分图的最大权匹配。 现在明白转换原理了吧，就是先找出问题对应的相等子图，然后求相等子图的完美匹配即可。现在的问题是，这个定义成立吗？答案是只要在算法过程中始终满足“A[i]+B[j] ≥ weight[i,j]”这个条件，这个定理就成立。因为对于二分图的任意一个匹配，如果这个匹配是相等子图的匹配，那么它的边权重之和等于所有顶点的顶标之和（显然这是最大的）；如果这个匹配不是相等子图的匹配（它的某些边不属于相等子图），那么它的边权重之和小于所有顶点的顶标和。所以只要始终满足“A[i]+B[j]≥ weight[i,j]”条件的相等子图的完备匹配一定是二分图的最大权匹配。 根据以上分析可知，Kuhn-Munkres算法的实现流程大致如下所示。 (1) 初始化各个顶点的顶标值； (2) 找出符合“A[i]+B[j] = weight[i,j]”条件的边构成相等子图，使用匈牙利算法寻找相等子图的完美匹配； (3) 如果找到相等子图的完美匹配，则算法结束，否则调整相关顶点的顶标值； (4) 重复步骤(2)(3)，直到找到完美匹配为止。 第(1)步初始化顶点顶标值可采用式(7-1)计算： (7-1) 因为A[i]总是取与之相邻的边中最大的权值作为初始值，因此初始阶段能保证满足“A[i]+B[j] ≥ weight[i,j]”条件。如果在第(2)步的相等子图中没有找到完美匹配，说明相等子图中某个顶点出发的增广路径不能覆盖所有顶点。此时需要调整各个顶点的顶标值，然后重新在相等子图中寻找完美匹配。调整顶标的目的是为了扩大相等子图，使得更多的边进入相等子图，并最终能够找到一个完美匹配。设当前增广路径上所有属于X集合的顶点构成一个子集S，所有属于Y集合的顶点构成一个子集T，dx为顶标调整的变化量，则dx可采用式(7-2)给出的方法计算： (7-2) 由dx的计算公式可知，如果把S集合中所有顶点的顶标都减少dx，一定会有一条一端在S中，另一端不在T中的边因满足“A[i]+B[j] = weight[i,j]”的条件而进入相等子图，这就扩大了相等子图。对S集合中所有顶点的顶标都减少dx之后，为了使原来已经在相等子图中的边继续留在相等子图中，需要将T集合中所有顶点的顶标值增加dx，使A[i]+B[j]之和不变。 现在总结一下顶标调整的方法，首先采用式(7-2)计算出调整变化量dx，然后将S集合中所有顶点的顶标值减少dx，同时将T集合中所有顶点的顶标值增加dx，这样的调整，对整个图上的所有顶点会产生如下四种结果。  对于两端点都在当前相等子图的增广路径上的边（xi, yj），其顶标值A[i]+B[j]的和没有变化。也就是说，原来属于相等子图的边，调整后仍然属于相等子图。  对于两端点都不在当前相等子图的增广路径上的边（xi, yj），其顶标值A[i]和B[j]的值没有变化。也就是说，此边与相等子图的隶属关系没有变化，原来属于相等子图的边现在仍然属于相等子图，原来不属于相等子图的边现在仍然不属于相等子图。  对于xi在当前相等子图的增广路径上，yj不在当前相等子图的增广路径上的边（xi, yj），其顶标值A[i]+B[j]的和略有减小，原来不属于相等子图，现在有可能属于相等子图，使得相等子图有机会得到扩大。  如果xi不在当前相等子图的增广路径上，yj在当前相等子图的增广路径上的边（xi, yj），其顶标值A[i]+B[j]的和略有增加，原来不属于相等子图，现在仍不属于相等子图。 由此可见，每次调整顶标，都能在图的基本状态保持不变的情况下扩大相等子图，使得相等子图有机会找到一个完美匹配，这就是顶标值调整在算法中的意义所在。 现在，我们就结合一个带权最大匹配问题的实例，给出这个算法在实际应用中的一个实现。问题是这样的，某公司有5名技术工人，他们都可以完成公司流程中的5种工作，但是每个工人的技术侧重点不同，熟练程度也不同，因此他们完成同样的工作所产生的经济效益也不相同。如果用0～5范围内的值对每个工人完成每种工作所产生的经济效益进行评价，可得到如表7-1所示的经济效益矩阵。假如你是这家公司的负责人，你需要找到一种工人和工作之间的匹配关系，使得这种匹配关系能产生的经济效益最大。根据之前对Kuhn-Munkres算法的分析，我们针对这个问题设计了KM_MATCH匹配数据结构，如下所示： typedef struct tagKmMatch { int edge[UNIT_COUNT][UNIT_COUNT]; //Xi与Yj对应的边的权重 bool sub_map[UNIT_COUNT][UNIT_COUNT];// 二分图的相等子图, sub_map[i][j] = 1 代表Xi与Yj有边 bool x_on_path[UNIT_COUNT]; // 标记在一次寻找增广路径的过程中，Xi是否在增广路径上 bool y_on_path[UNIT_COUNT]; // 标记在一次寻找增广路径的过程中，Yi是否在增广路径上 int path[UNIT_COUNT]; // 匹配信息，其中i为Y中的顶点标号，path[i]为X中顶点标号 }KM_MATCH; 相对于匈牙利算法中的GRAPH_MATCH定义，KM_MATCH的主要变化是增加了sub_map作为相等子图定义和标识yi是否在增广路径上的y_on_path标识。相对于我们前面对Kuhn-Munkres算法的分析，edge对应于边的权重表weight，sub_map对应于算法执行过程中的相等子图，x_on_path和y_on_path分别用于标识X集合和Y集合中的顶点是否属于增广路径上的S集合和T集合，path就是最后匹配的结果。 表7-1　不同工人完成不同工作的经济效益 工作1 工作2 工作3 工作4 工作5 工人1 3 5 5 4 1 工人2 2 2 0 2 2 工人3 2 4 4 1 0 工人4 0 1 1 0 0 工人5 1 2 1 3 3 下面给出Kuhn-Munkres算法的具体实现代码，Kuhn_Munkres_Match()函数虽然很长，但是并不难理解，因为这个代码是严格按照之前给出的Kuhn-Munkres算法的流程实现的。包括顶标的初始化、使用匈牙利算法求完美匹配和顶标调整在内的三个主要算法步骤在Kuhn_Munkres_Match()函数中都得到体现，并且界定非常清晰。其中寻找增广路径的FindAugmentPath()函数与之前介绍匈牙利算法时给出的FindAugmentPath()函数实现非常类似，区别就是使用sub_map而不是直接使用edge，并且额外记录了x_on_path标识。ResetMatchPath()函数负责每次开始寻找相等子图之前清除上一次搜寻产生的临时增广路径，ClearOnPathSign()函数负责在每次搜寻增广路径之前清除顶点是否属于S集合或T集合的标识，大家可以从本书的配套代码中找到此函数的代码。 bool Kuhn_Munkres_Match(KM_MATCH *km) { int i, j; int A[UNIT_COUNT], B[UNIT_COUNT]; // 初始化Xi与Yi的顶标 for(i = 0; i < UNIT_COUNT; i++) { B[i] = 0; A[i] = -INFINITE; for(j = 0; j < UNIT_COUNT; j++) { A[i] = std::max(A[i], km->edge[i][j]); } } while(true) { // 初始化带权二分图的相等子图 for(i = 0; i < UNIT_COUNT; i++) { for(j = 0; j < UNIT_COUNT; j++) { km->sub_map[i][j] = ((A[i]+B[j]) == km->edge[i][j]); } } //使用匈牙利算法寻找相等子图的完备匹配 int match = 0; ResetMatchPath(km); for(int xi = 0; xi < UNIT_COUNT; xi++) { ClearOnPathSign(km); if(FindAugmentPath(km, xi)) match++; else { km->x_on_path[xi] = true; break; } } //如果找到完备匹配就返回结果 if(match == UNIT_COUNT) { return true; } //调整顶标，继续算法 int dx = INFINITE; for(i = 0; i < UNIT_COUNT; i++) { if(km->x_on_path[i]) { for(j = 0; j < UNIT_COUNT; j++) { if(!km->y_on_path[j]) dx = std::min(dx, A[i] + B[j] - km->edge[i][j]); } } } for(i = 0; i < UNIT_COUNT; i++) { if(km->x_on_path[i]) A[i] -= dx; if(km->y_on_path[i]) B[i] += dx; } } return false; } 根据表7-1提供的数据初始化KM_MATCH数据结构，然后调用Kuhn_Munkres_Match()函数，得到一个最大权匹配的结果，因为原数据存在完美匹配，因此这个结果就是最佳匹配结果： 工人 1 分配 工作 3 （经济效益评价是5） 工人 2 分配 工作 1 （经济效益评价是2） 工人 3 分配 工作 2 （经济效益评价是4） 工人 4 分配 工作 5 （经济效益评价是0） 工人 5 分配 工作 4 （经济效益评价是3） 最后获得最大经济效益评价是14。需要说明的是，对于同一个问题，其最大权匹配的结果可能不唯一，也就是说，存在多个匹配的权重之和同为最大值的情况。Kuhn-Munkres算法可以找出其中的一个，但是无法找到全部匹配结果。