Hungarian algorithm: Difference between revisions

Revision as of 00:14, 24 November 2023

Latest revision as of 15:38, 23 May 2025

(28 intermediate revisions by 21 users not shown)

Line 2:

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 Hungarian mathematicians, [[Dénes Kőnig]] and [[Jenő Egerváry]].<ref name="kuhn1955">Harold W. Kuhn, "The Hungarian Method for the assignment problem", ''[[Naval Research Logistics Quarterly]]'', '''2''': 83–97, 1955. Kuhn's original publication.</ref><ref name="kuhn1956">Harold W. Kuhn, "Variants of the Hungarian method for assignment problems", ''Naval Research Logistics Quarterly'', '''3''': 253–258, 1956.</ref> 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.<ref>{{Cite web|url=http://www.lix.polytechnique.fr/~ollivier/JACOBI/presentationlEngl.htm|archive-url = https://web.archive.org/web/20151016182619/http://www.lix.polytechnique.fr/~ollivier/JACOBI/presentationlEngl.htm|archive-date = 16 October 2015|title = Presentation}}</ref>

[[James Munkres]] reviewed the algorithm in 1957 and observed that it is [[Time complexity#Strongly and weakly -polynomial time|(strongly) polynomial]].<ref name="munkres">J. Munkres, "Algorithms for the Assignment and Transportation Problems", ''[[Journal of the Society for Industrial and Applied Mathematics]]'', '''5'''(1):32–38, 1957 March.</ref> Since then the algorithm has been known also as the '''Kuhn–Munkres algorithm''' or '''Munkres assignment algorithm'''. The [[Computational complexity theory#Time and space complexity|time complexity]] of the original algorithm was <math>O(n^4)</math>, however [[Jack Edmonds|Edmonds]] and [[Richard Karp|Karp]], and independently Tomizawa, noticed that it can be modified to achieve an <math>O(n^3)</math> running time.<ref>{{Cite journal|last1=Edmonds|first1=Jack|last2=Karp|first2=Richard M.|date=1972-04-01|title=Theoretical Improvements in Algorithmic Efficiency for Network Flow Problems|journal=Journal of the ACM|volume=19|issue=2|pages=248–264|language=EN|doi=10.1145/321694.321699|s2cid=6375478|doi-access=free}}</ref><ref>{{Cite journal|last=Tomizawa|first=N.|date=1971|title=On some techniques useful for solution of transportation network problems|journal=Networks|language=en|volume=1|issue=2|pages=173–194|doi=10.1002/net.3230010206|issn=1097-0037}}</ref> One of the most popular{{Citation needed|date=November 2019}} <math>O(n^3)</math> variants is the Jonker–Volgenant algorithm.<ref name="JVAlg">{{cite journal |last1=Jonker |first1=R. |last2=Volgenant |first2=A. |title=A shortest augmenting path algorithm for dense and sparse linear assignment problems |journal=Computing |date=December 1987 |volume=38 |issue=4 |pages=325–340 |doi=10.1007/BF02278710|s2cid=7806079 }}</ref> [[L. R. Ford, Jr.|Ford]] and [[D. R. Fulkerson|Fulkerson]] extended the method to general maximum flow problems in form of the [[Ford–Fulkerson algorithm]].

==The problem==

Line 8:

=== Example ===

In this simple example, there are three workers: Alice, Bob and DoraCarol. 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 (mathematics)|matrix]] of the costs of the workers doing the jobs. For example:

:{| class="wikitable" style="text-align:center;"

! {{diagonal split header|Worker &nbsp; &nbsp;|Task}}

Line 25:

| '''$3'''

|-

! DoraCarol

| $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, DoraCarol sweep the floors, and Bob wash the windows. This can be confirmed using brute force:

:{| class="wikitable" style="text-align:center;"

! {{diagonal split header|Sweep|Clean}}

! Alice

! Bob

! DoraCarol

|-

! Alice

Line 48:

| $18

|-

! DoraCarol

| '''$15'''

| $16

Line 55:

=== Matrix formulation ===

In the matrix formulation, we are given a nonnegativean ''n''×''n'' [[Matrix (mathematics)|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 (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 nonnegative cost <math>c(i,j)</math>. We want to find a [[perfect matching]] with a minimum total cost.

==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: <math>\sum_{v\in S\cup T} y(v)</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|&Delta;}} 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|&Delta;}} on the vertices of <math>Z \cap S</math> and decrease {{mvar|y}} by {{math|&Delta;}} 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|&Delta;}} 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:

*/

#includeimport <cassert>;

#include <iostream>

import std;

#include <limits>

#include <vector>

template <typename T, typename U>

using namespace std;

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 <classtypename T>
constexpr bool ckmin(T & a, const T & b) {
return b < a ? a = b, 1true : 0false;
}

/**

* @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>

template <class T> vectorVector<T> hungarian(const vectorVector<vectorVector<T>>& &C) {

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

vectorVector<int> job(W + 1, -1);

vectorVector<T> ys(J),;
Vector<T> yt(W + 1); // potentials

// -yt[W] will equal the sum of all deltas

vectorVector<T> answers;

const T inf = numeric_limitsNumericLimits<T>::max();

for (int j_curjCur = 0; j_curjCur < J; ++j_curjCur) { // assign j_curjCur-th job

int w_curwCur = W;

job[w_curwCur] = j_curjCur;

// min reduced cost over edges from Z to worker w

vectorVector<T> min_tominTo(W + 1, inf);

vectorVector<int> prvprev(W + 1, -1); // previous worker on alternating path

vectorVector<bool> in_ZinZ(W + 1); // whether worker is in Z

while (job[w_curwCur] != -1) { // runs at most j_curjCur + 1 times

in_ZinZ[w_curwCur] = true;

const int j = job[w_curwCur];

T delta = inf;

int w_nextwNext;

for (int w = 0; w < W; ++w) {

if (!in_ZinZ[w]) {

if (ckmin(min_tominTo[w], C[j][w] - ys[j] - yt[w]))

prvprev[w] = w_curwCur;

if (ckmin(delta, min_tominTo[w])) w_next = w;

wNext = w;

}

}

// delta will always be non-negativenonnegative,

// except possibly during the first time this loop runs

// if any entries of C[j_curjCur] are negative

for (int w = 0; w <= W; ++w) {

if (in_ZinZ[w]) ys[job[w]] += delta, yt[w] -= delta;{

else min_to ys[job[w]] -+= delta;

yt[w] -= delta;

} else {

minTo[w] -= delta;

}

}

w_curwCur = w_nextwNext;

}

// update assignments along alternating path

for (int w; w_curwCur != W; w_curwCur = w) job[w_cur] = job[w = prv[w_cur]];

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: DoraCarol -> 4

* wash windows: Bob -> 3

*/

void sanity_check_hungariansanityCheckHungarian() {

vectorVector<vectorVector<int>> costs{{8, 5, 9}, {4, 2, 4}, {7, 3, 8}};

assert((hungarian(costs) == vectorVector<int>{5, 9, 15}));

cerr <<std::println(stderr, "Sanity check passed.\n");

}

/**

// * @brief solves https://open.kattis.com/problems/cordonbleu

void cordon_bleu() {

*/

int N, M;

void cordon_bleucordonBleu() {

cin >> N >> M;

vector<pair<int, int>> B(N), C(M);

int N, M;

vector<pair<int, int>> bottles(N), couriers(M);

for (auto &b std:: bottles) cin >> b.firstN >> b.secondM;

Vector<Pair<int, int>> B(N);

for (auto &c : couriers) cin >> c.first >> c.second;

pairVector<Pair<int, int>> restC(M);

vectorVector<pairPair<int, int>> bottles(N), couriers(M);

cin >> rest.first >> rest.second;

vectorVector<vectorPair<int>> costs(N, vector<int>> couriers(N + M - 1));

autofor dist(const = [&](pairPair<int, int>& x, pair<int, int>b: ybottles) {

returnstd::cin abs(x.first>> - yb.first) + abs(x.second ->> yb.second);

for (const Pair<int, int>& c: couriers)

for (auto &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

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);

}

}

cout << std::println(hungarian(costs).back() << "\n");

}

/**

* @brief Entry point into the program.

*

* @return The return code of the program.

*/

int main() {

sanity_check_hungariansanityCheckHungarian();

cordon_bleucordonBleu();

}

</syntaxhighlight>

Line 286 ⟶ 326:

<syntaxhighlight lang="cpp" line="1">

template <typename T>

template <class T> vectorVector<T> hungarian(const vectorVector<vectorVector<T>>& &C) {

const int J = (int)size(C), W = (int)size(C[0]);

const int J = (static_cast<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

vectorVector<int> job(W + 1, -1);

vectorVector<T> h(W); // Johnson potentials

vectorVector<T> answers;

T ans_curansCur = 0;

const T inf = numeric_limitsNumericLimits<T>::max();

// assign j_curjCur-th job using Dijkstra with potentials

for (int j_curjCur = 0; j_curjCur < J; ++j_curjCur) {

int w_curwCur = W; // unvisited worker with minimum distance

job[w_curwCur] = j_curjCur;

vectorVector<T> dist(W + 1, inf); // Johnson-reduced distances

dist[W] = 0;

vectorVector<bool> vis(W + 1); // whether visited yet

vectorVector<int> prvprev(W + 1, -1); // previous worker on shortest path

while (job[w_curwCur] != -1) { // Dijkstra step: pop min worker from heap

T min_distminDist = inf;

vis[w_curwCur] = true;

int w_nextwNext = -1; // next unvisited worker with minimum distance

// consider extending shortest path by w_curwCur -> job[w_curwCur] -> w

for (int w = 0; w < W; ++w) {

if (!vis[w]) {

// sum of reduced edge weights w_curwCur -> job[w_curwCur] -> w

T edge = C[job[w_curwCur]][w] - h[w];

if (w_curwCur != W) {

edge -= C[job[w_curwCur]][w_curwCur] - h[w_curwCur];

assert(edge >= 0); // consequence of Johnson potentials

}

if (ckmin(dist[w], dist[w_curwCur] + edge)) prv[w] = w_cur;

if (ckmin(min_dist, dist prev[w])) w_next = wwCur;

if (ckmin(minDist, dist[w]))

wNext = w;

}

}

w_curwCur = w_nextwNext;

}

for (int w = 0; w < W; ++w) { // update potentials

ckmin(dist[w], dist[w_curwCur]);

h[w] += dist[w];

}

ans_cur += h[w_curwCur];

for (int w; w_curwCur != W; w_curwCur = w) job[w_cur] = job[w = prv[w_cur]];

answers.push_back(ans_cur) job[wCur] = job[w = prev[wCur]];

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 non-negativenonnegative 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.

Line 416 ⟶ 467:

|}

Find a non-covered zero and prime it. (mark it with a [[Prime (symbol)|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.

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.

{{Confusing|section|date=May 2023}}

'''FindOtherwise, find the lowest uncovered value. Subtract this from every unmarkeduncovered 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===

Repeat steps 4–5 until an assignment is possible; this is when the minimum number of lines used to cover all the 0s is equal to min(number of people, number of assignments), assuming dummy variables (usually the max cost) are used to fill in when the number of people is greater than the number of assignments.

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 557 ⟶ 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]