Example

In this simple example, there are three workers: Alice, Bob and Dora. 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:

Task Worker	Clean bathroom	Sweep floors	Wash windows
Alice	$8	$4	$7
Bob	$5	$2	$3
Dora	$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, Dora sweep the floors, and Bob wash the windows. This can be confirmed using brute force:

Clean Sweep	Alice	Bob	Dora
Alice	—	$17	$16
Bob	$18	—	$18
Dora	$15	$16	—

(the unassigned person washes the windows)

Matrix formulation

In the matrix formulation, we are given a nonnegative 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,

${\displaystyle \min _{P}\operatorname {Tr} (PC)\;,}$

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

The algorithm can equivalently be described by formulating the problem using a bipartite graph. We have a complete bipartite graph ${\displaystyle G=(S,T;E)}$ with n worker vertices (S) and n job vertices (T), and the edges (E) each have a nonnegative cost ${\displaystyle c(i,j)}$. We want to find a perfect matching with a minimum total cost.

Proof that the algorithm makes progress

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.
• ${\displaystyle G_{y}}$ contains an augmenting path.
• G contains a loose-tailed path: a path from some vertex in ${\displaystyle R_{S}}$ to a vertex in ${\displaystyle T\setminus Z}$ 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 ${\displaystyle Z\cap S}$, 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 ${\displaystyle G_{y}}$: Although every even-numbered edge in P is tight by the definition of M, odd-numbered edges may be loose and thus absent from ${\displaystyle G_{y}}$. One endpoint of P is in ${\displaystyle R_{S}}$, the other in ${\displaystyle R_{T}}$; w.l.o.g., suppose it begins in ${\displaystyle R_{S}}$. If every edge on P is tight, then it remains an augmenting path in ${\displaystyle G_{y}}$ and we are done. Otherwise, let ${\displaystyle uv}$ be the first loose edge on P. If ${\displaystyle v\notin Z}$ 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 ${\displaystyle R_{S}}$. Let ${\displaystyle P_{v}}$ be the subpath of P beginning at v and continuing to the end, and let ${\displaystyle P'}$ be the path formed by traveling along Q until a vertex on ${\displaystyle P_{v}}$ is reached, and then continuing to the end of ${\displaystyle P_{v}}$. Observe that ${\displaystyle P'}$ is an augmenting path in G with at least one fewer loose edge than P. P can be replaced with ${\displaystyle P'}$ and this reasoning process iterated (formally, using induction on the number of loose edges) until either an augmenting path in ${\displaystyle G_{y}}$ or a loose-tailed path in G is found.

Proof that adjusting the potential y leaves M unchanged

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 ${\displaystyle vu}$ 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 ${\displaystyle u\in Z}$ but ${\displaystyle v\notin Z}$. u itself cannot be in ${\displaystyle R_{S}}$ because it is the endpoint of a matched edge, so there must be some directed path of tight edges from a vertex in ${\displaystyle R_{S}}$ 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 ${\displaystyle v'\in T}$. ${\displaystyle v'u}$ 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

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 ${\displaystyle y(u)}$ is increased by Δ, then either ${\displaystyle v\in Z\cap T}$, in which case ${\displaystyle y(v)}$ is decreased by Δ, leaving the total potential of the edge unchanged, or ${\displaystyle v\in T\setminus Z}$, in which case the definition of Δ guarantees that ${\displaystyle y(u)+y(v)+\Delta \leq c(u,v)}$. Thus y remains a potential.

Adding the j-th job in O(jW) time

We use the same notation as the previous section, though we modify their definitions as necessary. Let ${\displaystyle S_{j}}$ denote the set of the first ${\displaystyle j}$ jobs and ${\displaystyle T}$ denote the set of all workers.

Before the ${\displaystyle j}$th step of the algorithm, we assume that we have a matching on ${\displaystyle S_{j-1}\cup T}$ that matches all jobs in ${\displaystyle S_{j-1}}$ and potentials ${\displaystyle y}$ 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 ${\displaystyle j}$th step, we add the ${\displaystyle j}$th job to ${\displaystyle S_{j-1}}$ to form ${\displaystyle S_{j}}$ and initialize ${\displaystyle Z=\{j\}}$. At all times, every vertex in ${\displaystyle Z}$ will be reachable from the ${\displaystyle j}$th job in ${\displaystyle G_{y}}$. While ${\displaystyle Z}$ does not contain a worker that has not been assigned a job, let

${\displaystyle \Delta$ :=\min\{c(j,w)-y(j)-y(w):j\in Z\cap S_{j},w\in T\setminus Z\}}

and ${\displaystyle w_{\text{next}}}$ denote any ${\displaystyle w}$ 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 ${\displaystyle Z}$ to ${\displaystyle w_{\text{next}}}$.

• If ${\displaystyle w_{\text{next}}}$ is unmatched, then we have an augmenting path in the subgraph of tight edges from ${\displaystyle j}$ to ${\displaystyle w_{\text{next}}}$. After toggling the matching along this path, we have now matched the first ${\displaystyle j}$ jobs, and this procedure terminates.
• Otherwise, we add ${\displaystyle w_{\text{next}}}$ and the job matched with it to ${\displaystyle Z}$.

Adjusting potentials takes ${\displaystyle O(W)}$ time. Recomputing ${\displaystyle \Delta }$ and ${\displaystyle w_{\text{next}}}$ after changing the potentials and ${\displaystyle Z}$ also can be done in ${\displaystyle O(W)}$ time. Case 1 can occur at most ${\displaystyle j-1}$ times before case 2 occurs and the procedure terminates, yielding the overall time complexity of ${\displaystyle O(jW)}$.

Implementation in C++

For convenience of implementation, the code below adds an additional worker ${\displaystyle w_{W}}$ such that ${\displaystyle y(w_{W})}$ stores the negation of the sum of all ${\displaystyle \Delta }$ computed so far. After the ${\displaystyle j}$th job is added and the matching updated, the cost of the current matching equals the sum of all ${\displaystyle \Delta }$ computed so far, or ${\displaystyle -y(w_{W})}$.

This code is adapted from e-maxx :: algo.^[8]

/**
 * Solution to https://open.kattis.com/problems/cordonbleu using Hungarian
 * algorithm.
 */

#include <cassert>
#include <iostream>
#include <limits>
#include <vector>
using namespace std;

/**
 * Sets a = min(a, b)
 * @return true if b < a
 */
template <class T> bool ckmin(T &a, const T &b) { return b < a ? a = b, 1 : 0; }

/**
 * 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 <class T> vector<T> hungarian(const vector<vector<T>> &C) {
    const int J = (int)size(C), W = (int)size(C[0]);
    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), yt(W + 1);  // potentials
    // -yt[W] will equal the sum of all deltas
    vector<T> answers;
    const T inf = numeric_limits<T>::max();
    for (int j_cur = 0; j_cur < J; ++j_cur) {  // assign j_cur-th job
        int w_cur = W;
        job[w_cur] = j_cur;
        // min reduced cost over edges from Z to worker w
        vector<T> min_to(W + 1, inf);
        vector<int> prv(W + 1, -1);  // previous worker on alternating path
        vector<bool> in_Z(W + 1);    // whether worker is in Z
        while (job[w_cur] != -1) {   // runs at most j_cur + 1 times
            in_Z[w_cur] = true;
            const int j = job[w_cur];
            T delta = inf;
            int w_next;
            for (int w = 0; w < W; ++w) {
                if (!in_Z[w]) {
                    if (ckmin(min_to[w], C[j][w] - ys[j] - yt[w]))
                        prv[w] = w_cur;
                    if (ckmin(delta, min_to[w])) w_next = w;
                }
            }
            // delta will always be non-negative,
            // except possibly during the first time this loop runs
            // if any entries of C[j_cur] are negative
            for (int w = 0; w <= W; ++w) {
                if (in_Z[w]) ys[job[w]] += delta, yt[w] -= delta;
                else min_to[w] -= delta;
            }
            w_cur = w_next;
        }
        // update assignments along alternating path
        for (int w; w_cur != W; w_cur = w) job[w_cur] = job[w = prv[w_cur]];
        answers.push_back(-yt[W]);
    }
    return answers;
}

/**
 * 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: Dora -> 4
 *   wash windows: Bob -> 3
 */
void sanity_check_hungarian() {
    vector<vector<int>> costs{{8, 5, 9}, {4, 2, 4}, {7, 3, 8}};
    assert((hungarian(costs) == vector<int>{5, 9, 15}));
    cerr << "Sanity check passed.\n";
}

// solves https://open.kattis.com/problems/cordonbleu
void cordon_bleu() {
    int N, M;
    cin >> N >> M;
    vector<pair<int, int>> B(N), C(M);
    vector<pair<int, int>> bottles(N), couriers(M);
    for (auto &b : bottles) cin >> b.first >> b.second;
    for (auto &c : couriers) cin >> c.first >> c.second;
    pair<int, int> rest;
    cin >> rest.first >> rest.second;
    vector<vector<int>> costs(N, vector<int>(N + M - 1));
    auto dist = [&](pair<int, int> x, pair<int, int> y) {
        return abs(x.first - y.first) + 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);
        }
    }
    cout << hungarian(costs).back() << "\n";
}

