Recent from talks
Contribute something
Nothing was collected or created yet.
Hungarian algorithm
View on WikipediaThe Hungarian method is a combinatorial optimization algorithm that solves the assignment problem in polynomial time and which anticipated later 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 Hungarian mathematicians, Dénes Kőnig and Jenő Egerváry.[1][2] However, in 2006 it was discovered that Carl Gustav Jacobi had solved the assignment problem in the 19th century, and the solution had been published posthumously in 1890 in Latin.[3]
James Munkres reviewed the algorithm in 1957 and observed that it is (strongly) polynomial.[4] Since then the algorithm has been known also as the Kuhn–Munkres algorithm or Munkres assignment algorithm. The time complexity of the original algorithm was , however Edmonds and Karp, and independently Tomizawa, noticed that it can be modified to achieve an running time.[5][6] Ford and Fulkerson extended the method to general maximum flow problems in form of the Ford–Fulkerson algorithm.
The problem
[edit]Example
[edit]In this simple example, there are three workers: Alice, Bob and Carol. One of them has to clean the bathroom, another sweep the floors and the third washes the windows, but they each demand different pay for the various tasks. The problem is to find the lowest-cost way to assign the jobs. The problem can be represented in a matrix of the costs of the workers doing the jobs. For example:
- TaskWorker
Clean
bathroomSweep
floorsWash
windowsAlice $8 $4 $7 Bob $5 $2 $3 Carol $9 $4 $8
The Hungarian method, when applied to the above table, would give the minimum cost: this is $15, achieved by having Alice clean the bathroom, Carol sweep the floors, and Bob wash the windows. This can be confirmed using brute force:
(the unassigned person washes the windows)CleanSweepAlice Bob Carol Alice — $17 $16 Bob $18 — $18 Carol $15 $16 —
Matrix formulation
[edit]In the matrix formulation, we are given an n×n matrix, where the element in the i-th row and j-th column represents the cost of assigning the j-th job to the i-th worker. We have to find an assignment of the jobs to the workers, such that each job is assigned to one worker and each worker is assigned one job, such that the total cost of assignment is minimum.
This can be expressed as permuting the rows of a cost matrix C to minimize the trace of a matrix,
where P is a permutation matrix. (Equivalently, the columns can be permuted using CP.)
If the goal is to find the assignment that yields the maximum cost, the problem can be solved by negating the cost matrix C.
Bipartite graph formulation
[edit]The algorithm can equivalently be described by formulating the problem using a bipartite graph. We have a complete bipartite graph with n worker vertices (S) and n job vertices (T), and the edges (E) each have a cost . We want to find a perfect matching with a minimum total cost.
The algorithm in terms of bipartite graphs
[edit]Let us call a function a potential if for each .
The value of potential 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 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 is called tight for a potential y if . Let us denote the subgraph of tight edges by . The cost of a perfect matching in (if there is one) equals the value of y.
During the algorithm we maintain a potential y and an orientation of (denoted by ) which has the property that the edges oriented from T to S form a matching M. Initially, y is 0 everywhere, and all edges are oriented from S to T (so M is empty). In each step, either we modify 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 M are tight. We are done if M is a perfect matching.
In a general step, let and be the vertices not covered by M (so consists of the vertices in S with no incoming edge and consists of the vertices in T with no outgoing edge). Let Z be the set of vertices reachable in from by a directed path. This can be computed by breadth-first search.
If is nonempty, then reverse the orientation of all edges along a directed path in from to . Thus the size of the corresponding matching increases by 1.
If is empty, then let
Δ is well defined because at least one such edge 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 and . Increase y by Δ on the vertices of and decrease y by Δ on the vertices of . The resulting y is still a potential, and although the graph changes, it still contains M (see the next subsections). We orient the new edges from S to T. By the definition of Δ the set Z of vertices reachable from increases (note that the number of tight edges does not necessarily increase). If the vertex added to is unmatched (that is, it is also in ), then at the next iteration the graph will have an augmenting path.
We repeat these steps until M is a perfect matching, in which case it gives a minimum cost assignment. The running time of this version of the method is : M is augmented n times, and in a phase where M is unchanged, there are at most n potential changes (since Z increases every time). The time sufficient for a potential change is .
Proof that the algorithm makes progress
[edit]We must show that as long as the matching is not of maximum possible size, the algorithm is always able to make progress — that is, to either increase the number of matched edges, or tighten at least one edge. It suffices to show that at least one of the following holds at every step:
- M is of maximum possible size.
- contains an augmenting path.
- G contains a loose-tailed path: a path from some vertex in to a vertex in that consists of any number (possibly zero) of tight edges followed by a single loose edge. The trailing loose edge of a loose-tailed path is thus from , guaranteeing that Δ is well defined.
If M is of maximum possible size, we are of course finished. Otherwise, by Berge's lemma, there must exist an augmenting path P with respect to M in the underlying graph G. However, this path may not exist in : Although every even-numbered edge in P is tight by the definition of M, odd-numbered edges may be loose and thus absent from . One endpoint of P is in , the other in ; w.l.o.g., suppose it begins in . If every edge on P is tight, then it remains an augmenting path in and we are done. Otherwise, let be the first loose edge on P. If then we have found a loose-tailed path and we are done. Otherwise, v is reachable from some other path Q of tight edges from a vertex in . Let be the subpath of P beginning at v and continuing to the end, and let be the path formed by traveling along Q until a vertex on is reached, and then continuing to the end of . Observe that is an augmenting path in G with at least one fewer loose edge than P. P can be replaced with and this reasoning process iterated (formally, using induction on the number of loose edges) until either an augmenting path in or a loose-tailed path in G is found.
Proof that adjusting the potential y leaves M unchanged
[edit]To show that every edge in M remains after adjusting y, it suffices to show that for an arbitrary edge in M, either both of its endpoints, or neither of them, are in Z. To this end let be an edge in M from T to S. It is easy to see that if v is in Z then u must be too, since every edge in M is tight. Now suppose, toward contradiction, that but . u itself cannot be in because it is the endpoint of a matched edge, so there must be some directed path of tight edges from a vertex in to u. This path must avoid v, since that is by assumption not in Z, so the vertex immediately preceding u in this path is some other vertex . is a tight edge from T to S and is thus in M. But then M contains two edges that share the vertex u, contradicting the fact that M is a matching. Thus every edge in M has either both endpoints or neither endpoint in Z.
Proof that y remains a potential
[edit]To show that y remains a potential after being adjusted, it suffices to show that no edge has its total potential increased beyond its cost. This is already established for edges in M by the preceding paragraph, so consider an arbitrary edge uv from S to T. If is increased by Δ, then either , in which case is decreased by Δ, leaving the total potential of the edge unchanged, or , in which case the definition of Δ guarantees that . Thus y remains a potential.
The algorithm in O(n3) time
[edit]Suppose there are jobs and workers (). We describe how to compute for each prefix of jobs the minimum total cost to assign each of these jobs to distinct workers. Specifically, we add the th job and update the total cost in time , yielding an overall time complexity of . Note that this is better than when the number of jobs is small relative to the number of workers.
Adding the j-th job in O(jW) time
[edit]We use the same notation as the previous section, though we modify their definitions as necessary. Let denote the set of the first jobs and denote the set of all workers.
Before the th step of the algorithm, we assume that we have a matching on that matches all jobs in and potentials satisfying the following condition: the matching is tight with respect to the potentials, and the potentials of all unmatched workers are zero, and the potentials of all matched workers are non-positive. Note that such potentials certify the optimality of the matching.
During the th step, we add the th job to to form and initialize . At all times, every vertex in will be reachable from the th job in . While does not contain a worker that has not been assigned a job, let
and denote any at which the minimum is attained. After adjusting the potentials in the way described in the previous section, there is now a tight edge from to .
- If is unmatched, then we have an augmenting path in the subgraph of tight edges from to . After toggling the matching along this path, we have now matched the first jobs, and this procedure terminates.
- Otherwise, we add and the job matched with it to .
Adjusting potentials takes time. Recomputing and after changing the potentials and also can be done in time. Case 1 can occur at most times before case 2 occurs and the procedure terminates, yielding the overall time complexity of .
Implementation in C++
[edit]For convenience of implementation, the code below adds an additional worker such that stores the negation of the sum of all computed so far. After the th job is added and the matching updated, the cost of the current matching equals the sum of all computed so far, or .
This code is adapted from e-maxx :: algo.[7]
/**
* Solution to https://open.kattis.com/problems/cordonbleu using Hungarian
* algorithm.
*/
import <cassert>;
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 <typename T>
constexpr bool ckmin(T& a, const T& b) {
return b < a ? a = b, true : false;
}
/**
* @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.
*
* @tparam T a type large enough to represent integers on the order of J *
* max(|C|)
* @param C a matrix of dimensions JxW such that C[j][w] = cost to assign j-th
* job to w-th worker (possibly negative)
*
* @return a vector of length J, with the j-th entry equaling the minimum cost
* to assign the first (j+1) jobs to distinct workers
*/
template <typename T>
Vector<T> hungarian(const Vector<Vector<T>>& C) {
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<int> job(W + 1, -1);
Vector<T> ys(J);
Vector<T> yt(W + 1); // potentials
// -yt[W] will equal the sum of all deltas
Vector<T> answers;
const T inf = NumericLimits<T>::max();
for (int jCur = 0; jCur < J; ++jCur) { // assign jCur-th job
int wCur = W;
job[wCur] = jCur;
// min reduced cost over edges from Z to worker w
Vector<T> minTo(W + 1, inf);
Vector<int> prev(W + 1, -1); // previous worker on alternating path
Vector<bool> inZ(W + 1); // whether worker is in Z
while (job[wCur] != -1) { // runs at most jCur + 1 times
inZ[wCur] = true;
const int j = job[wCur];
T delta = inf;
int wNext;
for (int w = 0; w < W; ++w) {
if (!inZ[w]) {
if (ckmin(minTo[w], C[j][w] - ys[j] - yt[w]))
prev[w] = wCur;
if (ckmin(delta, minTo[w]))
wNext = w;
}
}
// delta will always be nonnegative,
// except possibly during the first time this loop runs
// if any entries of C[jCur] are negative
for (int w = 0; w <= W; ++w) {
if (inZ[w]) {
ys[job[w]] += delta;
yt[w] -= delta;
} else {
minTo[w] -= delta;
}
}
wCur = wNext;
}
// update assignments along alternating path
for (int w; wCur != W; wCur = w)
job[wCur] = job[w = prev[wCur]];
answers.push_back(-yt[W]);
}
return answers;
}
/**
* @brief Performs a sanity check for the Hungarian algorithm.
*
* Sanity check: https://en.wikipedia.org/wiki/Hungarian_algorithm#Example
* First job (5):
* clean bathroom: Bob -> 5
* First + second jobs (9):
* clean bathroom: Bob -> 5
* sweep floors: Alice -> 4
* First + second + third jobs (15):
* clean bathroom: Alice -> 8
* sweep floors: Carol -> 4
* wash windows: Bob -> 3
*/
void sanityCheckHungarian() {
Vector<Vector<int>> costs{{8, 5, 9}, {4, 2, 4}, {7, 3, 8}};
assert((hungarian(costs) == Vector<int>{5, 9, 15}));
std::println(stderr, "Sanity check passed.");
}
/**
* @brief solves https://open.kattis.com/problems/cordonbleu
*/
void cordonBleu() {
int N;
int M;
std::cin >> N >> M;
Vector<Pair<int, int>> B(N);
Vector<Pair<int, int>> C(M);
Vector<Pair<int, int>> bottles(N);
Vector<Pair<int, int>> couriers(M);
for (const Pair<int, int>& b: bottles)
std::cin >> b.first >> b.second;
for (const Pair<int, int>& c: couriers)
std::cin >> c.first >> c.second;
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) // courier -> bottle -> restaurant
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);
}
std::println(hungarian(costs).back());
}
/**
* @brief Entry point into the program.
*
* @return The return code of the program.
*/
int main() {
sanityCheckHungarian();
cordonBleu();
}
Connection to successive shortest paths
[edit]The Hungarian algorithm can be seen to be equivalent to the successive shortest path algorithm for minimum cost flow,[8][9] where the reweighting technique from Johnson's algorithm is used to find the shortest paths. 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 for workers are equal to the potentials from the previous solution up to a constant offset. When the graph is sparse (there are only allowed job, worker pairs), it is possible to optimize this algorithm to run in time by using a Fibonacci heap to determine instead of iterating over all workers to find the one with minimum distance (alluded to here).
template <typename T>
Vector<T> hungarian(const Vector<Vector<T>>& C) {
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<int> job(W + 1, -1);
Vector<T> h(W); // Johnson potentials
Vector<T> answers;
T ansCur = 0;
const T inf = NumericLimits<T>::max();
// assign jCur-th job using Dijkstra with potentials
for (int jCur = 0; jCur < J; ++jCur) {
int wCur = W; // unvisited worker with minimum distance
job[wCur] = jCur;
Vector<T> dist(W + 1, inf); // Johnson-reduced distances
dist[W] = 0;
Vector<bool> vis(W + 1); // whether visited yet
Vector<int> prev(W + 1, -1); // previous worker on shortest path
while (job[wCur] != -1) { // Dijkstra step: pop min worker from heap
T minDist = inf;
vis[wCur] = true;
int wNext = -1; // next unvisited worker with minimum distance
// consider extending shortest path by wCur -> job[wCur] -> w
for (int w = 0; w < W; ++w) {
if (!vis[w]) {
// sum of reduced edge weights wCur -> job[wCur] -> w
T edge = C[job[wCur]][w] - h[w];
if (wCur != W) {
edge -= C[job[wCur]][wCur] - h[wCur];
assert(edge >= 0); // consequence of Johnson potentials
}
if (ckmin(dist[w], dist[wCur] + edge))
prev[w] = wCur;
if (ckmin(minDist, dist[w]))
wNext = w;
}
}
wCur = wNext;
}
for (int w = 0; w < W; ++w) { // update potentials
ckmin(dist[w], dist[wCur]);
h[w] += dist[w];
}
ans_cur += h[wCur];
for (int w; wCur != W; wCur = w)
job[wCur] = job[w = prev[wCur]];
answers.push_back(ansCur);
}
return answers;
}
Matrix interpretation
[edit]This variant of the algorithm follows the formulation given by Flood,[10] and later described more explicitly by Munkres, who proved it runs in time.[4] Instead of keeping track of the potentials of the vertices, the algorithm operates only on a matrix:
where is the original cost matrix and 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 . As such, it can be viewed as taking the original cost matrix and modifying it.
Given n workers and tasks, the problem is written in the form of an n×n cost matrix
a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 d3 d4
where a, b, c and d are workers who have to perform tasks 1, 2, 3 and 4. a1, a2, a3, and a4 denote the penalties incurred when worker "a" does task 1, 2, 3, and 4 respectively.
The problem is equivalent to assigning each worker a unique task such that the total penalty is minimized. Note that each task can only be worked on by one worker.
Step 1
[edit]For each row, its minimum element is subtracted from every element in that row. This causes all elements to have nonnegative values. Therefore, an assignment with a total penalty of 0 is by definition a minimum assignment.
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.
0 a2 a3 a4 b1 b2 b3 0 c1 0 c3 c4 d1 d2 0 d4
The zeros above would be the assigned tasks.
Worst-case there are n! combinations to try, since multiple zeroes can appear in a row if multiple elements are the minimum. So at some point this naive algorithm should be short circuited.
Step 2
[edit]Sometimes it may turn out that the matrix at this stage cannot be used for assigning, as is the case for the matrix below.
0 a2 0 a4 b1 0 b3 0 0 c2 c3 c4 0 d2 d3 d4
To overcome this, we repeat the above procedure for all columns (i.e. the minimum element in each column is subtracted from all the elements in that column) and then check if an assignment with penalty 0 is possible.
In most situations this will give the result, but if it is still not possible then we need to keep going.
Step 3
[edit]All zeros in the matrix must be covered by marking as few rows and/or columns as possible. Steps 3 and 4 form one way to accomplish this.
For each row, try to assign an arbitrary zero. Assigned tasks are represented by starring a zero. Note that assignments can't be in the same row or column.
- We assign the first zero of Row 1. The second zero of Row 1 can't be assigned.
- We assign the first zero of Row 2. The second zero of Row 2 can't be assigned.
- Zeros on Row 3 and Row 4 can't be assigned, because they are on the same column as the zero assigned on Row 1.
We could end with another assignment if we choose another ordering of the rows and columns.
0* a2 0 a4 b1 0* b3 0 0 c2 c3 c4 0 d2 d3 d4
Step 4
[edit]Cover all columns containing a (starred) zero.
× × 0* a2 0 a4 b1 0* b3 0 0 c2 c3 c4 0 d2 d3 d4
Find a non-covered zero and prime it (mark it with a prime symbol). If no such zero can be found, meaning all zeroes are covered, skip to step 5.
- If the zero is on the same row as a starred zero, cover the corresponding row, and uncover the column of the starred zero.
- Then, GOTO "Find a non-covered zero and prime it."
- Here, the second zero of Row 1 is uncovered. Because there is another zero starred on Row 1, we cover Row 1 and uncover Column 1.
- Then, the second zero of Row 2 is uncovered. We cover Row 2 and uncover Column 2.
× 0* a2 0' a4 × b1 0* b3 0 0 c2 c3 c4 0 d2 d3 d4
0* a2 0' a4 × b1 0* b3 0' × 0 c2 c3 c4 0 d2 d3 d4
- Else the non-covered zero has no assigned zero on its row. We make a path starting from the zero by performing the following steps:
- Substep 1: Find a starred zero on the corresponding column. If there is one, go to Substep 2, else, stop.
- Substep 2: Find a primed zero on the corresponding row (there should always be one). Go to Substep 1.
The zero on Row 3 is uncovered. We add to the path the first zero of Row 1, then the second zero of Row 1, then we are done.
0* a2 0' a4 × b1 0* b3 0' × 0' c2 c3 c4 0 d2 d3 d4
- (Else branch continued) For all zeros encountered during the path, star primed zeros and unstar starred zeros.
- As the path begins and ends by a primed zero when swapping starred zeros, we have assigned one more zero.
0 a2 0* a4 b1 0* b3 0 0* c2 c3 c4 0 d2 d3 d4
- (Else branch continued) Unprime all primed zeroes and uncover all lines.
- Repeat the previous steps (continue looping until the above "skip to step 5" is reached).
- We cover columns 1, 2 and 3. The second zero on Row 2 is uncovered, so we cover Row 2 and uncover Column 2:
× × 0 a2 0* a4 b1 0* b3 0' × 0* c2 c3 c4 0 d2 d3 d4
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
[edit]If the number of starred zeros is n (or in the general case , where n is the number of people and m is the number of jobs), the algorithm terminates. See the Result subsection below on how to interpret the results.
Otherwise, find the lowest uncovered value. Subtract this from every uncovered element and add it to every element covered by two lines. Go back to step 4.
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
[edit]If following this specific version of the algorithm, the starred zeros form the minimum assignment.
From Kőnig's theorem,[11] the minimum number of lines (minimum vertex cover[12]) will be n (the size of maximum matching[13]). Thus, when n lines are required, minimum cost assignment can be found by looking at only zeroes in the matrix.
Bibliography
[edit]- R.E. Burkard, M. Dell'Amico, S. Martello: Assignment Problems (Revised reprint). SIAM, Philadelphia (PA.) 2012. ISBN 978-1-61197-222-1
- M. Fischetti, "Lezioni di Ricerca Operativa", Edizioni Libreria Progetto Padova, Italia, 1995.
- R. Ahuja, T. Magnanti, J. Orlin, "Network Flows", Prentice Hall, 1993.
- S. Martello, "Jeno Egerváry: from the origins of the Hungarian algorithm to satellite communication". Central European Journal of Operational Research 18, 47–58, 2010
References
[edit]- ^ Harold W. Kuhn, "The Hungarian Method for the assignment problem", Naval Research Logistics Quarterly, 2: 83–97, 1955. Kuhn's original publication.
- ^ Harold W. Kuhn, "Variants of the Hungarian method for assignment problems", Naval Research Logistics Quarterly, 3: 253–258, 1956.
- ^ "Presentation". Archived from the original on 16 October 2015.
- ^ a b J. Munkres, "Algorithms for the Assignment and Transportation Problems", Journal of the Society for Industrial and Applied Mathematics, 5(1):32–38, 1957 March.
- ^ Edmonds, Jack; Karp, Richard M. (1 April 1972). "Theoretical Improvements in Algorithmic Efficiency for Network Flow Problems". Journal of the ACM. 19 (2): 248–264. doi:10.1145/321694.321699. S2CID 6375478.
- ^ Tomizawa, N. (1971). "On some techniques useful for solution of transportation network problems". Networks. 1 (2): 173–194. doi:10.1002/net.3230010206. ISSN 1097-0037.
- ^ "Hungarian Algorithm for Solving the Assignment Problem". e-maxx :: algo. 23 August 2012. Retrieved 13 May 2023.
- ^ Jacob Kogler (20 December 2022). "Minimum-cost flow - Successive shortest path algorithm". Algorithms for Competitive Programming. Retrieved 14 May 2023.
- ^ "Solving assignment problem using min-cost-flow". Algorithms for Competitive Programming. 17 July 2022. Retrieved 14 May 2023.
- ^ Flood, Merrill M. (1956). "The Traveling-Salesman Problem". Operations Research. 4 (1): 61–75. doi:10.1287/opre.4.1.61. ISSN 0030-364X.
- ^ Kőnig's theorem (graph theory) Konig's theorem
- ^ Vertex cover minimum vertex cover
- ^ Matching (graph theory) matching
External links
[edit]- Bruff, Derek, The Assignment Problem and the Hungarian Method (matrix formalism).
- Mordecai J. Golin, Bipartite Matching and the Hungarian Method (bigraph formalism), Course Notes, Hong Kong University of Science and Technology.
- Hungarian maximum matching algorithm (both formalisms), in Brilliant website.
- 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.
Implementations
[edit]Note that not all of these satisfy the time complexity, even if they claim so. Some may contain errors, implement the slower 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.
- Lua and Python versions of R.A. Pilgrim's code (claiming time complexity)
- Julia implementation
- C implementation claiming time complexity
- Java implementation claiming time complexity
- Python implementation
- Ruby implementation with unit tests
- C# implementation claiming time complexity
- D implementation with unit tests (port of a Java version claiming ) Archived 30 December 2019 at the Wayback Machine
- Online interactive implementation
- Serial and parallel implementations.
- Matlab and C Archived 3 May 2008 at the Wayback Machine
- Perl implementation
- C++ implementation
- C++ implementation claiming time complexity (BSD style open source licensed)
- MATLAB implementation
- C implementation
- JavaScript implementation with unit tests (port of a Java version claiming time complexity)
- Clue R package proposes an implementation, solve_LSAP
- Node.js implementation on GitHub
- Python implementation in scipy package
Hungarian algorithm
View on GrokipediaThe Assignment Problem
Illustrative Example
Consider a simple assignment problem involving three workers (W1, W2, W3) and three tasks (T1, T2, T3), where the goal is to assign each worker to a unique task to minimize the total cost. The costs for assigning each worker to each task are given in the following matrix:| Worker | T1 | T2 | T3 |
|---|---|---|---|
| W1 | 25 | 44 | 36 |
| W2 | 28 | 41 | 40 |
| W3 | 23 | 50 | 35 |
- Row W1 minimum: 25, resulting row: [0, 19, 11]
- Row W2 minimum: 28, resulting row: [0, 13, 12]
- Row W3 minimum: 23, resulting row: [0, 27, 12]
| Worker | T1 | T2 | T3 |
|---|---|---|---|
| W1 | 0 | 19 | 11 |
| W2 | 0 | 13 | 12 |
| W3 | 0 | 27 | 12 |
- Column T1 minimum: 0, no change
- Column T2 minimum: 13, resulting column: [6, 0, 14]
- Column T3 minimum: 11, resulting column: [0, 1, 1]
| Worker | T1 | T2 | T3 |
|---|---|---|---|
| W1 | 0 | 6 | 0 |
| W2 | 0 | 0 | 1 |
| W3 | 0 | 14 | 1 |
| Worker | T1 | T2 | T3 |
|---|---|---|---|
| W1 | 25 | 44 | 36 |
| W2 | 28 | 41 | 40 |
| W3 | 23 | 50 | 35 |
Cost Matrix Formulation
The assignment problem can be formally defined using an cost matrix , where represents the cost of assigning the -th row entity (such as an agent or worker) to the -th column entity (such as a task or job), for .[1] The objective is to find a one-to-one assignment that minimizes the total cost, formulated as the integer linear program: subject to the constraints where if row is assigned to column , and otherwise.[9][1] These constraints ensure that each row and each column is assigned exactly once, corresponding to a perfect matching in the underlying structure.[1] This formulation assumes a balanced assignment problem, where the number of rows equals the number of columns ( agents and tasks), resulting in a square cost matrix.[1] In the unbalanced case, where the numbers differ (e.g., more tasks than agents), the problem can be converted to a balanced one by introducing dummy rows or columns with zero costs to balance the matrix and account for unassigned entities at no additional cost, though the standard Hungarian algorithm targets the balanced scenario.[10] A brute-force approach to solving this problem requires evaluating all possible permutations of assignments to identify the minimum-cost one, yielding exponential time complexity that becomes impractical for even moderate .[11] This motivates the development of polynomial-time algorithms like the Hungarian method to efficiently compute the optimal solution.[1]Bipartite Matching Formulation
The assignment problem is equivalently formulated as the problem of finding a minimum-weight perfect matching in a complete bipartite graph , where and are disjoint vertex sets each of size , and consists of all possible edges between and .[1] Each edge is assigned a weight , corresponding to the cost of assigning the element from to the element from ; for maximization objectives, weights can be negated to convert to a minimum-cost problem.[1] The objective is to select a matching such that every vertex in is incident to exactly one edge in (a perfect matching, required since the graph is balanced with ), minimizing the total weight . This graph-theoretic perspective bridges the algebraic cost matrix representation to broader matching algorithms in combinatorial optimization. The algorithm solving this problem is commonly known as the Hungarian algorithm, or alternatively the Kuhn-Munkres algorithm, named after Harold W. Kuhn, who formalized it in 1955, and James Munkres, who provided a simplified version in 1957; the "Hungarian" moniker honors the foundational contributions of Hungarian mathematicians Dénes Kőnig and Jenő Egerváry to related matching theory.[1][12]Core Algorithm Concepts
Feasibility and Dual Variables
In the Hungarian algorithm, which solves the minimum cost perfect matching problem in a complete bipartite graph with vertex partitions and and nonnegative edge costs , feasible potentials are defined as real-valued labels for each and for each such that for all , . These potentials induce reduced costs for every edge, transforming the original cost matrix into a nonnegative equivalent while preserving optimality of solutions.[13] The equality subgraph associated with a set of feasible potentials includes only those edges where the reduced cost vanishes, i.e., or . A perfect matching within this subgraph identifies candidate assignments that are tight with respect to the potentials, forming the basis for primal feasibility in the algorithm's iterative process.[13] Feasible potentials correspond directly to the dual variables of the linear programming formulation of the assignment problem. The primal problem is subject to Its dual is subject to with unrestricted in sign; feasible potentials thus provide a feasible dual solution whose objective value lower-bounds the primal minimum cost. By the strong duality theorem for linear programs, optimal potentials achieve equality with the optimal assignment cost.[13][1] A simple initial feasible solution sets for each and for each , ensuring the dual constraints hold since the minimum over is at most any specific . This initialization creates a starting equality subgraph with at least one zero per row, facilitating the algorithm's early iterations.[13]Equality Subgraph and Augmenting Paths
In the Hungarian algorithm for solving the minimum cost assignment problem on a bipartite graph with partitions and , the equality subgraph is constructed using a set of feasible dual variables or potentials for and for , where for all edges . This subgraph includes all vertices and only those edges where the reduced cost , ensuring the edges represent zero-cost opportunities relative to the current potentials.[1][14] Within , alternating paths are defined with respect to a partial matching . An alternating path begins at an unmatched vertex in and proceeds as a sequence of edges that alternate between non-matching edges (not in ) and matching edges (in ), all contained in . These paths explore the structure of feasible zero-cost edges while respecting the current assignment constraints.[1] An augmenting path is a specific type of alternating path in that terminates at an unmatched vertex in . Augmenting along such a path involves taking the symmetric difference , which flips the matched and unmatched edges along , thereby increasing the size of the matching by one without increasing the total cost, as all edges in have zero reduced cost.[1][14] The Hungarian algorithm builds a maximum cardinality matching—and ultimately an optimal assignment—by repeatedly identifying and augmenting along these paths in the equality subgraph, adjusting the potentials when no further augmenting paths exist to enlarge the subgraph and create new zero-reduced-cost edges, until a perfect matching is obtained. This process leverages the zero-reduced-cost structure of to ensure feasibility and optimality in the primal-dual framework.[1]Step-by-Step Algorithm in Bipartite Graphs
Initialization and Potential Adjustment
The Hungarian algorithm for solving the minimum cost assignment problem initializes the dual variables, or potentials, to establish a feasible starting point with non-negative reduced costs. The potentials for vertices in the left partition U are set to for each , while the potentials for vertices in the right partition V are set to for all . This ensures that the reduced costs for every edge , satisfying the dual feasibility constraints . The initial matching is empty.[15] When no augmenting path is found in the equality subgraph during a search phase, the potentials are adjusted to enlarge the equality subgraph and restore the possibility of finding augmenting paths. The search builds a tree of reachable vertices from free U vertices. The adjustment computes as the minimum reduced cost over edges from reached (visited) vertices in U to unreached vertices in V. The potentials are then updated by adding to for all reached vertices and subtracting from for all reached vertices . This uniform shift tightens edges in the search tree, making some previously non-tight edges tight (zero reduced cost), while preserving non-negativity of reduced costs elsewhere. The equality subgraph consists of tight edges (zero reduced cost) under the current potentials.[15] This update mechanism increases the dual objective by (or adjusted for reached sets in some implementations), where is the number of reached vertices, but ensures progress toward optimality by incorporating the slack. The process repeats until a perfect matching is found.[15]Finding and Augmenting Along Shortest Paths
In the Hungarian algorithm, finding an augmenting path occurs within the equality subgraph, which comprises all edges where the reduced cost , with and denoting the dual potentials (also called labels) for vertices and , respectively.[15] This subgraph encodes feasible alternating paths relative to the current matching , and the search begins from all free (unmatched) vertices in .[15] To locate the shortest augmenting path—defined by the minimal number of edges, as all edges in the equality subgraph carry zero reduced cost—a breadth-first search (BFS) is employed, treating the subgraph as unweighted for path length purposes.[15] The BFS operates in layers: the initial layer includes all free vertices, followed by reachable unmatched edges to vertices, then backward along matched edges to , and so forth, alternating directions to maintain the path's feasibility.[15] This layered approach handles multiple free vertices simultaneously by initiating the search from all of them as multi-sources, ensuring the first reachable free vertex yields a shortest path without prioritizing one source over others.[15] Alternatively, for implementations emphasizing the weighted structure, a modified Dijkstra's algorithm can be applied directly on the residual graph with reduced costs as edge weights, where potentials serve as node labels to guarantee non-negativity ( for forward edges and the negated for backward).[16] In this setup, the priority queue orders vertices by cumulative reduced cost distance from the sources (free vertices), and the algorithm identifies the minimal-cost augmenting path. Potential adjustments from prior steps ensure no negative cycles or weights disrupt the search.[16] Upon discovering an augmenting path from a free vertex to a free vertex, augmentation proceeds by symmetric difference: unmatched edges in are added to , and matched edges are removed, resulting in a new matching with cardinality .[15] This flip operation updates the assignment without altering the total cost under current potentials, as the path's reduced cost sums to zero. The process repeats in the same equality subgraph until no more augmenting paths exist, at which point potentials are adjusted if necessary.[15]Progress and Termination Guarantees
The Hungarian algorithm progresses by iteratively identifying and augmenting along shortest augmenting paths in the equality subgraph defined by the current dual variables, with each such augmentation strictly increasing the cardinality of the matching by one.[15] In a balanced bipartite graph with vertices per partition, this process performs at most augmentations before attaining a perfect matching of size . Multiple augmentations may occur in a single phase with fixed potentials until no more paths are found in the current equality subgraph.[15] When no augmenting path exists in the equality subgraph for the current phase, the dual variables for left vertices and for right vertices are updated using the minimum slack from the search tree to maintain feasibility of the reduced costs for all edges. This adjustment raises the dual objective by an amount related to , which is positive and ensures new tight edges are created while preserving nonnegativity elsewhere.[15] The algorithm terminates when no augmenting path exists in the equality subgraph and all vertices are matched. At this point, is a maximum matching by Berge's lemma, which asserts that a matching in a graph is maximum if and only if no augmenting path exists with respect to it. Moreover, the dual solution remains feasible, and complementary slackness holds: reduced costs are zero for all edges in , and there are no zero-reduced-cost edges incident to unmatched vertices (though in perfect matching, none unmatched), equating the primal and dual objectives and confirming 's optimality.[15] Termination occurs after at most augmentations, with phases bounded by the use of shortest paths to limit phase length. Each phase requires time to detect and trace augmenting paths via breadth-first search in the equality subgraph.[15] This bounded iteration count, combined with progress per phase, guarantees finite completion with an optimal assignment.[15]Time Complexity and Implementation
O(n^3) Implementation Details
The O(n³) implementation of the Hungarian algorithm operates on a complete bipartite graph with n vertices on each side, using dual variables (potentials) to maintain feasibility and compute reduced costs for edges. The algorithm consists of up to n phases, each finding a shortest augmenting path in the residual graph with respect to reduced costs and augmenting the matching along it; each phase requires O(n²) time for path finding and updates, yielding the overall O(n³) time complexity.[17] Key data structures include arrays for row potentials and column potentials (both of size n), a matching array and (size n) to track assignments, and temporary arrays such as (size n) to reconstruct paths, (size n) to mark visited columns, and (size n) for slack values. The cost matrix (n × n) stores original edge weights, while reduced costs are computed on-the-fly as to ensure non-negativity, allowing efficient shortest path computations without explicit adjacency lists for dense graphs.[17] Slack variables in track the minimum reduced cost from the current search tree's frontier (unmatched rows) to each column j in the right partition V, updated during the path search to identify the bottleneck delta for potential adjustments. This enables a specialized O(n²) shortest path algorithm akin to Dijkstra, but implemented via iterative minimum selection over slacks rather than a full priority queue, avoiding logarithmic factors.[17] The main loop initializes feasible potentials (e.g., , ) and an empty matching, then iterates until the matching size reaches n. In each phase, starting from a free row, the algorithm builds an alternating tree in the equality subgraph (edges with zero reduced cost), iteratively adjusting potentials by the minimum slack to expand the tree until an augmenting path to a free column is found, followed by augmentation and potential updates to ensure progress.[17] Potential updates occur after each augmentation by the minimum slack along the path, ensuring the equality subgraph grows and progress is made.[17]Incremental Addition of Assignments
The incremental variant of the Hungarian algorithm addresses dynamic scenarios where the assignment problem evolves by adding one new job (or worker) at a time to an existing bipartite graph, allowing reuse of the prior optimal matching and dual variables to avoid recomputing the entire solution from scratch. This approach begins with a solved k × k assignment problem, which provides an optimal matching and feasible dual potentials (labels) for the vertices. When adding the (k+1)th row and column—representing the new job and its associated costs to existing workers—the algorithm initializes the dual potentials for the new vertices based on maximum weight differences relative to the existing structure, ensuring initial feasibility. The previous matching remains partially valid, as it saturates the first k vertices on both sides, leaving only the new vertices unsaturated. To update the matching, the algorithm focuses computation on the new elements by constructing the equality subgraph using the reused potentials and searching for a single augmenting path that incorporates the new job, typically via a Hungarian tree rooted at the unsaturated new vertex. If no immediate path exists, potentials are adjusted by computing the minimum slack (the smallest difference between dual sums and edge costs) over relevant unmatched edges, which expands the equality subgraph and enables path discovery without altering the entire dual solution. This adjustment preserves the optimality of the prior matching while integrating the new assignment, resulting in a cardinality increase of one. The process re-runs only aspects involving the new column, leveraging the structure of the previous solution to limit exploration. The time complexity of this incremental addition is O(n²) per update in the worst case, where n is the current graph size, achieved through efficient slack maintenance and tree growth bounded by the graph's density; over n sequential additions starting from a 1 × 1 problem, the total complexity amortizes to O(n³), matching the standard Hungarian algorithm but offering efficiency for partial updates. Earlier implementations may exhibit O(n W) behavior if using unscaled shortest path methods with integer costs up to W, but polynomial-time scaling techniques, such as those employing Dijkstra's algorithm with potentials, ensure consistent O(n³) overall without dependence on cost magnitude. An improved variant further reduces average-case time to O(n) for certain configurations by optimizing dual initialization and path detection, though it falls back to O(n²) generally.[18] This incremental approach finds applications in online assignment scenarios, such as dynamic resource allocation in radar tracking systems where targets appear sequentially, or interactive topic modeling where mappings update based on user feedback without full recomputation. By enabling efficient handling of evolving inputs, it supports real-time decision-making in environments like scheduling or matching markets with arriving entities.[19][20]Theoretical Connections
Relation to Successive Shortest Paths
The assignment problem can be modeled as a minimum-cost flow problem in a bipartite flow network with unit capacities. A source vertex is connected to each vertex in the left partition by an edge of capacity 1 and cost 0; each vertex is connected to each vertex by an edge of capacity 1 and cost ; and each vertex in the right partition is connected to a sink vertex by an edge of capacity 1 and cost 0. The goal is to compute a minimum-cost flow of value from to . The successive shortest paths algorithm addresses this min-cost flow instance by beginning with a zero flow and repeatedly identifying the shortest - path in the residual graph, where edge costs are adjusted using node potentials to obtain non-negative reduced costs . Flow is augmented by 1 along this path, potentials are updated based on the path distances to preserve reduced cost non-negativity, and the process repeats for iterations until the maximum flow is achieved. The Hungarian algorithm is equivalent to the successive shortest paths algorithm when specialized to this unit-capacity bipartite network for the assignment problem. The row and column potentials in the Hungarian method function as the node potentials for reduced cost computations; the Hungarian reduced costs match the flow network's reduced costs on forward and backward residual edges. Furthermore, the equality subgraph in the Hungarian algorithm, formed by edges with zero reduced cost, corresponds exactly to the subgraph of zero-cost residual edges used to find shortest augmenting paths via breadth-first search.[21] This equivalence highlights how the Hungarian algorithm optimizes the successive shortest paths framework by exploiting the bipartite structure and unit capacities, thereby eliminating the need for general-purpose shortest path computations like Dijkstra's algorithm and focusing instead on efficient augmenting path searches within the equality subgraph to achieve time complexity without broader min-cost flow overhead.Linear Programming Duality Interpretation
The Hungarian algorithm for solving the assignment problem can be interpreted through the framework of linear programming (LP) duality, where the primal problem corresponds to finding a minimum-cost perfect matching in a bipartite graph, and the dual provides bounds via node potentials. The primal LP formulation is as follows: where is the cost of assigning row to column , and indicates the assignment (with for integral solutions). This LP has an integral optimal solution because the constraint matrix (the incidence matrix of the complete bipartite graph ) is totally unimodular, ensuring that the polytope defined by the constraints has integer vertices.[22][23] The associated dual LP is: where and represent unrestricted potentials (or labels) for rows and columns, respectively. By weak duality, any feasible primal solution provides an upper bound on the dual objective, and vice versa.[1][23] The Hungarian algorithm operates as a primal-dual method, maintaining a feasible primal solution (a partial matching, which is integral) and a feasible dual solution (feasible potentials ) throughout its execution. It iteratively augments the matching while adjusting the potentials to preserve dual feasibility, ensuring reduced costs for all edges. Upon termination with a perfect matching, strong duality holds with zero duality gap, as the primal and dual objectives equalize. This optimality is certified by complementary slackness conditions: for every assigned edge with , the dual constraint is tight (), and for every dual variable with positive contribution, the corresponding primal constraint is tight. These conditions align the algorithm's potential adjustments with the economic interpretation of dual variables as shadow prices in assignment markets.[1][23] This duality perspective traces back to foundational work by John von Neumann, whose 1953 analysis of expanding economies in the second edition of Theory of Games and Economic Behavior (co-authored with Oskar Morgenstern) applied similar dual formulations to assignment-like problems in resource allocation, influencing the development of primal-dual algorithms for combinatorial optimization.[24]Matrix-Based Algorithm Description
Initial Reduction Steps
The initial reduction steps of the Hungarian algorithm prepare the cost matrix for the assignment problem, which involves assigning n workers to n tasks to minimize total cost, represented by an n × n matrix C where c_{ij} denotes the cost of assigning worker i to task j. These steps transform the original matrix into a reduced form with non-negative entries and zeros indicating potential optimal assignments, without altering the optimal solution set.[25] In the first step, row reduction is performed by identifying the minimum entry in each row of the cost matrix and subtracting this minimum value from every element in that row. This operation ensures that at least one zero appears in each row, as the smallest entry becomes zero while all others remain unchanged or increase relative to it. For example, consider a 3 × 3 cost matrix:| Task 1 | Task 2 | Task 3 | |
|---|---|---|---|
| Worker 1 | 5 | 9 | 2 |
| Worker 2 | 6 | 7 | 4 |
| Worker 3 | 8 | 3 | 1 |
| Task 1 | Task 2 | Task 3 | |
|---|---|---|---|
| Worker 1 | 3 | 7 | 0 |
| Worker 2 | 2 | 3 | 0 |
| Worker 3 | 7 | 2 | 0 |
| Task 1 | Task 2 | Task 3 | |
|---|---|---|---|
| Worker 1 | 1 | 5 | 0 |
| Worker 2 | 0 | 1 | 0 |
| Worker 3 | 5 | 0 | 0 |
