Content deleted Content added
m Removed unnecessary conditional "only following edges that are tight". |
|||
| (34 intermediate revisions by 25 users not shown) | |||
Line 1:
{{Short description|Polynomial-time algorithm for the assignment problem}}
The '''Hungarian method''' is a [[combinatorial optimization]] [[algorithm]] that solves the [[assignment problem]] in [[polynomial time]] and which anticipated later [[Duality (optimization)|primal–dual methods]]. It was developed and published in 1955 by [[Harold Kuhn]], who gave it the name "Hungarian method" because the algorithm was largely based on the earlier works of two
[[James Munkres]] reviewed the algorithm in 1957 and observed that it is [[
==The problem==
Line 8:
=== Example ===
In this simple example, there are three workers: Alice, Bob and
:{| class="wikitable" style="text-align:center;"
! {{diagonal split header|Worker |Task}}
Line 25:
| '''$3'''
|-
!
| $9
| '''$4'''
Line 31:
|}
The Hungarian method, when applied to the above table, would give the minimum cost: this is $15, achieved by having Alice clean the bathroom,
:{| class="wikitable" style="text-align:center;"
! {{diagonal split header|Sweep|Clean}}
! Alice
! Bob
!
|-
! Alice
Line 48:
| $18
|-
!
| '''$15'''
| $16
Line 55:
=== Matrix formulation ===
In the matrix formulation, we are given
This can be expressed as permuting the rows of a cost matrix ''C'' to minimize the [[Trace (linear algebra)|trace]] of a matrix,
Line 68:
=== Bipartite graph formulation ===
The algorithm can equivalently be described by formulating the problem using a bipartite graph. We have a [[complete bipartite graph]] <math>G=(S, T; E)</math> with {{mvar|n}} worker vertices ({{mvar|S}}) and {{mvar|n}} job vertices ({{mvar|T}}), and the edges ({{mvar|E}}) each have a
==The algorithm in terms of bipartite graphs==
Let us call a function <math>y: (S\cup T) \to \mathbb{R}</math> a '''potential''' if <math>y(i)+y(j) \leq c(i, j)</math> for each <math>i \in S, j \in T</math>.
The ''value'' of potential {{mvar|y}} is the sum of the potential over all vertices:
The cost of each perfect matching is at least the value of each potential: the total cost of the matching is the sum of costs of all edges; the cost of each edge is at least the sum of potentials of its endpoints; since the matching is perfect, each vertex is an endpoint of exactly one edge; hence the total cost is at least the total potential.▼
:<math>\sum_{v\in S\cup T} y(v)</math>.
The Hungarian method finds a perfect matching and a potential such that the matching cost equals the potential value. This proves that both of them are optimal. In fact, the Hungarian method finds a perfect matching of '''tight edges''': an edge <math>ij</math> is called tight for a potential {{mvar|y}} if <math>y(i)+y(j) = c(i, j)</math>. Let us denote the [[Glossary of graph theory#Subgraphs|subgraph]] of tight edges by <math>G_y</math>. The cost of a perfect matching in <math>G_y</math> (if there is one) equals the value of {{mvar|y}}.▼
▲The cost of each perfect matching is at least the value of each potential: the total cost of the matching is the sum of costs of all edges it contains; the cost of each edge is at least the sum of potentials of its endpoints; since the matching is perfect, each vertex is an endpoint of exactly one edge; hence the total cost is at least the total potential.
▲The Hungarian method finds a perfect matching and a potential such that the matching cost equals the potential value. This proves that both of them are optimal. In fact, the Hungarian method finds a perfect matching of '''tight edges''': an edge <math>ij</math> is called tight for a potential {{mvar|y}} if <math>y(i)+y(j) = c(i, j)</math>. Let us denote the [[Glossary of graph theory#Subgraphs|subgraph]] of '''tight''' edges by <math>G_y</math>. The cost of a perfect matching in <math>G_y</math> (if there is one) equals the value of {{mvar|y}}.
During the algorithm we maintain a potential {{mvar|y}} and an [[Glossary of graph theory#orientation|orientation]] of <math>G_y</math> (denoted by <math>\overrightarrow{G_y}</math>) which has the property that the edges oriented from {{mvar|T}} to {{mvar|S}} form a matching {{mvar|M}}. Initially, {{mvar|y}} is 0 everywhere, and all edges are oriented from {{mvar|S}} to {{mvar|T}} (so {{mvar|M}} is empty). In each step, either we modify {{mvar|y}} so that its value increases, or modify the orientation to obtain a matching with more edges. We maintain the invariant that all the edges of {{mvar|M}} are tight. We are done if {{mvar|M}} is a perfect matching.
Line 81 ⟶ 85:
In a general step, let <math>R_S \subseteq S</math> and <math>R_T \subseteq T</math> be the vertices not covered by {{mvar|M}} (so <math>R_S</math> consists of the vertices in {{mvar|S}} with no incoming edge and <math>R_T</math> consists of the vertices in {{mvar|T}} with no outgoing edge). Let {{mvar|Z}} be the set of vertices reachable in <math>\overrightarrow{G_y}</math> from <math>R_S</math> by a directed path. This can be computed by [[breadth-first search]].
If <math>R_T \cap Z</math> is nonempty, then reverse the orientation of all edges along a directed path in <math>\overrightarrow{G_y}</math> from <math>R_S</math> to <math>R_T</math>. Thus the size of the corresponding matching increases by 1.
If <math>R_T \cap Z</math> is empty, then let
Line 87 ⟶ 91:
:<math>\Delta := \min \{c(i,j)-y(i)-y(j): i \in Z \cap S, j \in T \setminus Z\}.</math>
{{math|Δ}} is well defined because at least one such edge <math>ij</math> must exist whenever the matching is not yet of maximum possible size (see the following section); it is positive because there are no tight edges between <math>Z \cap S</math> and <math>T \setminus Z</math>. Increase {{mvar|y}} by {{math|Δ}} on the vertices of <math>Z \cap S</math> and decrease {{mvar|y}} by {{math|Δ}} on the vertices of <math>Z \cap T</math>. The resulting {{mvar|y}} is still a potential, and although the graph <math>G_y</math> changes, it still contains {{mvar|M}} (see the next subsections). We orient the new edges from {{mvar|S}} to {{mvar|T}}. By the definition of {{math|Δ}} the set {{mvar|Z}} of vertices reachable from <math>R_S</math> increases (note that the number of tight edges does not necessarily increase). If the vertex added to <math>Z \cap T</math> is unmatched (that is, it is also in <math>R_T</math>), then at the next iteration the graph will have an augmenting path.
We repeat these steps until {{mvar|M}} is a perfect matching, in which case it gives a minimum cost assignment. The running time of this version of the method is <math>O(n^4)</math>: {{mvar|M}} is augmented {{mvar|n}} times, and in a phase where {{mvar|M}} is unchanged, there are at most {{mvar|n}} potential changes (since {{mvar|Z}} increases every time). The time sufficient for a potential change is <math>O(n^2)</math>.
Line 144 ⟶ 148:
*/
import std;
template <typename T, typename U>
using Pair = std::pair<T, U>;
template <typename T>
using Vector = std::vector<T>;
template <typename T>
using NumericLimits = std::numeric_limits<T>;
/**
* @brief Checks if b < a
*
* Sets a = min(a, b)
* @param a The first parameter to check
* @param b The second parameter to check
* @tparam The type to perform the check on
* @return true if b < a
*/
template <
constexpr bool ckmin(T return b < a ? a = b, } /**
* @brief Performs the Hungarian algorithm.
*
* Given J jobs and W workers (J <= W), computes the minimum cost to assign each
* prefix of jobs to distinct workers.
Line 168 ⟶ 188:
* to assign the first (j+1) jobs to distinct workers
*/
template <typename T>
const int J = (int)size(C), W = (int)size(C[0]);▼
const int J = static_cast<int>(C.size());
const int W = static_cast<int>(C[0].size());
assert(J <= W);
// job[w] = job assigned to w-th worker, or -1 if no job assigned
// note: a W-th worker was added for convenience
Vector<T> yt(W + 1); // potentials // -yt[W] will equal the sum of all deltas
const T inf =
for (int
int
job[
// min reduced cost over edges from Z to worker w
while (job[
const int j = job[
T delta = inf;
int
for (int w = 0; w < W; ++w) {
if (!
if (ckmin(
if (ckmin(delta,
wNext = w;
}
}
// delta will always be
// except possibly during the first time this loop runs
// if any entries of C[
for (int w = 0; w <= W; ++w) {
if (
yt[w] -= delta;
} else {
minTo[w] -= delta;
}▼
}
}
// update assignments along alternating path
for (int w;
job[wCur] = job[w = prev[wCur]];
answers.push_back(-yt[W]);
}
Line 214 ⟶ 243:
/**
* @brief Performs a sanity check for the Hungarian algorithm.
*
* Sanity check: https://en.wikipedia.org/wiki/Hungarian_algorithm#Example
* First job (5):
Line 222 ⟶ 253:
* First + second + third jobs (15):
* clean bathroom: Alice -> 8
* sweep floors:
* wash windows: Bob -> 3
*/
void
assert((hungarian(costs) ==
}
/**
void cordon_bleu() {▼
*/
int N, M;▼
vector<pair<int, int>> bottles(N), couriers(M);▼
Vector<Pair<int, int>> B(N);
for (auto &c : couriers) cin >> c.first >> c.second;▼
cin >> rest.first >> rest.second;▼
for (const Pair<int, int>& c: couriers)
Pair<int, int> rest;
▲ std::cin >> rest.first >> rest.second;
Vector<Vector<int>> costs(N, Vector<int>(N + M - 1));
auto dist = [&](const Pair<int, int>& x, const Pair<int, int>& y) -> int {
return std::abs(x.first - y.first) + std::abs(x.second - y.second);
};
for (int b = 0; b < N; ++b) {
for (int c = 0; c < M; ++c)
costs[b][c] = dist(couriers[c], bottles[b]) + dist(bottles[b], rest);
▲ }
▲ for (int _ = 0; _ < N - 1; ++_) { // restaurant -> bottle -> restaurant
costs[b][_ + M] = 2 * dist(bottles[b], rest);
}
}
/**
* @brief Entry point into the program.
*
* @return The return code of the program.
*/
int main() {
}
</syntaxhighlight>
Line 272 ⟶ 312:
| website = Algorithms for Competitive Programming
| access-date = May 14, 2023
}}</ref>
| url = https://cp-algorithms.com/graph/Assignment-problem-min-flow.html
| title = Solving assignment problem using min-cost-flow
Line 279 ⟶ 319:
| website = Algorithms for Competitive Programming
| access-date = May 14, 2023
}}</ref> where the reweighting technique from [[
The implementation from the previous section is rewritten below in such a way as to emphasize this
connection; it can be checked that the potentials <math>h</math> for workers <math>0\dots W-1</math> are equal to the potentials <math>y</math> from the previous solution up to a constant offset. When the graph is sparse (there are only <math>M</math> allowed job, worker pairs), it is possible
Line 286 ⟶ 326:
<syntaxhighlight lang="cpp" line="1">
template <typename T>
▲ const int J = (int)size(C), W = (int)size(C[0]);
const int W = static_cast<int>(C[0].size());
assert(J <= W);
// job[w] = job assigned to w-th worker, or -1 if no job assigned
// note: a W-th worker was added for convenience
T
const T inf =
// assign
for (int
int
job[
dist[W] = 0;
while (job[
T
vis[
int
// consider extending shortest path by
for (int w = 0; w < W; ++w) {
if (!vis[w]) {
// sum of reduced edge weights
T edge = C[job[
if (
edge -= C[job[
assert(edge >= 0); // consequence of Johnson potentials
}
if (ckmin(dist[w], dist[
if (ckmin(minDist, dist[w]))
wNext = w;
}
}
}
for (int w = 0; w < W; ++w) { // update potentials
ckmin(dist[w], dist[
h[w] += dist[w];
}
ans_cur += h[
for (int w;
answers.push_back(ansCur);
}
return answers;
Line 337 ⟶ 382:
==Matrix interpretation==
This variant of the algorithm follows the formulation given by Flood,<ref name="flood">{{cite journal | last=Flood | first=Merrill M. | title=The Traveling-Salesman Problem | journal=Operations Research | volume=4 | issue=1 | date=1956 | issn=0030-364X | doi=10.1287/opre.4.1.61 | pages=61–75}}</ref> and later described more explicitly by Munkres, who proved it runs in <math>\mathcal{O}(n^4)</math> time.<ref name="munkres"/> Instead of keeping track of the potentials of the vertices, the algorithm operates only on a matrix:
Given {{mvar|n}} workers and tasks, the problem is written in the form of an {{mvar|n}}×{{mvar|n}} matrix▼
:<math>a_{ij}:=c(i,j)-y(i)-y(j)</math>
where <math>c(i,j)</math> is the original cost matrix and <math>y(i), y(j)</math> are the potentials from the graph interpretation. Changing the potentials corresponds to adding or subtracting from rows or columns of this matrix. The algorithm starts with <math>a_{ij} = c(i, j)</math>. As such, it can be viewed as taking the original cost matrix and modifying it.
▲Given {{mvar|n}} workers and tasks, the problem is written in the form of an {{mvar|n}}×{{mvar|n}} cost matrix
:{{aligned table|cols=4|class=wikitable
Line 351 ⟶ 402:
===Step 1===
'''For each row, its minimum element is subtracted from every element in that row.''' This causes all elements to have
This also leads to at least one zero in each row. As such, a naive greedy algorithm can attempt to assign all workers a task with a penalty of zero. This is illustrated below.
Line 416 ⟶ 467:
|}
Find a non-covered zero and prime it
* If the zero is on the same row as a starred zero, cover the corresponding row, and uncover the column of the starred zero.
Line 513 ⟶ 564:
All zeros are now covered with a minimal number of rows and columns.
The aforementioned detailed description is ''just one way'' to draw the minimum number of lines to cover all the 0s. Other methods work as well.
===Step 5===
If the number of starred zeros is {{mvar|n}} (or in the general case <math>min(n,m)</math>, where {{mvar|n}} is the number of people and {{mvar|m}} is the number of jobs), the algorithm terminates. See the [[#Result|Result]] subsection below on how to interpret the results.
This is equivalent to subtracting a number from all rows which are not covered and adding the same number to all columns which are covered. These operations do not change optimal assignments.
===Result===
If following this specific version of the algorithm, the starred zeros form the minimum assignment.
From Kőnig's theorem,<ref>[[K%C5%91nig%27s theorem (graph theory)]] Konig's theorem</ref> the minimum number of lines (minimum vertex cover<ref>[[Vertex cover]] minimum vertex cover</ref>) will be {{mvar|n}} (the size of maximum matching<ref>[[Matching (graph theory)]] matching</ref>). Thus, when {{mvar|n}} lines are required, minimum cost assignment can be found by looking at only zeroes in the matrix.
==Bibliography==
* R.E. Burkard, M. Dell'Amico, S. Martello: ''Assignment Problems'' (Revised reprint). SIAM, Philadelphia (PA.) 2012. {{ISBN|978-1-61197-222-1}}
Line 550 ⟶ 602:
=== Implementations ===
Note that not all of these satisfy the <math>O(n^3)</math> time complexity, even if they claim so. Some may contain errors, implement the slower <math>O(n^4)</math> algorithm, or have other inefficiencies. In the worst case, a code example linked from Wikipedia could later be modified to include exploit code. Verification and benchmarking is necessary when using such code examples from unknown authors.
* [https://sourceforge.net/projects/tas/ Lua and Python versions of R.A. Pilgrim's code (claiming <math>O(n^3)</math> time complexity)]
* [https://github.com/Jonah-Heyl/Hungarian-algorithm/blob/main/HungarianAlgorithminJulia.ipynb Julia implementation]
* [https://github.com/maandree/hungarian-algorithm-n3/blob/master/hungarian.c C implementation claiming <math>O(n^3)</math> time complexity]
* [https://github.com/KevinStern/software-and-algorithms/blob/master/src/main/java/blogspot/software_and_algorithms/stern_library/optimization/HungarianAlgorithm.java Java implementation claiming <math>O(n^3)</math> time complexity]
Line 556 ⟶ 609:
* [https://github.com/evansenter/gene/blob/f515fd73cb9d6a22b4d4b146d70b6c2ec6a5125b/objects/extensions/hungarian.rb Ruby implementation with unit tests]
* [https://github.com/antifriz/hungarian-algorithm-n3 C# implementation claiming <math>O(n^3)</math> time complexity]
* [http://www.fantascienza.net/leonardo/so/hungarian.d D implementation with unit tests (port of a Java version claiming <math>O(n^3)</math>)] {{Webarchive|url=https://web.archive.org/web/20191230174110/http://www.fantascienza.net/leonardo/so/hungarian.d |date=30 December 2019 }}
* [http://www.ifors.ms.unimelb.edu.au/tutorial/hungarian/welcome_frame.html Online interactive implementation]<!--Please note that this implements a variant of the algorithm as described above. -->
* [http://www.netlib.org/utk/lsi/pcwLSI/text/node220.html Serial and parallel implementations.]
* [http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6543 Matlab and C] {{Webarchive|url=https://web.archive.org/web/20080503063209/http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6543 |date=3 May 2008 }}
* [https://metacpan.org/module/Algorithm::Munkres Perl implementation]
* [https://github.com/saebyn/munkres-cpp C++ implementation]
| |||