KM算法解决不同行不同列求和最大值问题

  昨天在大学专业的交流群里看到某个学弟发了张自己打算学的算法列表图，然后瞟了一眼，看到了KM算法，而之前写的一篇深信服在线笔试-2018.9.22，的编程题的第四题：求一个二维数组或者说矩阵的不同行不同列的合的最大值，我采用的是穷举方法，而看到别人说最好的方式是采用KM算法，穷举计算量实在过大了，所以来学习一下KM算法。
  （PS：这段话是闲扯）然后就去百度了一下KM算法，看到了很多篇大多都是先介绍匈牙利算法，然后KM算法，匈牙利算法让我突然想起了大四上学期学习的运筹学（忘光了），然后脑子抽了一下，突然想起可不可以用动态规划去解决之前的那个问题，即把问题模拟成一个最长路径单向图问题，这样的话就需要建立一个路径表，但是这是不可行的，因为动态规划有一个性质，如果某阶段状态给定后，则在这个阶段以后过程的发展不受这个阶段以前各段状态的影响，而不同行不同列的值计算，把每一列看成一个阶段，那么每个阶段会依赖于之前的各个段影响（受不同行的影响），那么就不适用于动态规划，所以还是老老实实学习KM算法怎么用吧。
  
  
  由于写的有些多，在这里大概介绍文章结构分布：
    1、KM算法步骤以及计算过程
    2、代码实现思路
    3、代码
  
  个人理解：KM算法，用通俗的话说其过程，其实就是对二分图中的两个顶点集合X、Y中的每个元素设置标杆（初始值），然后给其中的一个集合的所有元素标杆初始化为连接到另一个顶点的边的最大（小）权值，然后另一个集合的所有标杆初始化为0，这样二分图中的每个顶点与顶点的一开始的连接都是最大权值的边，为了让每个顶点都只匹配一个顶点（在我们解决的问题中就是解决不同行不同列），我们就需要对连接的边进行调适，让冲突的边（X集合中多个顶点连接Y集合中的一个顶点）进行换边，使换边过程造成的权值变化最小的顶点更换连接然后更新标杆值，通过这样的操作，将每个点与点进行匹配，然后不断的进行换边，从而达到全局最优（最大或最小）的匹配结果。

KM算法：

  用于求二分图匹配的最佳匹配。何为最佳匹配？就是带权二分图的权值最大的完备匹配称为最佳匹配。 那么何为完备匹配？X部中的每一个顶点都与Y部中的一个顶点匹配，或者Y部中的每一个顶点也与X部中的一个顶点匹配，则该匹配为完备匹配。
  二分图：就是一个图中的顶点集V可分割为两个互不相交的子集，并且图中每条边依附的两个顶点都分属于这两个互不相交的子集（如X、Y），两个子集内的顶点不相邻。
  这里我会一步一步的说明的，相关的名词也会慢慢介绍，对于增广路这个概念，我个人觉得在代码中没体现出来，所以就没打算说这个事，如果使用增广路，然后取反得到完备匹配的话，那么个人感觉是需要设置两种状态，一种是可匹配（未匹配），另一种是已匹配，然后取反操作，似乎这样就变得有些难理解了，所以在这里并不打算说这个。

