Hungarian algorithm

Example 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: : 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) Matrix formulation 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, : \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 G=(S, T; E) with worker vertices () and job vertices (), and the edges () each have a cost c(i,j). We want to find a perfect matching with a minimum total cost. ==The algorithm in terms of bipartite graphs==

Let us call a function y: (S\cup T) \to \mathbb{R} a potential if y(i)+y(j) \leq c(i, j) for each i \in S, j \in T. The value of potential is the sum of the potential over all vertices: :\sum_{v\in S\cup T} y(v). The cost of each perfect matching is at least the value of each potential. This can be seen by first noticing that 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 ij is called tight for a potential if y(i)+y(j) = c(i, j). Let us denote the subgraph of tight edges by G_y. The cost of a perfect matching in G_y (if there is one) equals the value of . During the algorithm we maintain a potential and an orientation of G_y (denoted by \overrightarrow{G_y}) which has the property that the edges oriented from to form a matching . Initially, is 0 everywhere, and all edges are oriented from to (so is empty). In each step, either we modify 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 are tight. We are done if is a perfect matching. In a general step, let R_S \subseteq S and R_T \subseteq T be the vertices not covered by (so R_S consists of the vertices in with no incoming edge and R_T consists of the vertices in with no outgoing edge). Let be the set of vertices reachable in \overrightarrow{G_y} from R_S by a directed path. This can be computed by breadth-first search. If R_T \cap Z is nonempty, then reverse the orientation of all edges along a directed path in \overrightarrow{G_y} from R_S to R_T. Thus the size of the corresponding matching increases by 1. If R_T \cap Z is empty, then let :\Delta := \min \{c(i,j)-y(i)-y(j): i \in Z \cap S, j \in T \setminus Z\}. is well defined because at least one such edge ij 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 Z \cap S and T \setminus Z. Increase by on the vertices of Z \cap S and decrease by on the vertices of Z \cap T. The resulting is still a potential, and although the graph G_y changes, it still contains (see the next subsections). We orient the new edges from to . By the definition of the set of vertices reachable from R_S increases (note that the number of tight edges does not necessarily increase). If the vertex added to Z \cap T is unmatched (that is, it is also in R_T), then at the next iteration the graph will have an augmenting path. We repeat these steps until is a perfect matching, in which case it gives a minimum cost assignment. The running time of this version of the method is O(n^4): is augmented times, and in a phase where is unchanged, there are at most potential changes (since increases every time). The time sufficient for a potential change is O(n^2). 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: • is of maximum possible size. • G_y contains an augmenting path. • contains a loose-tailed path: a path from some vertex in R_S to a vertex in 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 Z \cap S, guaranteeing that is well defined. If is of maximum possible size, we are of course finished. Otherwise, by Berge's lemma, there must exist an augmenting path with respect to in the underlying graph . However, this path may not exist in G_y: Although every even-numbered edge in is tight by the definition of , odd-numbered edges may be loose and thus absent from G_y. One endpoint of is in R_S, the other in R_T; w.l.o.g., suppose it begins in R_S. If every edge on is tight, then it remains an augmenting path in G_y and we are done. Otherwise, let uv be the first loose edge on . If v \notin Z then we have found a loose-tailed path and we are done. Otherwise, is reachable from some other path of tight edges from a vertex in R_S. Let P_v be the subpath of beginning at and continuing to the end, and let P' be the path formed by traveling along until a vertex on P_v is reached, and then continuing to the end of P_v. Observe that P' is an augmenting path in with at least one fewer loose edge than . can be replaced with P' and this reasoning process iterated (formally, using induction on the number of loose edges) until either an augmenting path in G_y or a loose-tailed path in is found. Proof that adjusting the potential y leaves M unchanged To show that every edge in remains after adjusting , it suffices to show that for an arbitrary edge in , either both of its endpoints, or neither of them, are in . To this end let vu be an edge in from to . It is easy to see that if is in then must be too, since every edge in is tight. Now suppose, toward contradiction, that u \in Z but v \notin Z. itself cannot be in 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 R_S to . This path must avoid , since that is by assumption not in , so the vertex immediately preceding in this path is some other vertex v' \in T. v'u is a tight edge from to and is thus in . But then contains two edges that share the vertex , contradicting the fact that is a matching. Thus every edge in has either both endpoints or neither endpoint in . Proof that remains a potential To show that 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 by the preceding paragraph, so consider an arbitrary edge from to . If y(u) is increased by , then either v \in Z \cap T, in which case y(v) is decreased by , leaving the total potential of the edge unchanged, or v \in T \setminus Z, in which case the definition of guarantees that y(u)+y(v)+\Delta \leq c(u,v). Thus remains a potential. ==The algorithm in O(n3) time==

