TheHungarian method is acombinatorial optimizationalgorithm that solves theassignment problem inpolynomial time and which anticipated laterprimal–dual methods. It was developed and published in 1955 byHarold 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 andJenő Egerváry.[1][2] However, in 2006 it was discovered thatCarl 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 theKuhn–Munkres algorithm orMunkres assignment algorithm. Thetime complexity of the original algorithm was, howeverEdmonds andKarp, and independently Tomizawa, noticed that it can be modified to achieve an running time.[5][6]Ford andFulkerson extended the method to general maximum flow problems in form of theFord–Fulkerson algorithm.
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 amatrix 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 |
| 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:
Clean Sweep | Alice | Bob | Carol |
|---|---|---|---|
| Alice | — | $17 | $16 |
| Bob | $18 | — | $18 |
| Carol | $15 | $16 | — |
In the matrix formulation, we are given ann×nmatrix, where the element in thei-th row andj-th column represents the cost of assigning thej-th job to thei-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 matrixC to minimize thetrace of a matrix,
whereP is apermutation matrix. (Equivalently, the columns can be permuted usingCP.)
If the goal is to find the assignment that yields themaximum cost, the problem can be solved by negating the cost matrixC.
The algorithm can equivalently be described by formulating the problem using a bipartite graph. We have acomplete bipartite graph withn worker vertices (S) andn job vertices (T), and the edges (E) each have a cost. We want to find aperfect matching with a minimum total cost.
Let us call a function apotential if for each.
Thevalue of potentialy 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 oftight edges: an edge is called tight for a potentialy if. Let us denote thesubgraph oftight edges by. The cost of a perfect matching in (if there is one) equals the value ofy.
During the algorithm we maintain a potentialy and anorientation of (denoted by) which has the property that the edges oriented fromT toS form a matchingM. Initially,y is 0 everywhere, and all edges are oriented fromS toT (soM is empty). In each step, either we modifyy so that its value increases, or modify the orientation to obtain a matching with more edges. We maintain the invariant that all the edges ofM are tight. We are done ifM is a perfect matching.
In a general step, let and be the vertices not covered byM (so consists of the vertices inS with no incoming edge and consists of the vertices inT with no outgoing edge). LetZ be the set of vertices reachable in from by a directed path. This can be computed bybreadth-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. Increasey byΔ on the vertices of and decreasey byΔ on the vertices of. The resultingy is still a potential, and although the graph changes, it still containsM (see the next subsections). We orient the new edges fromS toT. By the definition ofΔ the setZ 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 untilM 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 augmentedn times, and in a phase whereM is unchanged, there are at mostn potential changes (sinceZ increases every time). The time sufficient for a potential change is.
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:
IfM is of maximum possible size, we are of course finished. Otherwise, byBerge's lemma, there must exist an augmenting pathP with respect toM in the underlying graphG. However, this path may not exist in: Although every even-numbered edge inP is tight by the definition ofM, odd-numbered edges may be loose and thus absent from. One endpoint ofP is in, the other in; w.l.o.g., suppose it begins in. If every edge onP is tight, then it remains an augmenting path in and we are done. Otherwise, let be the first loose edge onP. If then we have found a loose-tailed path and we are done. Otherwise,v is reachable from some other pathQ of tight edges from a vertex in. Let be the subpath ofP beginning atv and continuing to the end, and let be the path formed by traveling alongQ until a vertex on is reached, and then continuing to the end of. Observe that is an augmenting path inG with at least one fewer loose edge thanP.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 inG is found.
To show that every edge inM remains after adjustingy, it suffices to show that for an arbitrary edge inM, either both of its endpoints, or neither of them, are inZ. To this end let be an edge inM fromT toS. It is easy to see that ifv is inZ thenu must be too, since every edge inM 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 tou. This path must avoidv, since that is by assumption not inZ, so the vertex immediately precedingu in this path is some other vertex. is a tight edge fromT toS and is thus inM. But thenM contains two edges that share the vertexu, contradicting the fact thatM is a matching. Thus every edge inM has either both endpoints or neither endpoint inZ.
To show thaty 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 inM by the preceding paragraph, so consider an arbitrary edgeuv fromS toT. 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. Thusy remains a potential.
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 theth 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.
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 theth 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 theth step, we add theth job to to form and initialize. At all times, every vertex in will be reachable from theth 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.
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.
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 theth 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>;importstd;template<typenameT,typenameU>usingPair=std::pair<T,U>;template<typenameT>usingVector=std::vector<T>;template<typenameT>usingNumericLimits=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<typenameT>constexprboolckmin(T&a,constT&b){returnb<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<typenameT>Vector<T>hungarian(constVector<Vector<T>>&C){constintJ=static_cast<int>(C.size());constintW=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 convenienceVector<int>job(W+1,-1);Vector<T>ys(J);Vector<T>yt(W+1);// potentials// -yt[W] will equal the sum of all deltasVector<T>answers;constTinf=NumericLimits<T>::max();for(intjCur=0;jCur<J;++jCur){// assign jCur-th jobintwCur=W;job[wCur]=jCur;// min reduced cost over edges from Z to worker wVector<T>minTo(W+1,inf);Vector<int>prev(W+1,-1);// previous worker on alternating pathVector<bool>inZ(W+1);// whether worker is in Zwhile(job[wCur]!=-1){// runs at most jCur + 1 timesinZ[wCur]=true;constintj=job[wCur];Tdelta=inf;intwNext;for(intw=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 negativefor(intw=0;w<=W;++w){if(inZ[w]){ys[job[w]]+=delta;yt[w]-=delta;}else{minTo[w]-=delta;}}wCur=wNext;}// update assignments along alternating pathfor(intw;wCur!=W;wCur=w)job[wCur]=job[w=prev[wCur]];answers.push_back(-yt[W]);}returnanswers;}/** * @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 */voidsanityCheckHungarian(){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 */voidcordonBleu(){intN;intM;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(constPair<int,int>&b:bottles)std::cin>>b.first>>b.second;for(constPair<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));autodist=[&](constPair<int,int>&x,constPair<int,int>&y)->int{returnstd::abs(x.first-y.first)+std::abs(x.second-y.second);};for(intb=0;b<N;++b){for(intc=0;c<M;++c)// courier -> bottle -> restaurantcosts[b][c]=dist(couriers[c],bottles[b])+dist(bottles[b],rest);for(int_=0;_<N-1;++_)// restaurant -> bottle -> restaurantcosts[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. */intmain(){sanityCheckHungarian();cordonBleu();}
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 fromJohnson'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 possibleto optimize this algorithm to run in time by using aFibonacci heap to determine instead of iterating over all workers to find the one with minimum distance (alluded tohere).
template<typenameT>Vector<T>hungarian(constVector<Vector<T>>&C){constintJ=static_cast<int>(C.size());constintW=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 convenienceVector<int>job(W+1,-1);Vector<T>h(W);// Johnson potentialsVector<T>answers;TansCur=0;constTinf=NumericLimits<T>::max();// assign jCur-th job using Dijkstra with potentialsfor(intjCur=0;jCur<J;++jCur){intwCur=W;// unvisited worker with minimum distancejob[wCur]=jCur;Vector<T>dist(W+1,inf);// Johnson-reduced distancesdist[W]=0;Vector<bool>vis(W+1);// whether visited yetVector<int>prev(W+1,-1);// previous worker on shortest pathwhile(job[wCur]!=-1){// Dijkstra step: pop min worker from heapTminDist=inf;vis[wCur]=true;intwNext=-1;// next unvisited worker with minimum distance// consider extending shortest path by wCur -> job[wCur] -> wfor(intw=0;w<W;++w){if(!vis[w]){// sum of reduced edge weights wCur -> job[wCur] -> wTedge=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(intw=0;w<W;++w){// update potentialsckmin(dist[w],dist[wCur]);h[w]+=dist[w];}ans_cur+=h[wCur];for(intw;wCur!=W;wCur=w)job[wCur]=job[w=prev[wCur]];answers.push_back(ansCur);}returnanswers;}
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.
Givenn workers and tasks, the problem is written in the form of ann×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.
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 naivegreedy 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 aren! 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.
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.
All zeros in the matrix must be covered by marking as few rows and/or columns as possible. Steps 3 and 4 formone 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 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 |
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 aprime symbol). If no such zero can be found, meaning all zeroes are covered, skip to step 5.
| × | ||||
| 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 |
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 |
| 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 |
All zeros are now covered with a minimal number of rows and columns.
The aforementioned detailed description isjust one way to draw the minimum number of lines to cover all the 0s. Other methods work as well.
If the number of starred zeros isn (or in the general case, wheren is the number of people andm is the number of jobs), the algorithm terminates. See theResult 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.
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 ben (the size of maximum matching[13]). Thus, whenn lines are required, minimum cost assignment can be found by looking at only zeroes in the matrix.
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.