匈牙利算法

<image>

匈牙利算法是一种在多项式时间内求解任务分配问题的组合优化算法，并推动了后来的原始对偶方法。美国数学家哈罗德·W·库恩于1955年提出该算法。此算法之所以被称作匈牙利算法，是因为算法很大一部分是基于以前匈牙利数学家科尼格·德内什和艾盖瓦里·耶内的工作之上创建起来的。

詹姆士·芒克勒斯在1957年回顾了该算法，并发现它的时间复杂度为（强）多项式时间。此后该算法被称为库恩-芒克勒斯算法或芒克勒斯分配算法。原始算法的时间复杂度为$O(n^4)$，但杰克·爱德蒙斯与理查德·卡普发现可以修改算法达到$O(n^3)$运行时间。小莱斯特·伦道夫·福特和德尔伯特·雷·富尔克森将该方法推广到一般运输问题，称作福特-富尔克森算法。2006年发现卡尔·雅可比在19世纪就解决了指派问题，该解法在他死后在1890年以拉丁文发表。

## 问题

你有三个工人：吉姆，史提夫和艾伦。

你需要其中一个清洁浴室，另一个打扫地板，第三个洗窗，但他们每个人对各项任务要求不同数目数量的钱。

以最低成本的分配工作的方式是什么？

可以用工人做工的成本矩阵来表示该问题。例如：

清洁浴室｜打扫地板｜洗窗

吉姆｜$2｜$3｜$3

史提夫｜$3｜$2｜$3

艾伦｜$3｜$3｜$2

当把匈牙利方法应用于上面的表格时，会给出最低成本：为6美元，让吉姆清洁浴室、史提夫打扫地板、艾伦清洗窗户就可以达到这一结果。

## 设定

给定一个 $n\times n$ 的非负矩阵，其中第 $i$ 行第 $j$ 列元素表示把第 $i$ 个工人派到第 $j$个工作的成本。我们必须找到成本最低的工人工作分配。如果目标是找到最高成本的分配，该问题可以将每个成本都换为最高一个成本减去该成本以适应题目。

如果用二分图来阐述该问题可以更容易描述这个算法。对于一个有 $n$个工人节点（$S$）与 $n$个工作节点（$T$）的完全二分图 $G=(S\cup T, E)$，每条边都有 $c(i, j)$ 的非负成本。我们要找到最低成本的完美匹配。

如果函数 $y: (S \cup T) \mapsto \mathbb{R}$ 满足对于每个 $i \in S, j \in T$ 都有 $y(i)+y(j) \leq c(i, j)$，则把该函数叫做势。势 $y$ 的值为 $\sum_{v\in S\cup T} y(v)$。可以看出，每个完美匹配的成本最低是每个势的值。匈牙利算法找出了完美匹配及与之成本/值相等的势，这证明了两者的最优性。实际上它找到了紧边集的完美匹配：紧边 $ij$ 是指对于势 $y$ 满足 $y(i)+y(j) = c(i, j)$。我们将紧边子图表示为 $G_y$。$G_y$ 中的完美匹配的成本（如果存在）就等于 $y$的值。

## 用二分图描述此算法

在算法中，我们维护一个势函数 $y$、其对应的紧边子图 $G_y$、当前匹配 $M$，以及与 $M$ 相容的 $G_y$ 的一个定向（记为 $\overrightarrow{G_y}$）。这里约定：属于当前匹配 $M$ 的紧边定向为从 $T$ 到 $S$，其余紧边定向为从 $S$ 到 $T$。初始时，$y$ 处处为 0。

所有边都由 $S$ 指向 $T$（因此 $M$ 为空）。每一步中，我们或者改变 $y$ 使其值增加。

或者改变边的方向以增大匹配 $M$。我们始终保持 $M$ 的所有边都是紧边。当 $M$ 为完美匹配时结束。

在一般的步骤中，令 $R_S \subseteq S$ 和 $R_T \subseteq T$ 为 $M$ 未覆盖的节点（则 $R_S$ 包含 $S$ 中入度为零的节点，而 $R_T$ 中包含 $T$ 中出度为零的节点）。

令 $Z$ 为在 $\overrightarrow{G_y}$ 中从 $R_S$ 出发，只沿有向紧边可达的节点集。可由广度优先搜索求得。

若 $R_T \cap Z$ 非空，则将 $\overrightarrow{G_y}$ 中从 $R_S$ 到 $R_T$ 的有向路径反向。则相应匹配数增加1。

若 $R_T \cap Z$ 为空，则令 $\Delta := \min \{c(i,j)-y(i)-y(j): i \in Z \cap S, j \in T \setminus Z\}$。$\Delta$ 为正，因为 $Z \cap S$ 与 $T \setminus Z$ 之间没有紧边。

