二分圖帶權匹配（推箱子問題的思考）
轉載請附上原文連結： http://blog.csdn.net/u013351484/article/details/51598270
二分圖帶權匹配（也叫二分圖最優（最佳）匹配，Kuhn－Munkres 演算法）
二分圖帶權匹配，可以把它看作集合 X 的每個頂點到集合 Y 的每個頂點均有邊的二分圖（設權重均為正數，則原來分別在集合 X 與集合 Y 中沒邊的頂點 xi，yi 就可以新增一條權重為 0 的邊，即 w[xi,yi] = 0）；這樣的二分圖必定存在完備匹配。
接下來求解二分圖最大權值匹配過程的所有內容都將以下面這個例子展開：
設有二分圖 G：
最好把這個表畫在紙上 :)
首先要有頂標的概念：
二分圖最大權值匹配通過給集合 X 和集合 Y 的每個頂點定義一個值，稱為頂標；可以設集合 X 裡的頂點 xi 的頂標為 lx[i]，集合 Y 同理。演算法執行過程中確保 lx[i] + ly[j] >= w[xi, yj] 成立（只有這式子成立，才能證明 km 演算法的正確性）。
圖 G 在演算法剛開始執行時要滿足 lx[i] + ly[j] >= w[xi, yj] 的要求可以這樣設定：lx[i] 等於 xi 與集合 Y 中所有頂點連線的最大值；即 lx[i] = max{w[xi, yj]}，j = {0, ... , 4}，ly[i] 全為 0；因此，一開始的圖可以是這樣：
接著要有相等子圖
二分圖中滿足 lx[i] + ly[j] = w[xi,yj] 的邊稱為等邊，而由等邊構成的子圖稱為相等子圖。如上圖就是一個相等子圖；上圖去掉任意邊（至少保留一條）仍為相等子圖；但上圖增加任意邊就不是相等子圖了，因為增加任意邊都會使得新增加的邊的 lx[i] + ly[j] > w[xi, yj]。
接著就會有：
如果二分圖的相等子圖存在完備匹配，那這個完備匹配就是二分圖的最大權值匹配。因為對於二分圖的任意一個匹配，如果它包含於相等子圖，那麼它的邊權和等於所有頂點的頂標和；如果它有的邊不包含於相等子圖，那麼它的邊權和小於所有頂點的頂標和。所以相等子圖的完備匹配一定是二分圖的最大權匹配。看不懂沒關係，把整個過程搞一遍再回來看看也許就理解了。
因此，二分圖最大權值匹配的過程其實就是通過修改集合 X 與集合 Y 的頂標使得二分圖裡 lx[i] + ly[j] = w[xi, yi] 的邊不斷增多，也即相等子圖不斷擴大，直到相等子圖出現完備匹配為止；而此完備匹配就是所求。
km 演算法流程基本基於下圖：
再接著就會有：怎樣修改頂標？什麼情況下要修改頂標？帶著這個疑問往下走~
有了以上的概念（雖然可能還是模模糊糊），我們嘗試找出二分圖 G 的最大權值匹配。
為了理解方便，求解過程都是基於 DFS 的順序。
那~，開始了：
1）
為 X0 找交錯路，X0 -> Y1；這一步是沒有問題的。
2）
為 X1 找交錯路，X1 -> Y1 -> X0；終點不在集合 Y，說明目前的相等子圖不能為 X0，X1（目前遍歷到它們）都找到各自的匹配；接著自然想到需要往當前相等子圖裡面新增新的等邊使匹配成為可能。
哪些邊可以進入相等子圖呢？必須滿足這兩個條件：
1、任何時候必須確保二分圖裡所有頂點滿足 lx[i] + ly[j] >= w[xi, yj]。
2、新加入的邊必須是等邊，因為只有等邊才屬於相等子圖的邊集。
其實在不修改頂標的情況下，條件 1 是肯定滿足的；如果不考慮其他，直接把邊加進去，像這樣：
我們加入了邊 [x1, y2]，確實能使得 X0，X1 找到各自的匹配（X0 -> Y1；X1 -> Y2）但此邊 lx[1] + ly[2] = 5 > w[x1,y2] = 0；因此它並不包含於相等子圖中，這將導致求解錯誤（km 演算法原理就是在相等子圖裡找到完備匹配，若匹配裡含有非相等子圖裡的邊，自然錯誤）。因此，我們必須通過修改 lx[] 或 ly[] 的頂標使得合適的邊加入到相等子圖，使匹配成為可能，同時又能滿足上面兩個條件。
來一些新的概念：
當集合 X 裡的某個點 xi 在相等子圖裡匹配失敗的時候，從 xi 出發的所有交錯路可以組成一顆交錯樹（匈牙利樹）；我們把交錯樹裡屬於 X 集合的所有元素定義為集合 S；屬於 Y 集合的所有元素定義為集合 T。因此，把集合 X 裡不屬於集合 S 的所有元素定義為集合 notS；把集合 Y 裡不屬於集合 T 的所有元素定義為集合 notT。
回到 X1 匹配失敗的時候，有交錯路 X1 -> Y1 ->X0；所以 S = {X0, X1}，T = {Y1}。我們修改頂標的方法是通過把 S 裡的所有頂點的頂標減去一個值 d，且 d 為正數（先別管它怎樣來），把 T 裡的所有頂點的頂標加上 d。這樣做會使得原來相等子圖發生什麼變化？
1、對於在交錯樹（即 xi 屬於 S，yi 屬於 T）的邊，各自的頂標 (lx[i] - d) + (ly[j] + d) = lx[i] + ly[j] 沒有變；也就是這些邊原來在相等子圖裡，修改了頂標後仍在相等子圖裡。如圖中的邊 X0 - Y1 和邊 X1 - Y1。
2、對於 xi 屬於 notS，yi 屬於 notT 的邊，它們各自的頂標並沒有改變，lx[i] + ly[j] 沒有變；也就是這些邊原來不在相等子圖裡，修改了頂標後仍然不在相等子圖裡，即這些邊不會被加入相等子圖。
3、對於 xi 屬於 notS，yi 屬於 T 的邊，各自的頂標 lx[i] + (ly[j] + d) 變大，這說明如果這些邊原來在相等子圖裡的話，修改了頂標後，這些邊會退出相等子圖，比如圖中的邊 Y1 - X2（因為事實是這樣，我是這樣理解的~）。如果這些邊原來就不在相等子圖裡的話，修改了頂標後，就更不可能被加入相等子圖。
4、對於 xi 屬於 S，yi 屬於 notT 的邊，這就不一樣了，這會有 (lx[i] - d) + ly[j]，而修改頂標前為 lx[i] + ly[j] > w[xi,yj]；也即是修改頂標後這些邊各自的頂標會減小！我們令新頂標 lx[i] = lx[i] - d，ly[j] = ly[j]；這意味著這些邊各自新的頂標可能會使得 lx[i] + ly[j] = w[xi,yj] 成立，這表示使等式成立的邊可以加入到相等子圖裡！使匹配成功成為可能。
接下來，d 怎樣來？km 演算法裡是這樣算的：
d = min{lx[i] + ly[j] - w[xi,yj]}，其中 xi 屬於 S，yj 屬於 notT
so why ? 為什麼要這樣算？先跳過，算一遍再說：
d = lx[0] + ly[0] - w[x0,y0] = 4 d = lx[1] + ly[0] - w[x1,y0] = 1
d = lx[0] + ly[2] - w[x0,y2] = 3 d = lx[1] + ly[2] - w[x1,y2] = 5
d = lx[0] + ly[3] - w[x0,y3] = 1 d = lx[1] + ly[3] - w[x1,y3] = 5
d = lx[0] + ly[4] - w[x0,y4] = 7 d = lx[1] + ly[4] - w[x1,y4] = 2
嗯，最小的是 d = 1；修改頂標後的相等子圖為：
可以發現多了兩條邊（X0 - Y3 和 X1 - Y0），同時少了一條邊（X2 - Y1）。
多的兩條邊恰好是上面 8 條式子中滿足 d = 1 的兩條式子的邊；其實也不是那麼巧，是必然！因為在修改頂標的過程中，對於邊 (X0 - Y3)，(lx[0] - d) + ly[3] 把 d = lx[0] + ly[3] - w[x0,y3] 代入可得 (lx[0] - d) + ly[3] = w[x0,y3]，也即是新頂標滿足 lx[0] + ly[3] = w[x0,y3]，因此它是等邊，自然也是相等子圖的邊；邊 (X1 - Y0) 類似。
而少的那條邊是因為修改頂標後 lx[2] + ly[1] = 9 > w[x2,y1] = 8，它不是等邊，因此不是屬於相等子圖裡的邊，把它從相等子圖裡除掉。
d = min{lx[i] + ly[j] - w[xi,yj]}，其中 xi 屬於 S，yj 屬於 notT。
再看一下這個式子，對於上面 8 條式子，如果 d 不取 1，而使 d > 1，比如 d = 2
這會使 (X1 - Y4) 進入相等子圖，(X2 - Y1) 退出相等子圖。看起來可以有 X1 -> Y4 成立，但實際上這不滿足：任何時候必須確保二分圖裡所有頂點滿足 lx[i] + ly[j] >= w[xi, yj] 這個條件，比如 lx[0] + ly[3] = 5 < w[x0,y3] = 6，這將導致 km 演算法錯誤。
而如果 d < 1，比如 d = 0.5，那將不會有新的邊加入相等子圖，因為對於唯一可能加入相等子圖的邊；即xi 屬於 S，yi 屬於 notT 的邊都將有 (lx[i] - d) + ly[j] > w[xi,yj]（因為 d = 1 時才使等式成立），因此不可能有等邊產生。
所以有 d = min{lx[i] + ly[j] - w[xi,yj]}，其中 xi 屬於 S，yj 屬於 notT。
3）
成功為 X1 找到匹配後（X1 -> Y0），接著為 X2 找交錯路，也失敗了，它的交錯樹只有 X2 頂點本身；S = {X2}，T = {}。
d = lx[2] + ly[0] - w[x2,y0] = 6
d = lx[2] + ly[1] - w[x2,y1] = 1
d = lx[2] + ly[2] - w[x2,y2] = 1
d = lx[2] + ly[3] - w[x2,y3] = 5
d = lx[2] + ly[4] - w[x2,y4] = 3
可見 d = 1
下圖（左）可以看到，有 (X2 - Y1 和 X2 - Y2) 兩條邊加入相等子圖，沒有邊退出相等子圖；需要改變的頂標只有 lx[2]。
因而有交錯路 X2 -> Y1 -> X0 -> Y3，對路徑屬性取反即可得到上圖（右）的匹配。
4）
接著為 X3 找交錯路，還好，一次成功 :)
5）
最後為 X4 找交錯路，交錯樹為 X4 -> Y2 -> X3；S = {X4, X3}，T = {Y2}。
d = lx[4] + ly[0] - w[x4,y0] = 8 d = lx[3] + ly[0] - w[x3,y0] = 3
d = lx[4] + ly[1] - w[x4,y1] = 3 d = lx[3] + ly[1] - w[x3,y1] = 3
d = lx[4] + ly[3] - w[x4,y3] = 1 d = lx[3] + ly[3] - w[x3,y3] = 1
d = lx[4] + ly[4] - w[x4,y4] = 9 d = lx[3] + ly[4] - w[x3,y4] = 4
d = 1
下圖可以看到，有 (X4 - Y3 和 X3 - Y3) 兩條邊加入相等子圖，邊 (X2 - Y2) 退出相等子圖。
因而有交錯路 X4 -> Y3 -> X0 ->Y1 -> X2，gg~，終點不屬於集合 Y；不慌 :)，還有另一條 X4 -> Y2 -> X3 -> Y3 -> X0 -> Y1 -> X2；完蛋，兩條路都沒有。只能按步驟再做一遍；注意這顆以 X4 為起點的交錯樹有兩個分支，要把兩條路的都算進來，即 S = {X0, X2, X3, X4}，T = {Y1, Y2, Y3};
d = lx[0] + ly[0] - w[x0,y0] = 3 d = lx[3] + ly[0] - w[x3,y0] = 2
d = lx[0] + ly[4] - w[x0,y4] = 6 d = lx[3] + ly[4] - w[x3,y4] = 3
d = lx[2] + ly[0] - w[x2,y0] = 5 d = lx[4] + ly[0] - w[x4,y0] = 7
d = lx[2] + ly[4] - w[x2,y4] = 2 d = lx[4] + ly[4] - w[x4,y4] = 8
d = 2
下圖可以看到，有 (X2 - Y4 和 X3 - Y0) 兩條邊加入相等子圖，邊 (X1, Y1) 退出相等子圖。
在上面（左）圖中再來為 X4 找交錯路 X4 -> Y2 -> X3 -> Y0 -> X1，沒找到？心驚~；還有一條 X4 -> Y2 -> X3 -> Y3 -> X0 -> Y1 -> X2 -> Y4，看到了可愛的終點！路徑屬性取反就變成了上面（右）圖。至此，集合 X 裡的所有頂點都已經有對應的匹配，也就是完備匹配！也即此二分圖的最大權值匹配！
X0 -> Y1
X1 -> Y0
X2 -> Y4
X3 -> Y3
X4 -> Y2
最大權值為 30
那要求最小權值匹配怎麼辦？很簡單，在求解前把所有權值取相反數，求得的結果後再對結果取相反數就可以了。
接著談一下複雜度，比較權威的說法是這樣（如果理解了以上講的內容，複雜度就好理解）：
普通做法時間複雜度為 O(n^4) —— 需要找 O(n) 次增廣路，每次增廣最多需要修改 O(n) 次頂標，每次修改頂標時由於要列舉邊來求 d 值，複雜度為 O(n^2)。實際上 Kuhn-Munkers 演算法的複雜度是可以做到 O(n^3) 的。我們給每個 Y 頂點一個“鬆弛量”函式 slack ，每次開始找增廣路時初始化為無窮大。在尋找增廣路的過程中，檢查邊 (xi,yj) 時，如果它不在相等子圖中，則讓 slack[j] 變成原值與 lx[i]
+ ly[j] - w[xi,yj] 的較小值。這樣，在修改頂標時，取所有不在交錯樹中的 Y 頂點的 slack 值中的最小值作為 d 值即可。但還要注意一點：修改頂標後，要把所有的 slack 值都減去 d。
同樣，如果理解了以上講的內容，程式碼也能輕鬆理解。
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define N 110
const int inf = (1<<30);
int m[N][N];
int lx[N];
int ly[N];
int vx[N];
int vy[N];
int slack[N];
int link[N];
int m_count = 5;
int h_count = 5;
int row, col;
void cal_dis()
{
int i, j;
m[0][0] = 3;
m[0][1] = 7;
m[0][2] = 4;
m[0][3] = 6;
m[0][4] = 0;
m[1][0] = 4;
m[1][1] = 5;
m[1][2] = 0;
m[1][3] = 0;
m[1][4] = 3;
m[2][0] = 2;
m[2][1] = 8;
m[2][2] = 7;
m[2][3] = 3;
m[2][4] = 5;
m[3][0] = 3;
m[3][1] = 4;
m[3][2] = 6;
m[3][3] = 5;
m[3][4] = 2;
m[4][0] = 1;
m[4][1] = 7;
m[4][2] = 9;
m[4][3] = 8;
m[4][4] = 0;
/* 初始化頂標 */
for(i = 0; i < N; i++)
lx[i] = -inf;
for(i = 0; i < m_count; i++){
for(j = 0; j < h_count; j++){
if(m[i][j] > lx[i]){
lx[i] = m[i][j];
}
}
}
memset(ly, 0, sizeof(ly));
}
int dfs(int u)
{
int v;
vx[u] = 1;
for(v = 0; v < h_count; v++){
if(!vy[v]){
int t = lx[u] + ly[v] - m[u][v];
if(!t){
vy[v] = 1;
if(link[v] == -1 || dfs(link[v])){
link[v] = u;
return 1;
}
}
else
if(t < slack[v])
slack[v] = t;
}
}
return 0;
}
int km()
{
int i, j;
int ans = 0;
memset(link, -1, sizeof(link));
for(i = 0; i < m_count; i++){
while(1){
int d = inf;
for(j = 0; j < h_count; j++){
slack[j] = inf;
}
memset(vx, 0, sizeof(vx));
memset(vy, 0, sizeof(vy));
if(dfs(i))
break;
for(j = 0; j < h_count; j++){
if(!vy[j] && slack[j] < d)
d = slack[j];
}
for(j = 0; j < h_count; j++){
if(vx[j])
lx[j] -= d;
}
for(j = 0; j < h_count; j++){
if(vy[j])
ly[j] += d;
}
}
}
for(i = 0; i < h_count; i++)
ans += m[link[i]][i];
return ans;
}
int main()
{
int i;
int ans;
cal_dis();
ans = km();
printf("ans: %d\n", ans);
for(i = 0; i < 5; i++){
if(link[i] != -1){
printf("x%d -> y%d\n", link[i], i);
}
}
return 0;
}
因為程式使用 DFS 求解，可見程式得到的結果跟我們的推算是一致的。
心得：
km 演算法我在去年也理解過了一遍，但當時忙於刷 OJ，那時的急切心態使自己沒有好好做總結。時隔一年，最近要為推箱子求解程式的 IDA* 演算法設計一個比較準確的估值函式，很明顯二分圖最小權值匹配就能活脫脫的運用上去。然而 km 演算法幾乎全忘~；咬牙在網上翻了 2 天資料，成功運用在程式中後，取得了比較滿意的效果；雖然對 km 演算法還有一些不理解的地方，而且可能有理解不對的地方，以後有時間再補全。迅速把自己理解的過程記錄下來。其實很早就有這樣的意識，但很多時候都沒有這樣去做：對自己剛掌握的知識，如果認為這些知識對自己很有用，同時又不容易理解；那最好把自己理解它們的過程記錄下來，方便以後需要的時候查閱，不然，就算一次弄明白了，往後需要就得重新在網上一陣惡啃（說真的，網上覆制貼上黨還是很多的，搜尋到自己需要的東西不容易）~，那樣反而更耗時間，為何不稍微花點心血趁熱打鐵呢？
儘管 n 為箱子數，但 O(n^3) 的複雜度還是使得求解的很多時間用在二分圖匹配上面~，原因是用到 km 演算法的 H() 函式和 dfs() 函式呼叫太頻繁，加起來都 4,000,000+ 次，這只是一個很小的關卡~；如果 n^3 那個 3 能減少會非常地不錯。同時發現了 Linux 下為 gcc/g++ 編譯器提供效能評測的 gprof 命令 及其類似命令，說真的，非常地有用。
