前面介绍了二分图上求最大匹配的匈牙利算法，这里介绍在带权二分图上求最优匹配的KM算法。
最优匹配
对于赋权完全偶图，具有最大权值的完美匹配就是 最优匹配。有时候权值表示相似度之类的东西，像是DETR里面根据IoU来分配检测框，就要求最大。有时候权值的含义是成本、开销，这个时候是求最小权值的完美匹配，这里两种情况都会介绍，不过以求最大的为主。为了区分求最小和最大的情形，有时候也会把最优匹配称为最小权匹配或最大权匹配。
本文用$G$表示二分图，二分图的两个顶点集合用$X$和$Y$表示，权重矩阵用$W$表示，$w_{xy}$表示边$xy$的权重。在求最优匹配的时候只用考虑完全二分图，这个时候最优匹配一定存在。除此之外假设$|X|=|Y|$，因为不等的情况可以通过添加一些顶点以及权为0的边转化为相等的情况。
KM算法
Kuhn-Munkres算法，也叫匈牙利算法，二分图里求最大匹配、最优匹配的算法都叫匈牙利算法。。。不过我个人更喜欢叫KM算法。求最大匹配的匈牙利算法是最大流问题的特例，求最小权匹配的KM算法则是最小费用流的特例。
相等子图
介绍KM算法需要引入两个概念，一个概念是 可行顶点标号，它指的是满足$\forall x\in X,y\in Y,l(x)+l(y)≥w_{xy}$的、关于顶点的函数$l$，本文中简写为 标号。可行顶点标号也叫顶标、期望度、势函数。另一个概念是 相等子图，它是一个根据$l$得到的原图$G$的边导出子图。记$E_l=\lbrace xy\in E(G)|l(x)+l(y)=w_{xy}\rbrace$，称$G_l=G(E_l)$是二分图$G$对应于标号$l$的相等子图，这里$G(E_l)$表示$E_l$这个边集在$G$里对应的子图。
这两个概念看起来很抽象，但其实结合例子来看非常直观。考虑下面这样一个完全二分图：
它的可行顶点标号和相等子图分别是：
可行顶点标号就是对顶点赋一个值，每个图都存在标号，最常见的标号是把$X$里的顶点赋值为邻接边里面最大的权值，而$Y$里的全部设为0，即
\[\begin{aligned} l(x) &\leftarrow \min_{y} w_{xy} \\ l(y) &\leftarrow 0 \end{aligned}\]相等子图就是对边做了筛选，选出来的是权重比较大的那些边。其实标号和相等子图是通过对偶理论推出来的，标号对应对偶变量，相等子图是等式约束的一种体现。
Kuhn-Munkres定理
在相等子图的基础上有一个KM定理，它是KM算法的基础：
若相等子图$G_l$存在完美匹配$M^* $，则 $M^* $是$G$的最优匹配。
这个定理很容易证明，因为标号和是任意一个完美匹配的权重之和的上界，记一个匹配的权重之和为${\rm val}(M)$：
\[{\rm val}(M)=\sum_{x,y\in M} w_{xy} \le \sum_{x,y\in M} [l(x)+l(y)] \le \sum_{x}\sum_{y} [l(x)+l(y)]\]对于相等子图上匹配$M$，第一个不等号取等；如果$M$还是完美匹配，那么第二个等号也取等，所以这个匹配等于最小标号和，因而也是所有完美匹配的上界，所以是最优匹配。
基于这个定理，我们可以把求最优匹配转化为求完美匹配，这样就可以利用匈牙利算法求解了。但是这个前提是我们要找到一个合适的标号，使得它的相等子图存在完美匹配才行。这个可以通过初始化一个$l$，然后不断调整$l$实现。初始化就用上面例子里那种取法就可以了，那么如何调整呢？考虑Hall定理，不存在完美匹配是因为存在冲突，而冲突的地方是找不到交错路的子图。如果在冲突的地方添加更多的边，就可以解决冲突。KM算法通过对冲突的顶点集合$S$和$T$的标号进行调整，使得调整后的$l$的相等子图上比原来多出一些边来解决冲突。通过不断地调整，就能得到一个$l$，它的相等子图$G_l$存在完美匹配。
Kuhn-Munkres算法
下面给出KM算法的流程：
初始化标号lx和ly
for x in range(n):
while True: # while循环用来给顶点x找增广路，找不到就调整标号，直到可以找到为止，再去处理下一个顶点
# O(n^2)
在相等子图上，从x开始找一条增广路P # 可以根据lx[x] + ly[y] == W[x][y]判断一条边是否属于相等子图
if P存在:
扩大匹配
break
# O(n^2)
更新标号 # 根据 S和 T更新lx和ly
整体和最最大匹配类似，但是多了一个在找不到增广路时更新标号的步骤。在找不到增广路的时候，根据Hall定理，存在冲突的部分$S$，$T$，并且有$N_{G_l}(S)=T$。基于$S$和$T$按以下方法调整标号：
\[\begin{aligned} \delta_l & ← \min\{l(x)+l(y)-w_{xy}|x\in S,y\in Y-T\} \\ l(x) & ← \left\{ \begin{array}{} l(x)-\delta_l & x \in S\\ l(x)+\delta_l & x \in T\\ \end{array} \right. \end{aligned}\]简单来说就是先求$\delta_l$，然后$S$里的顶点都减小$\delta_l$，而$T$里的顶点都增大$\delta_l$。由于$|T|$比$|S|$小1，$T$里顶点标号增大而$S$里减小会导致整体减小$\delta_l$（所有顶点标号的和会减小$\delta_l$）。
这里给出一个例子帮助理解如何调整标号以及调整标号的作用（后面会给出详细的分析）：
一开始利用初始化的标号$l$得到一个相等子图，上面最大匹配数量为3，不存在完美匹配。根据Hall定理，原因在于顶点1,3都与6相连，冲突了。调整了标号后，相等子图多了2条边，冲突得到解决，可以得到完美匹配。
理解调整标号的作用
那么为什么需要这样调整标号？这样调整有什么影响？我们可以根据边所属的集合进行分类讨论来分析（一共有4种边）：
-
由于$S$里的顶点标号减小，$T$里的增加，但是和不变，所以它们之间连的边不受影响
-
$X-S$和$Y-T$之间的边显然也不受影响，因为对$S$和$T$的调整根本不影响它们
-
$X-S$和$T$之间的边会可能被移出$G_l$，因为$T$里的顶点标号增大了
有边被移出看起来会更加容易冲突，但是注意$X-S$里的顶点要么已经匹配了，要么还没处理，假设当前正在处理$x$，那么$x$之前的顶点必然已经匹配，而$x$之后的都还没处理。调整标号后，被移出的边肯定不是匹配边，所以已经匹配的边不受影响，而还没处理的顶点也压根不用考虑，因为后面才会给它们找匹配。
-
$S$和$Y-T$之间至少有一条边会被添加到$G_l$中，因为$S$里的顶点标号减小了而$Y-T$里的不变，这意味着$l(x)+l(y)$变小。由于$\delta_l$就是从$S$和$Y-T$之间的边里面取最小边权得到的，调整标号后，取到$\min$的边会使$l(x)+l(y)=w_{xy}$，因而被添加到相等子图里。这个结合前面的例子很好理解。
考虑增加的边，如果连接的是非饱和点，那下一次匈牙利算法就能找到可扩路了；如果是饱和点，他一定与$X-S$中的点匹配，下一次交错路会被延长，还需要继续调整标号才能找到增广路。
综上，调整标号后，一些边会被移除，但移除的边没什么影响；一些边会被添加，被添加的边一定是$S$和$Y-T$之间的边，添加边后是否能解决冲突，使得增广路出现，这取决于被添加的边是否连接一个非饱和点（这是一个可以优化的点）。
结合上面的理解，$G_l$其实在限制每个$x$能被分配到的$y$，一开始贪心地选最大的边，但是这样很容易找不到匹配。通过减去$\delta_l$就可以允许更多的边进来，这些新边的权值会更低一些，但是添加新的边可以解决冲突。$\delta_l$取$\min$保证了最优性，标号的调整是最小的。标号$l$的和可以理解为我们对$G_l$上匹配的预期，如果预期太高就没有饱和匹配，预期刚好等于最优匹配的话，就能得到最优匹配，这是上面那个定理的证明。
最后总结一下KM算法的逻辑，就是遍历$X$，不断给$x$找匹配，找不到说明需要给相等子图添加新的边，于是就去调整标号，调整后必有新的边加入，然后继续给$x$找匹配，可能需要添加多次才能给$x$找到匹配。等所有$x$都找到匹配了，就得到了最大匹配。和求最大匹配的匈牙利算法相比，只是多了调整标号这一步。
算法复杂度与最小权匹配
KM算法跟匈牙利算法类似，最外层循环遍历$X$，$X$有$n$个顶点要处理，每次都要找增广路和调整标号。找增广路最坏$O(n^2)$，最坏的情况下每个顶点调整标号$O(n)$次，相应的也要重复找增广路，所以总复杂度是$O(n^4)$。 标号的调整次数的多少这取决于二分图里是不是高权重的边都集中到几个顶点之间，每个顶点最多调整$O(n)$次的原因是每次调整必定使$T$增加一个点，最多运行$n$次$T$就等于$Y$了，这个时候肯定有非饱和点被添加到相等子图里，于是就有增广路了。
对于最小权值的情况，可以把令$W=\max(W)-W$转换成求最大的情况，也可以对算法稍作修改来求：把标号初始化的$\max$改成$\min$，然后调整标号的时候，是$S$里的增大，$T$里的减小。
KM算法的实现优化
计算delta需要遍历$S$和$T$，复杂度$O(n^2)$。我们可以把delta的计算从$O(n^2)$降低到$O(n)$，方法是引入一个slack数组，在找增广路的时候记录$\min_{x\in S} l(x)+l(y)-w_{xy}$，这样求delta就只需要对一维数组求$\min$了。
vector<int> lx, ly;
vector<int> matchY;
bool findAugPath(vector<vector<int>> &W, int x, vector<int> &visitedX, vector<int> &visitedY) {
int n = W.size();
// 从x出发找一条增广路
visitedX[x] = true;
for(int y=0; y<n; y++) {
if(visitedY[y]) continue;
if(lx[x] + ly[y] == W[x][y]) {
visitedY[y] = true;
if(matchY[y] == -1 || findAugPath(W, matchY[y], visitedX, visitedY)) {
// 这里有短路求值，合并了两种情况
matchY[y] = x; // 隐式的根据交错树修改匹配
return true;
}
} else {
// 在这里加刚好保证了 x属于S，y属于Y-T
slack[y] = min(slack[y], lx[x] + ly[y] - W[x][y]); // x邻接的y都要调整
}
}
return false;
}
int KM(vector<vector<int>> &W) {
int n = W.size(); // W是n*n的方阵
// 初始化变量
matchY = vector<int>(n, -1);
lx = vector<int>(n);
ly = vector<int>(n);
for(int i=0; i<n; i++) {
for(int j=0; j<n; j++) {
lx[i] = max(lx[i], W[i][j]);
ly[j] = 0;
}
}
// 开始迭代
for(int x=0; x<n; x++) {
while(true) {
// 从x开始迭代找增广路，O(n^2)，如果找到了就修改匹配
vector<int> visitedX(n), visitedY(n);
slack = vector<int>(n, numeric_limits<int>::max());
if(findAugPath(W, x, visitedX, visitedY)) {
break; // 找下一个x的增广路
}
// 计算delta, O(n)
int delta = numeric_limits<int>::max();
for(int y=0; y<n; y++) {
delta = min(delta, slack[y]);
}
// 调整标号, O(n)
for(int i=0; i<n; i++) {
if(visitedX[i]) lx[i] -= delta;
if(visitedY[i]) ly[i] += delta;
}
}
}
// 根据匹配计算结果
int ans = 0;
for(int y=0; y<n; y++) {
ans += W[matchY[y]][y];
}
return ans;
}
其实常见的slack不是每次迭代都重新初始化的，但是这么写更好理解，也不怎么影响时间复杂度。
经过这两个优化后就可以得到常见的KM算法实现，这个算法复杂度是$O(n^4)$，不知道为什么网上都说这个是$O(n^3)$复杂度。如果要达到$O(n^3)$的复杂度，还需要别的优化，比如搜索增广路的时候接着上一次去搜，比如选择更好地计算delta的方式，使得只要调整一步就能引入非饱和点，不过我懒得去细究这个了。Leetcode上有一道最优匹配的题1947. 最大兼容性评分和，可以尝试做一下。
KM算法的矩阵形式
KM算法还有一个矩阵形式，用来求最小权匹配（稍微修改一下也能求最大权匹配），它是对矩阵$W$进行一些处理来求解的，一般是打草稿的时候用这种画矩阵的方法来求。
这个算法基于两个比较显然的定理：
$W$的任意一行或列减去任意一个数$t$，问题的最优解不变，但最优值减少$t$。
如果$W$的每个元素非负，如果存在总报酬为0的分配方法，则该分配方法是最优的。
这个算法通过对$W$进行一些操作，最后得到$n$个不同行不同列的0元素，这样0元素的位置就反映了如何分配。
# 第一步，先让每行每列都出现0元素
# 这个相当于初始化顶点标号，但实现方式是把边权减去min值，权为0的边就属于相等子图
每一行，减去这一行中的最小值 # 每行都有0元素
每一列，减去这一列中的最小值 # 有些列会减去0
# 第二步，继续处理矩阵，让矩阵存在总报酬为0的分配方法
while True:
用最少的线覆盖所有存在0的行和列 # 找到一个最小点覆盖
if 标记了n个行和列: # 已经得到最优匹配
break
else: # 调整顶点标号
找到矩阵中没被划线的最小元素
所有未被标记的行都减去该元素 # 会出现负数
所有被标记的列都加上该元素 # 会把负数加回来
# 第三步，根据 0元素的位置分配任务，略
如何用最少的线数覆盖存在0的行和列，这个其实是求最小点覆盖，打草稿的话，这一步就靠肉眼观察。找到矩阵中没被划线的最小元素是在求$\delta_l$，所有未被标记的行减去该元素是在更新$S$里的顶点标号（等价于给顶点标号加$\delta_l$），所有被标记的列都加上该元素是在更新$T$里面的顶点标号（等价于给顶点标号减$\delta_l$）。总的来说，就是把对标号的调整转换为对矩阵的调整，使得矩阵上0元素对应的边属于相等子图。
KM算法的推导
这里利用线性规划和KKT条件推导出KM算法。求最优匹配是一个线性规划问题，每条边对应一个变量$x_{ij}$：
\[\begin{aligned} \max \sum_{i=1}^n\sum_{j=1}^n &w_{ij}x_{ij}\\ s.t.\ \sum_{i=1}^{n} x_{ij}&=1 & j=1,2,\cdots,n \\ \sum_{j=1}^{n} x_{ij}&=1 & i=1,2,\cdots,n \\ x \ge 0 \end{aligned}\]两个等式约束限制了每个顶点都只能匹配一个顶点。虽然最优匹配这种选边的问题应该被建模成01整数规划问题，但是这里比较特殊，把问题松弛成上面的问题也可以证明最优解的每个$x$取值为0或1，所以两个问题是等价的，解这个松弛了的问题肯定比01规划容易一点，也更好推对偶问题。
该问题的对偶问题是：
\[\begin{aligned} \min \sum_{i=1}^n \alpha_i+\sum_{j=1}^n \beta_j\\ s.t.\ \alpha_i+\beta_j\ge w_{ij} \end{aligned}\]这里$\alpha$和$\beta$是等式约束对应的变量，$\alpha_i+\beta_i\ge w_{ij}$是从$\lambda≥0$得到的约束。如果从更直观的角度来理解，原问题的变量$x_{ij}$对应边，而这里的$\alpha_i$和$\beta_j$分别对应二部图$X$和$Y$里的顶点（KM算法里起了个名字叫顶点标号），我们要在满足约束的情况下最小化顶点对应的标号，约束的限制来源于边权$w_{ij}$。
有了原问题和对偶问题，就可以从KKT条件分析最优解的性质，并且其实只用考虑最核心的互补松弛条件。也就是说$x$以及$\alpha$和$\beta$要满足：
\[(\alpha_i+\beta_j-w_{ij})x_{ij}=0\]即对于最优解而言，要么$w_{ij}=\alpha_i+\beta_j$，要么$x_{ij}=0$。然后利用 原始对偶方法（Primal Dual Method）去设计算法。这里稍微介绍一下这个方法，因为我也不是特别了解，所以可能有些错误。Primal Dual Method是一种给线性规划问题设计算法的方法，在一些问题上（比如这里的最优匹配）可以得到比单纯形法复杂度低很多的算法。它的核心就是利用互补松弛条件，通过给定对偶变量，得到一个关于$x$的、更简单的优化问题，如果求解出的$x$和给定的对偶变量共同满足KKT条件，那就求出最优解了。那么关键的有两点：一是怎么取对偶变量，二是得到一个什么样的更简单的优化问题。
这里先说这个更简单的问题是什么。如果给定$\alpha$和$\beta$，这个问题是（称为限定问题，Restricted Primal Problem，简写为RP问题）：
\[\begin{aligned} \min \sum_{i=1}^n s_i + \sum_{j=1}^n s_{j}\\ s.t.\ \sum_{j=1}^n x_{ij}+s_i=&1, i=1,2,\cdots,n \\ \sum_{i=1}^n x_{ij}+s_j=&1, j=1,2,\cdots,n \\ x_{ij}\ge& 0,\alpha_i+\beta_j= w_{ij} \\ x_{ij}=& 0, \alpha_i+\beta_j\ne w_{ij}\\ s \ge& 0 \end{aligned}\]这里优化变量除了$x$，还引入了额外的变量$s$，$s$相当于单纯形法里的人工变量。如果最优解的$s$都是0，说明我们找到了一个$x$，从约束可以看出它和给定的对偶变量满足互补松弛条件，因而找到了最优解；如果没有找到，就可以调整$\alpha$和$\beta$，得到新的问RP问题进行求解，具体怎么调整的后面说。
再来看看这个优化问题的含义是什么。$\Sigma_j x_{ij}=1$的约束说的是顶点$i$连接的边只能选1个且必须选一个，作为匹配边。而加上人工变量$s_i$意味着不一定需要选择一条连接$i$的匹配边了，也就是不一定要给每个顶点做一个匹配。如果说不加$s_i$约束了是完美匹配，那么加上了$s_i$的约束就变成了求最大匹配。后面两个约束要求在相等子图上求解这个问题，而不是原图上。优化目标指出要尽可能多的让$s_i$更小，也就是说还是尽可能多的要让顶点匹配。如果利用等式约束把变量$s$都消掉会更直观：
\[\begin{aligned} \min\ & 2n-\sum_{i=1}^n\sum_{j=1}^n x_{ij}\\ s.t.\ & x_{ij} \ge 0,\alpha_i+\beta_j= w_{ij} \\ & x_{ij} = 0, \alpha_i+\beta_j\ne w_{ij}\\ & s \ge 0 \end{aligned}\]再变形一下就是：
\[\begin{aligned} \max\ & \sum_{i=1}^n\sum_{j=1}^n x_{ij}\\ s.t.\ & x_{ij} \ge 0, \alpha_i+\beta_j= w_{ij} \\ & x_{ij} = 0, \alpha_i+\beta_j\ne w_{ij}\\ & s \ge 0 \end{aligned}\]这相当于在相等子图上求最大匹配，就样可以理解为什么KM算法要在给定的顶点标号下要在相等子图上求最大匹配了。如果求出来的$x$满足原问题约束，那么问题就解决了，如果没有满足，就需要调整$\alpha$和$\beta$得到新的RP问题去求解。
那具体怎么调整？原始对偶方法给出了一般的线性规划问题（优化目标是最小化）调整对偶变量$y$的方法是：
\[\begin{aligned} \delta &\leftarrow \min \{-\frac{(y^TA_{\cdot j}-c_j)}{(s^TA_{\cdot j})}|s^TA_{\cdot j}>0 \} \\ y &\leftarrow y + \delta s \end{aligned}\]我也不知道怎么推的，看起来跟单纯形法里面的换基有点像。对应到这里的二分图匹配，$A$是关联矩阵，$\alpha$和$\beta$拼到一起是$-y$，$W$是$-c$，所以：
\[\begin{aligned} \delta &\leftarrow \min \{\alpha_i+\beta_j-w_{ij}|\alpha_i+\beta_j > w_{ij}\} \\ \alpha &\leftarrow \alpha+\theta s\\ \beta &\leftarrow \beta+\theta s \end{aligned}\]公式能得到简化是因为$s$的最优解只能取0或1，系数矩阵$A$里面也都是0或者1，而且$s^TA$里面的元素要么是0要么是2，所以第一行公式的分母就可以去掉。。
$s$取1的维度对应存在冲突、找不到匹配的顶点，对应了KM算法里需要调整标号的顶点集合$S$和$T$。$s$只能取1，所以调整的大小取决于$\delta$，这就对应了KM算法里用$\delta_l$来调整标号。
总结
匈牙利算法有2个，分别用来求最大匹配和最优匹配。前者是Berge定理和Hall定理的直接推论，也是最大流算法的特例。后者或者说KM算法利用标号和相等子图把求最优匹配转化为求最大匹配，是利用线性规划的原始对偶方法推出来的一个算法。它还有一个矩阵形式，可以通过在矩阵上划线的方式对比较小的矩阵手动求解，原理涉及到匹配和点覆盖的关系。在实现层面上，2个匈牙利算法都做了一些优化，优化了复杂度，也简化了代码实现。