在 $Z \cap T$ 中的节点将$y$增加$\Delta$ 并在 $Z \cap S$ 中节点将 $y$减小 $\Delta$。

得到的$y$仍然是势。图 $G_y$ 改变了，但它仍包含$M$。我们把新的边从$S$指向$T$。

由 $\Delta$ 的定义，$R_S$ 可达的节点集 $Z$ 增大（注意到紧边的数量不一定增加）。

我们重复这些步骤直到$M$为完美匹配，该情形下给出的是最小成本（即时间消耗）的匹配。此版本的运行时间为 $O(n^4)$：$M$ 增广 $n$次。

在 $M$ 不改变的一个阶段中，势最多改变 $n$ 次（因为 $Z$ 每次都增加）。改变势所需的时间在 $O(n^2)$。

### 证明：改变势函数y不改变M

证明改变$y$后$M$中每条边不发生改变，这等价于证明对于$M$中任意边，它的两个顶点要么都在$Z$中，要么都不在$Z$中。为此，定义$vu$为$M$中一条从$T$到$S$的边，则若$v\in Z$，那么有$u\in Z$。（反证）假设$u \in Z$但$v \notin Z$；由于$u$是一个匹配边的末端点，$u\notin R_S$，因此存在从$R_S$到$u$的有向路径。这条路径必不能经过$v$（根据假设），因此这条路径上紧邻$u$的点是其他点$v' \in T$。$v'u$是一条从$T$到$S$的紧边因此也是$M$中的一个元素。但因此$M$包含两条有共点的边，与$M$是匹配的定义矛盾。因此$M$中每条边的两个顶点要么都在$Z$中，要么都不在$Z$中。

### 证明：y仍然是势函数

证明$y$在更改之后是势函数，等价于证明不存在总势超过成本的边。这已经在之前的论述中已经为$M$中的边建立这个概念，因此我们考虑任意从$S$到$T$的边$uv$。如果$y(u)$升高了$\Delta$，那么：(1)$v \in Z \cap T$，在这种情况下，$y(v)$减小了$\Delta$，使总体的势未发生改变；(2)$v \in T \setminus Z$，这种情况下$\Delta$的定义保证了$y(u)+y(v)+\Delta \leq c(u,v)$。因此$y$仍然是势函数。

## 矩阵解释

给定 $n$ 个工人和任务，以及一个包含分配给每个工人一个任务的成本的 $n\times n$矩阵，寻找成本最小化分配。

首先把问题写成下面的矩阵形式

$$\begin{bmatrix} a_1 & a_2 & a_3 & a_4\\ b_1 & b_2 & b_3 & b_4\\ c_1 & c_2 & c_3 & c_4\\ d_1 & d_2 & d_3 & d_4\end{bmatrix}$$

其中 $a,b,c,d$是执行任务 1, 2, 3, 4 的工人。$a_1, a_2, a_3, a_4$ 分别表示当工人 $a$ 做任务 1, 2, 3, 4 时的时间损失（成本）。对于其他符号也同样适用。该矩阵是方阵，所以每个工人只能执行一个任务。

第 1 步

接下来我们对矩阵的行进行操作。将所有 $a_i$ （$i$ 从 1 到 4）中最小的元素取走，并将该行每个元素都减去刚刚取走的元素。这会让该行至少出现一个零（当一行有两个相等的最小元素时会得到多个零）。对此过程为所有行重复。我们现在得到一个每行至少有一个零的矩阵。现在我们尝试给工人指派任务，以使每个工人只做一项任务，并且每个情况的耗散都为零。说明如下。

$0$｜$a_2'$｜$a_3'$｜$a_4'$

$b_1'$｜$b_2'$｜$b_3'$｜$0$

$c_1'$｜$0$｜$c_3'$｜$c_4'$

$d_1'$｜$d_2'$｜$0$｜$d_4'$

每一行的 0 的个数可以超过 1。

若有出现行中 0 的个数超过 1，在最坏的情况下，会需要在 $n!$ 的组合中找到一个成本最小的配对。

第 2 步

有时此阶段的该矩阵不能符合指派的要求，例如下面所示矩阵。

$0$｜$a_2'$｜$a_3'$｜$a_4'$

$b_1'$｜$b_2'$｜$b_3'$｜$0$

$0$｜$c_2'$｜$c_3'$｜$c_4'$

$d_1'$｜$0$｜$d_3'$｜$d_4'$

在上述情形下，不能做出指派。注意到任务 1 由工人 $a$ 和 $c$ 做都很高效 ($a_1'$ 和 $c_1'$ 为 0 )。但两个工人无法分配到同一个任务。

还注意到，没有任何一个工人能有效地做任务 3 ($a_3'$、$b_3'$、$c_3'$、$d_3'$ 皆不为 0)。

为了克服这个问题，我们对所有列重复上述流程（即每一列所有元素都减去该列最小元素）并检查是否可以完成分配。

