Loading [Contrib]/a11y/accessibility-menu.js
$\newcommand{\O}{\mathrm{O}}$ My Algorithm : kopricky アルゴリズムライブラリ

kopricky アルゴリズムライブラリ

Hungarian

コードについての説明

割当問題を効率よく解くアルゴリズム. 割当問題とは $n$ 人の人と $n$ 個の仕事があり, また各(人 $i$, 仕事 $j$) のペアについて人 $i$ が仕事 $j$ を行うのにかかる時間が与えられているとする. このとき $n$ 人の人をちょうど $1$ つずつ別々の仕事に割り当てることで仕事にかかるのべ時間を最小にしたいという問題のことを言う.
二部グラフのマッチングの重みつきの場合のアルゴリズムである. 最小費用流のアルゴリズムでも解くことができるが, 計算量としては $\O (n m \log n)$ かかる (正確には二分ヒープによる Dijkstra 法を用いて最短経路に沿ってフローを増加させるアルゴリズムの場合の話で フィボナッチヒープによる Dijkstra 法 を用いると $\O (n m + n^2 \log n)$ で可能になる). 一方 Hungarian 法では $\O (n^3)$ の計算量で解くことが可能である. また学問的に割当問題は重み付き(線形)マトロイド交差問題というより一般的な枠組みで捉えられることが多い.
ちなみに Hungarian はこの問題に対する重要な貢献をした $2$ 人のハンガリー人数学者 Kőnig, Egerváry にちなんで呼ばれている.
Hungarian 法の説明は有名な Goemans 先生の Lecture note が分かりやすかった. こちらにのっているのは $\O (n^4)$ までの話だが, ポテンシャルの更新の際にマッチングが増えない場合は集合 $L$ の更新は新しく $w_{ij} = 0$ となった辺による影響のみを追加で考えればよく、 結局のところマッチングが $1$ 増えるまでに見る辺の本数を $\O (n^2)$ にできる($\delta$ を求める部分も良い感じにやる). よって全体で $\O (n^3)$ のアルゴリズムとなる.
初めに双対問題の実行可能解から始めて相補性条件から導かれるグラフ上で増加道アルゴリズムを用いてマッチングの更新を行い, 完全マッチングが得られるまで続けるといった Primal-Dual 法を行う (この手順で完全マッチングが得られたとき主問題, 双対問題の実行可能解に入っていること, 相補性条件が成り立つことから最適性が保証される).
以上の操作を辺コストを要素とする行列上で行っている実装が多い気がするが、個人的に理解しづらかったので陽にグラフを構築する形で実装した. 速度的には遅くはないが、 高速化したいなら Jonker-Volgenant などを実装すると良い.
(追記) 人の数 $n$ が仕事の数 $m$ より大きい場合にも気をつけて実装することで $\O(n m^2)$ のアルゴリズムを実装することができたので以下はその実装となっている.

(関数)
solve$(cost)$ : $n \times m(n \geq m)$ の辺コストを要素とする配列を渡し, (総コスト, 割当) を返す.
$n > m$ の場合はちょうど $m$ 組の割り当てを構成し、 どの要素にも割り当てられないものは $-1$ が格納される.

時間計算量: $\O (n m^2)$