算法步骤如下：
  1.用邻接矩阵（或其他方法也行啦）来储存图，注意：如果只是想求最大权值匹配而不要求是完全匹配的话，请把各个不相连的边的权值设置为0。
  2.运用贪心算法初始化标杆。
  3.运用匈牙利算法找到完备匹配（遍历X集合的每个顶点）。
  4.如果找不到，则通过修改标杆，增加一些边。
  5.重复3，4的步骤，直到完全匹配时可结束。
  
  在这里设cx[i]为X集合的标杆，cy[i]为Y集合的标杆，w[i][j]为图的邻接矩阵表示X[i]–Y[j]边的权值。
  标杆：是为了求最大（小）权值而引入的一个概念，用于存储X集合每个顶点到Y顶点中的最大值，而Y集合的标杆值全初始化为0，如下图中，X集合为{A,B,C}，而Y集合为{D,E,F}，每个X的顶点标杆值为路径权值的最大值，而Y的顶点的标杆值初始化为0.

  重点：因为使用了标杆这个概念，我们在使用匈牙利算法进行边匹配的时候就不能是看边是否存在，也就是判断w[i][j]==0来判断边是否存在，引入标杆后，我们判断边是否存在（这里指的是匹配的两个标杆的值的合是否与路径值相等，不相等就不存在），即cx[i]+cy[j]-w[i][j]==0,如果该等式成立的话，就表明这两个顶点可以匹配（连接），但并不表明这两个顶点可以直接匹配（连接），因为可能还有其他点与其连接，所以在这种条件下，存在两种情况：
  1、如果这个成立且Y[j]边没有被X集合中其他顶点匹配。
  2、Y[j]边已经被X集合中的某一顶点X[k]匹配，但是X[k]可以匹配其他的顶点，或者是X[i]可以改选其他的匹配。
  只要满足以上两个条件之一，X[i]与Y[j]就可以匹配，然后这里问题又来了，如果碰上情况2的话，如何让X[k]去匹配其他顶点呢？这时候，我们就要靠d值去更新标杆值，d=min(cx[i]+cy[j]-w[i][j])，这时候标杆的作用就体现出来了，不太了解是什么意思的话，我这里举一个例子：
  在上图中，首先A与D匹配，记作A-D，然后B的标杆值为14，而14+0=14，所以B也要去匹配D，此时就碰上条件2的问题，那么我们该让A、B两点谁去匹配其他点呢？，这时候就得根据d值进行选择，从图上，我们可以看出d=cx[A]+cy[E]-w[A]][E]=3，是A到除了被匹配的点以外的其他点（即除了D点以外的点）权值变化最小的d值，而B到除了D点以外的点的变化最小的d值为cx[B]+cy[F]-w[B][F]=6，这里我们就可以看出，如果要改动B的话，需要变化（或者说是变动）的权值为6，相比起改动A的匹配点仅仅需要变化3来说，改动B是不划算（对于求最大权值来说是不划算），所以我们就更新各个标杆的值，让A可以去匹配E，即A-E，B-D，当然，这里并没完成匹配，只是更新标杆值。
  到这里，应该清楚如何添加新的边了，添加完新的边后，对应的标杆值也需要更新，即A、B的标杆值-d，D的标杆值+d，这里我们就需要清楚，到底是哪些值会受影响，如上面受影响的是集合X中的{A、B}与Y集合中的{D}，也就是当我们遍历B点匹配的时候，匹配形成的路线是B-D-A，是一条不合法（冲突）的路径，此时设置一个集合S存储不合法路径中的X集合元素A、B，集合T存储不合法路径中的Y集合元素D，这样每次更新d值的时候，遍历S集合进行-d，遍历T集合进行+d即可。
  上面算是对步骤3、4进行详解了，接下来可以来手工模拟上图的求值步骤。
  1、遍历X集合，首先，A可匹配D，D没有被匹配过，A匹配D。
  2、B可匹配D，而A已经匹配D，形成B-D-A不合法路径，此时S集合为{A、B}，T集合为{D}，那么从A、B两者中找匹配E、F变动权值最小的点进行更改匹配，经过d值比较得出A匹配E变动较小，d=3，所以更新S集合与T集合中的标杆值，值的变化使得A可以匹配E，也仍旧可以匹配D。
  3、更新完值后清空集合S、T，再次遍历B点匹配，然后B与D匹配，A与E匹配。
  4、C无可匹配的点，此时S集合为{C}，T集合为{}，那么从C匹配D、E、F中找d值变动最小的进行匹配，此时D标杆值为3，所以d1=cx[C]+cy[D]-w[C][D]=3，d2=cx[C]+cy[E]-w[C][E]=1，d3=3，选取最小的d2=1，然后更新集合S中元素的值，此时S集合只有C，所以C标杆值=cx[C]-d=12。
  5、因为C没有匹配到值，所以还得继续匹配，此时C可以匹配E，而A已经匹配E了，形成C-E-A-D-B不合法路径，此时S集合为{A、B、C}，T集合为{D、E}，这时候找A、B、C到点F变动最小得值d，然后再更新S集合与T集合各顶点的标杆，在这里d=2，即C可匹配F变动最小。
  6、更新完值后清空集合S、T，C标杆值为10，再次为C进行匹配，此时F可匹配，因为X集合点遍历完成，所以程序结束。
  7、经过上面的匹配，得到的结果是A-E，B-D，C-F，则最大值就是12+14+10=36.
  
  

代码实现思路：

  
  下面是伪代码，和实际代码其实差别不大…
  findpath函数，即找point到Y集合中每个点是否存在匹配路径，S、T则负责记录经过的路径用于在找不到匹配路径的情况下添加新的边，以及防止进入死递归。
  

//邻接矩阵的行与列值
int row,col;
//记录X-Y顶点匹配的数组，即link[y]=x，表示集合Y的第y个顶点与X集合的第x顶点匹配
int link[MAXLINE]={0}；
//存储图的权值的邻接矩阵
int w[MAXLINE][MAXLINE]={0};
//存储X、Y标杆权值的数组
int cx[MAXLINE]={0},cy[MAXLINE]={0};

