匈牙利算法

匈牙利算法是一種在多項式時間內求解任務分配問題的組合優化演算法，並推動了後來的原始對偶方法。美國數學家哈羅德·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