大多数情况下，这都会给出结果，但如果仍然是不可行，那么我们需要继续下去。

第 3 步

必须用尽可能少的列或行标记来覆盖矩阵中的所有零。下面的过程是完成这个要求的一种方法：

首先，尽可能多地分配任务，用 0' 表示的零为已指派的任务。

- 第 1 行有一个零，所以用 0' 表示为分配。第 3 行的 0 由于处于同一列而被划掉。
- 第 2 行有一个零，所以用 0' 表示为分配。
- 第 3 行只有一个已经划掉的零，所以不能分配。
- 第 4 行有两个未划掉的零。可以分配任何一个（都是最优），并将另一个零划去。

(在此阶断任何合法的分配都可以接受，因此若一开始分配的是第 3 行的 0，就会使第 1 行的 0 被划掉。)

$0'$｜$a_2'$｜$a_3'$｜$a_4'$

$b_1'$｜$b_2'$｜$b_3'$｜$0'$

$0$｜$c_2'$｜$c_3'$｜$c_4'$

$d_1'$｜$0'$｜$0$｜$d_4'$

现在到了画图的部分。

- 标记所有未分配的行（第 3 行）。
- 标记所有新标记的行中 0所在（且未标记）的对应列（第 1 列）。
- 标记所有在新标记的列中有分配的行（第 1 行）。
- 对所有未分配的行重复上述过程。

×

$0'$｜$a_2'$｜$a_3'$｜$a_4'$｜×

$b_1'$｜$b_2'$｜$b_3'$｜$0'$

$0$｜$c_2'$｜$c_3'$｜$c_4'$｜×

$d_1'$｜$0'$｜$0$｜$d_4'$

现在划掉所有已标记的列和未标记的行（第 1 列和第 2, 4 行）。

×

$0'$｜$a_2'$｜$a_3'$｜$a_4'$｜×

$b_1'$｜$b_2'$｜$b_3'$｜$0'$

$0$｜$c_2'$｜$c_3'$｜$c_4'$｜×

$d_1'$｜$0'$｜$0$｜$d_4'$

上述详细的描述只是画出覆盖所有 0 的（直、行）线的一种方法。也可以使用其他方法。

第 4 步

现在删除已画线的行和列。这将留下一个矩阵如下：

$$\begin{bmatrix} a_2 & a_3 & a_4\\ c_2 & c_3 & c_4\end{bmatrix}$$

返回到步骤 1，然后重复这个过程，直到矩阵是空的；当矩阵为空时，意即画出的覆盖所有 0 的（直、行）线的个数等于原本矩阵的行数(或列数) $n$，在上述的例子中 $n = 4$，此时代表有一个分配的最佳解存在于这些 0 的组合之中，我们可以在步骤 1 中从所有分配的组合中找到成本最小的解。

## 外部链接

- Bruff, Derek, "The Assignment Problem and the Hungarian Method",
- Mordecai J. Golin, Bipartite Matching and the Hungarian Method, Course Notes, Hong Kong University of Science and Technology.
- R. A. Pilgrim, Munkres' Assignment Algorithm. Modified for Rectangular Matrices , Course notes, Murray State University.
- Mike Dawes, The Optimal Assignment Problem, Course notes, University of Western Ontario.
- On Kuhn's Hungarian Method – A tribute from Hungary , András Frank, Egervary Research Group, Pazmany P. setany 1/C, H1117, Budapest, Hungary.
- Lecture: Fundamentals of Operations Research - Assignment Problem - Hungarian Algorithm , Prof. G. Srinivasan, Department of Management Studies, IIT Madras.
- Extension: Assignment sensitivity analysis (with O(n^4) time complexity) , Liu, Shell.
- Solve any Assignment Problem online , provides a step by step explanation of the Hungarian Algorithm.

### 实现

（请注意，并非所有这些都满足 $O(n^3)$ 时间约束。）

- C implementation with $O(n^3)$ time complexity
- Java implementation of $O(n^3)$ time variant
- Python implementation (see also here)
- Ruby implementation with unit tests
- C# implementation
- D implementation with unit tests (port of the Java $O(n^3)$ version)
- Online interactive implementation Please note that this implements a variant of the algorithm as described above.
- Graphical implementation with options (Java applet)
- Serial and parallel implementations.
- Implementation in Matlab and C
- Perl implementation
- Lisp implementation
- C++ (STL) implementation (multi-functional bipartite graph version)
- C++ implementation
- C++ implementation of the $O(n^3)$ algorithm (BSD style open source licensed)
- Another C++ implementation with unit tests
- Another Java implementation with JUnit tests (Apache 2.0)
- MATLAB implementation
- C implementation
- Javascript implementation
- The clue R package proposes an implementation, solve_LSAP