bool findpath(int point,char *S,char *T){
	//路径上X集合的点，加入S集合中
	S[point]=1;
	for y=1;y<=col;y++
		if T[y]==0 and cx[point]+cy[y]==w[point][y]{
		//point与y顶点可匹配，把y加入集合T中
			T[y]=1;
		//判断是否y顶点是否与其他顶点点k连接
		//如果存在点k那就试图找其他路径
			if link[y]==0 || findpath(link[y],S,T){
				link[y]==point
				return true;
			}
		}
	return false;	
}

int KM(){
	//初始化标杆
	for i =1;i<=row;i++
		cx[i]=w[i][1];
		for j =2;j<col;j++
			cx[i]=max(cx[i],w[i][j])
	//从X集合种的第一个顶点开始匹配
	int Xpoint=1; 
	int d=INF;
	while(Xpoint<=row){
		//初始化集合S与T，存放不合法路径上的顶点，以便修改标杆值
		char *S=new char[row](0);
		char *T=new char[col](0);
		//findpath为使用Xpoint遍历Y集合的顶点函数
		//return true表示Xpoint匹配成功，进行下一个点匹配
		//return false表示匹配不成功，需要加一条边或者说是让S集合中的某个顶点改选路径
		if (findpath(Xpoint,S,T)){
			Xpoint++；
			delete []S;
			delete []T;
			continue;
		}
		//找出集合S中匹配除了T集合中的顶点的其他顶点权值变动最小的值
		for x =1;x<=row;x++
			if S[x]==1:
				for y=1;y<=col;y++
					if Y[y]!=1:
						d=min(d,cx[x]+cy[y]-w[x][y])
		//没找到匹配点
		if d==INF
			return -1；
		//更新标杆值
		for x=1;x<=row;x++
			if S[x]==1:
				cx[x]-=d;
		for y=1;y<=col;y++
			if T[y]==1:
				cy[y]+=d;
		delete []S;
		delete []T;
	}
	int res=0;
	for i=1;i<=row;i++
		int flag=link[i];
		if flag>0:
			res+=w[flag][i]
	return res;
}

  
  

代码

  
实际代码采用的是别人写的，出自于：二分图（三）——KM算法
  个人仅仅修改了部分代码以及注释，因为逻辑挺清楚的，就不打算重写了，这里也挺佩服这些算法竞赛的人，初中高中就在学这些算法，想想我那时候还在网吧玩游戏…哈哈…

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;

const size_t SIZE = 1000;
int w[SIZE][SIZE];
const int inf = (1 << 20) - 1;
int m, n;
int cx[SIZE], cy[SIZE];
bool usex[SIZE], usey[SIZE];//用于记录每轮的集合S与T
int link[SIZE];//link[i]=x代表集合Y[i]中的顶点与x匹配或者连接

//dfs 1、return false，用于记录集合S与T，通过深度遍历，找到连接路径上的点。
//2、retrurn true，用于记录匹配，或者对已匹配的点进行再匹配。
bool dfs(int u) {
	usex[u] = 1;
	for (int i = 1; i <= m; i++)
		if (!usey[i] && cx[u] + cy[i] == w[u][i]) {
			usey[i] = 1;
			if (link[i] == 0 || dfs(link[i])) {
				link[i] = u;
				return 1;
			}
		}
	return 0;
}
int KM() {
	memset(cy, 0, sizeof(cy));
	memset(cx, 0, sizeof(cx));
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= m; j++)
			cx[i] = max(cx[i], w[i][j]);
	for (int i = 1; i <= n; i++) {
		while (1) {
			int d = inf;
			memset(usex, 0, sizeof(usex));
			memset(usey, 0, sizeof(usey));
			if (dfs(i))
				break;
			//以下为对S集合找变动最小的d值
			for (int i = 1; i <= n; i++) {
				if (usex[i]) {
					for (int j = 1; j <= m; j++)
						if (!usey[j])
							d = min(d, cx[i] + cy[j] - w[i][j]);
				}
			}
			if (d == inf)
				return -1;
			//更新S、T集合中的标杆值。
			for (int i = 1; i <= n; i++)
				if (usex[i])cx[i] -= d;
			for (int i = 1; i <= m; i++)
				if (usey[i])cy[i] += d;
		}
	}
	int ans = 0;
	for (int i = 1; i <= m; i++) {
		int flag = link[i];
		if (flag)
			ans += w[flag][i];
	}
	return ans;
}
int main() {
	while (scanf("%d%d", &n, &m)) {
		memset(w, 0, sizeof(w));
		memset(link, 0, sizeof(link));
		for (int i = 1; i <= n; i++)
			for (int j = 1; j <= m; j++)
				scanf("%d", &w[i][j]);
		printf("%d\n", KM());
	}
	return 0;
}