コード

  1. template<typename T> class Hungarian
  2. {
  3. private:
  4.     const int U, V;
  5.     vector<vector<int> > graph;
  6.     vector<T> dual;
  7.     vector<int> alloc, rev_alloc, prev;
  8.     const vector<vector<T> >& cost;
  9.     int matching_size;
  10.     T diff(const int i, const int j){
  11.         return cost[i][j] - dual[i] - dual[U + j];
  12.     }
  13.     void init_feasible_dual(){
  14.         for(int i = 0; i < U; ++i){
  15.             dual[i] = 0;
  16.             for(int j = 0; j < V; ++j){
  17.                 dual[U + j] = min(dual[U + j], cost[i][j]);
  18.             }
  19.         }
  20.     }
  21.     void construct_graph(){
  22.         for(int i = 0; i < U; ++i){
  23.             for(int j = 0; j < V; ++j){
  24.                 graph[i][j] = (diff(i, j) == 0 && rev_alloc[j] != i);
  25.             }
  26.         }
  27.     }
  28.     bool find_augmenting_path(const int cur, const int prv, int& pos){
  29.         prev[cur] = prv;
  30.         if(cur >= U){
  31.             if(rev_alloc[cur - U] < 0) return true;
  32.             if(find_augmenting_path(rev_alloc[cur - U], cur, pos)){
  33.                 graph[rev_alloc[cur - U]][cur - U] = 1;
  34.                 return true;
  35.             }
  36.         }else{
  37.             const int MX = (alloc[cur] < 0 && pos == U) ? U : V;
  38.             for(int i = 0; i < MX; ++i){
  39.                 if(graph[cur][i] && prev[U + i] < 0 && find_augmenting_path(U + i, cur, pos)){
  40.                     graph[cur][i] = 0, alloc[cur] = i, rev_alloc[i] = cur;
  41.                     return true;
  42.                 }
  43.             }
  44.             if(alloc[cur] < 0 && pos < U){
  45.                 graph[cur][pos] = 0, alloc[cur] = pos, rev_alloc[pos] = cur, prev[U + pos] = cur;
  46.                 return ++pos, true;
  47.             }
  48.         }
  49.         return false;
  50.     }
  51.     void update_dual(const T delta){
  52.         for(int i = 0; i < U; ++i) if(prev[i] >= 0) dual[i] += delta;
  53.         for(int i = U; i < U + V; ++i) if(prev[i] >= 0) dual[i] -= delta;
  54.     }
  55.     void maximum_matching(bool initial=false){
  56.         int pos = initial ? V : U;
  57.         for(bool update = false;; update = false){
  58.             fill(prev.begin(), prev.end(), -1);
  59.             for(int i = 0; i < U; ++i){
  60.                 if(alloc[i] < 0 && find_augmenting_path(i, 2 * U, pos)){
  61.                     update = true, ++matching_size;
  62.                     break;
  63.                 }
  64.             }
  65.             if(!update) break;
  66.         }
  67.     }
  68.     int dfs(const int cur, const int prv, vector<int>& new_ver){
  69.         prev[cur] = prv;
  70.         if(cur >= U){
  71.             if(rev_alloc[cur - U] < 0) return cur;
  72.             else return dfs(rev_alloc[cur - U], cur, new_ver);
  73.         }else{
  74.             new_ver.push_back(cur);
  75.             for(int i = 0; i < V; ++i){
  76.                 if(graph[cur][i] && prev[U + i] < 0){
  77.                     const int res = dfs(U + i, cur, new_ver);
  78.                     if(res >= U) return res;
  79.                 }
  80.             }
  81.         }
  82.         return -1;
  83.     }
  84.     int increase_matching(const vector<pair<int, int> >& vec, vector<int>& new_ver){
  85.         for(const auto& e : vec){
  86.             if(prev[e.first] < 0){
  87.                 const int res = dfs(e.first, e.second, new_ver);
  88.                 if(res >= U) return res;
  89.             }
  90.         }
  91.         return -1;
  92.     }
  93.     void hint_increment(int cur){
  94.         while(prev[cur] != 2 * U){
  95.             if(cur >= U){
  96.                 graph[prev[cur]][cur - U] = 0, alloc[prev[cur]] = cur - U, rev_alloc[cur - U] = prev[cur];
  97.             }else{
  98.                 graph[cur][prev[cur] - U] = 1;
  99.             }
  100.             cur = prev[cur];
  101.         }
  102.     }
  103. public:
  104.     Hungarian(const vector<vector<T> >& _cost)
  105.          : U((int)_cost.size()), V((int)_cost[0].size()), graph(U, vector<int>(U, 1)), dual(U + V, numeric_limits<T>::max()),
  106.             alloc(U, -1), rev_alloc(U, -1), prev(2 * U), cost{_cost}, matching_size(0){
  107.         assert(U >= V);
  108.     }
  109.     pair<T, vector<int> > solve(){
  110.         init_feasible_dual(), construct_graph();
  111.         bool end = false;
  112.         maximum_matching(true);
  113.         while(matching_size < U){
  114.             vector<pair<T, int> > cand(V, {numeric_limits<T>::max(), numeric_limits<int>::max()});
  115.             for(int i = 0; i < U; ++i){
  116.                 if(prev[i] < 0) continue;
  117.                 for(int j = 0; j < V; ++j){
  118.                     if(prev[U + j] >= 0) continue;
  119.                     cand[j] = min(cand[j], {diff(i, j), i});
  120.                 }
  121.             }
  122.             while(true){
  123.                 T delta = numeric_limits<T>::max();
  124.                 for(int i = 0; i < V; ++i){
  125.                     if(prev[U + i] >= 0) continue;
  126.                     delta = min(delta, cand[i].first);
  127.                 }
  128.                 update_dual(delta);
  129.                 vector<pair<int, int> > vec;
  130.                 vector<int> new_ver;
  131.                 for(int i = 0; i < V; ++i){
  132.                     if(prev[U + i] >= 0) continue;
  133.                     if((cand[i].first -= delta) == 0) vec.emplace_back(U + i, cand[i].second);
  134.                 }
  135.                 int res = increase_matching(vec, new_ver);
  136.                 if(res >= U){
  137.                     hint_increment(res);
  138.                     if(++matching_size == U) end = true;
  139.                     else construct_graph();
  140.                     break;
  141.                 }else{
  142.                     for(const int v : new_ver){
  143.                         for(int i = 0; i < V; ++i){
  144.                             if(prev[U + i] >= 0) continue;
  145.                             cand[i] = min(cand[i], {diff(v, i), v});
  146.                         }
  147.                     }
  148.                 }
  149.             }
  150.             if(!end) maximum_matching();
  151.         }
  152.         T total_cost = 0;
  153.         for(int i = 0; i < U; ++i){
  154.             if(alloc[i] < V) total_cost += cost[i][alloc[i]];
  155.             else alloc[i] = -1;
  156.         }
  157.         return make_pair(total_cost, alloc);
  158.     }
  159. };

verify 用の問題

AOJ : Bipartite Matching(重みなし 2 部マッチング) 提出コード
AOJ : marukaite 提出コード
yosupo さんの library checker : Assignment Problem 提出コード