Suppose there are J jobs and W workers (J\le W). 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 jth job and update the total cost in time O(jW), yielding an overall time complexity of O\left(\sum_{j=1}^JjW\right)=O(J^2W). Note that this is better than O(W^3) when the number of jobs is small relative to the number of workers. 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 S_j denote the set of the first j jobs and T denote the set of all workers. Before the jth step of the algorithm, we assume that we have a matching on S_{j-1}\cup T that matches all jobs in S_{j-1} and potentials 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 jth step, we add the jth job to S_{j-1} to form S_j and initialize Z=\{j\}. At all times, every vertex in Z will be reachable from the jth job in G_y. While Z does not contain a worker that has not been assigned a job, let :\Delta := \min\{c(j, w)-y(j)-y(w):j\in Z\cap S_j, w\in T\setminus Z\} and w_{\text{next}} denote any 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 Z to w_{\text{next}}. • If w_{\text{next}} is unmatched, then we have an augmenting path in the subgraph of tight edges from j to w_{\text{next}}. After toggling the matching along this path, we have now matched the first j jobs, and this procedure terminates. • Otherwise, we add w_{\text{next} } and the job matched with it to Z. Adjusting potentials takes O(W) time. Recomputing \Delta and w_{\text{next}} after changing the potentials and Z also can be done in O(W) time. Case 1 can occur at most j-1 times before case 2 occurs and the procedure terminates, yielding the overall time complexity of O(jW). Implementation in C++ For convenience of implementation, the code below adds an additional worker w_{W} such that y(w_{W}) stores the negation of the sum of all \Delta computed so far. After the jth job is added and the matching updated, the cost of the current matching equals the sum of all \Delta computed so far, or -y(w_{W}). This code is adapted from e-maxx :: algo. /** * Solution to https://open.kattis.com/problems/cordonbleu using Hungarian * algorithm. */ import ; import std; template using Pair = std::pair; template using Vector = std::vector; template using NumericLimits = std::numeric_limits; /** * @brief Checks if b constexpr bool ckmin(T& a, const T& b) { return b Vector hungarian(const Vector>& C) { const int J = static_cast(C.size()); const int W = static_cast(C[0].size()); assert(J job(W + 1, -1); Vector ys(J); Vector yt(W + 1); // potentials // -yt[W] will equal the sum of all deltas Vector answers; const T inf = NumericLimits::max(); for (int jCur = 0; jCur minTo(W + 1, inf); Vector prev(W + 1, -1); // previous worker on alternating path Vector 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 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> costs{{8, 5, 9}, {4, 2, 4}, {7, 3, 8}}; assert((hungarian(costs) == Vector{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> B(N); Vector> C(M); Vector> bottles(N); Vector> couriers(M); for (const Pair& b: bottles) std::cin >> b.first >> b.second; for (const Pair& c: couriers) std::cin >> c.first >> c.second; Pair rest; std::cin >> rest.first >> rest.second; Vector> costs(N, Vector(N + M - 1)); auto dist = [&](const Pair& x, const Pair& y) -> int { return std::abs(x.first - y.first) + std::abs(x.second - y.second); }; for (int b = 0; b bottle -> restaurant costs[b][c] = dist(couriers[c], bottles[b]) + dist(bottles[b], rest); for (int _ = 0; _ 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==

The Hungarian algorithm can be seen to be equivalent to the successive shortest path algorithm for minimum cost flow, 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 h for workers 0\dots W-1 are equal to the potentials y from the previous solution up to a constant offset. When the graph is sparse (there are only M allowed job, worker pairs), it is possible to optimize this algorithm to run in O(JM+J^2\log W) time by using a Fibonacci heap to determine w_{\text{next}} instead of iterating over all W workers to find the one with minimum distance (alluded to here). template Vector hungarian(const Vector>& C) { const int J = static_cast(C.size()); const int W = static_cast(C[0].size()); assert(J job(W + 1, -1); Vector h(W); // Johnson potentials Vector answers; T ansCur = 0; const T inf = NumericLimits::max(); // assign jCur-th job using Dijkstra with potentials for (int jCur = 0; jCur dist(W + 1, inf); // Johnson-reduced distances dist[W] = 0; Vector vis(W + 1); // whether visited yet Vector 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 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 ==Matrix interpretation==

This variant of the algorithm follows the formulation given by Flood, and later described more explicitly by Munkres, who proved it runs in \mathcal{O}(n^4) time. the minimum number of lines (minimum vertex cover) will be (the size of maximum matching). Thus, when lines are required, minimum cost assignment can be found by looking at only zeroes in the matrix. ==Bibliography==