int main() {
    sanity_check_hungarian();
    cordon_bleu();
}

Step 1

For each row, its minimum element is subtracted from every element in that row. This causes all elements to have non-negative 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

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

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

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. (If 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:
  1. Substep 1: Find a starred zero on the corresponding column. If there is one, go to Substep 2, else, stop.
  2. 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

Find the lowest uncovered value. Subtract this from every unmarked element and add it to every element covered by two lines.

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.

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,^[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.

Implementations

Note that not all of these satisfy the ${\displaystyle O(n^{3})}$ time complexity, even if they claim so. Some may contain errors, implement the slower ${\displaystyle O(n^{4})}$ 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 ${\displaystyle O(n^{3})}$ time complexity)
• Julia implementation
• C implementation claiming ${\displaystyle O(n^{3})}$ time complexity
• Java implementation claiming ${\displaystyle O(n^{3})}$ time complexity
• Python implementation
• Ruby implementation with unit tests
• C# implementation claiming ${\displaystyle O(n^{3})}$ time complexity
• D implementation with unit tests (port of a Java version claiming ${\displaystyle O(n^{3})}$)
• Online interactive implementation
• Serial and parallel implementations.
• Matlab and C
• Perl implementation
• C++ implementation
• C++ implementation claiming ${\displaystyle O(n^{3})}$ time complexity (BSD style open source licensed)
• MATLAB implementation
• C implementation
• JavaScript implementation with unit tests (port of a Java version claiming ${\displaystyle O(n^{3})}$ time complexity)
• Clue R package proposes an implementation, solve_LSAP
• Node.js implementation on GitHub
• Python implementation in scipy package