The Complexity Of The Hungarian Algorithm For Assignment Problem: O(N^3) Implementation Details

Title: The Complexity Of The Hungarian Algorithm For Assignment Problem: O(N^3) Implementation Details

Introduction

Imagine you are the operations manager for a disaster relief organization. A devastating earthquake has just struck a major metropolitan area, and you have a fleet of 500 emergency vehicles—ambulances, fire engines, and rescue trucks—scattered across the city. Simultaneously, you have 500 distinct distress calls, each requiring a specific class of vehicle and located at a unique coordinate. You have exactly one second to decide which vehicle goes to which call. The cost? In this context, “cost” isn’t measured in dollars; it’s measured in human lives. Your goal is to minimize the total response time across the entire fleet, ensuring that the most critical calls are serviced as quickly as the logistics allow.

This is the Assignment Problem in its most urgent, real-world form. On the surface, it seems deceptively simple: match n agents to n tasks, one-to-one, with a given cost matrix, to minimize total cost. But this simplicity is a lie. The combinatorial explosion is immediate. For just 500 workers and 500 jobs, the brute-force approach—checking every possible permutation—would require evaluating roughly 500! (factorial) possibilities. That number is larger than the number of atoms in the observable universe by a margin that defies comprehension. Even the most powerful supercomputer on Earth would require more time than the lifespan of the star it orbits to find an answer.

We need strategy. We need structure. We need a polynomial-time solution.

Enter the Hungarian Algorithm. Originally developed by Harold Kuhn in 1955 and refined by James Munkres, this algorithm is a masterpiece of combinatorial optimization. It solves the assignment problem in polynomial time, but here is the critical nuance that separates a theoretical academic exercise from a piece of production-grade software: the naive implementation runs in O(N^4) time, while the optimized version—the one that actually gets deployed in satellite scheduling, supply chain logistics, and autonomous vehicle coordination—runs in O(N^3). Understanding that difference, and implementing it correctly, can mean the difference between a program that chokes on 500×500 matrices and one that handles 5000×5000 with ease.

In this blog post, we will dissect the Hungarian Algorithm at every level: from the combinatorial intuition behind Kuhn’s original work, to the linear programming duality that underpins its correctness, to the nitty-gritty details of the O(N^3) implementation. We’ll walk through a full code example in Python, analyze its complexity in rigorous detail, and explore real-world applications where milliseconds matter. By the end, you will not only understand why the Hungarian algorithm is O(N^3), but you will be able to implement it yourself, debug it, and optimize it for your own high-performance systems.

1. The Assignment Problem: Definitions and Hardness

1.1 Formal Definition

Given a square cost matrix C of size (n \times n), where entry (C_{ij}) represents the cost of assigning agent (i) to task (j), the assignment problem asks for a permutation (\pi) of ({1, 2, \dots, n}) that minimizes

[ \sum*{i=1}^n C*{i, \pi(i)}. ]

Each agent gets exactly one task, and each task gets exactly one agent. This is a classic combinatorial optimization problem, and it appears under many guises: matching rows to columns in a bipartite graph, optimal transport in the discrete Monge–Kantorovich sense, and the linear sum assignment problem (LSAP).

1.2 Why Brute Force is Impossible

For (n = 500), the number of permutations is (500!). Approximating this using Stirling’s formula:

[ n! \approx \sqrt{2\pi n} \left( \frac{n}{e} \right)^n. ]

Plugging in (n=500):

[ 500! \approx \sqrt{1000\pi} \left( \frac{500}{e} \right)^{500} \approx 1.22 \times 10^{1134}. ]

That is a 1 followed by 1134 zeros. Compare this to the number of atoms in the observable universe, which is “only” about (10^{80}). Even a quantum computer running a hypothetical perfect Grover’s search (quadratic speedup) would need to evaluate (\sqrt{500!} \approx 10^{567}) possibilities—still impossible.

The Hungarian algorithm reduces this to a clean (O(n^3)) using a combination of linear programming duality and clever matrix manipulation. To understand why it works, we need a small detour into history.

2. Historical Background: Kuhn, Munkres, and the Dual Problem

2.1 Kuhn’s Original 1955 Paper

Harold Kuhn was inspired by the work of two Hungarian mathematicians, Dénes Kőnig and Jenő Egerváry, on the theory of bipartite matching and the assignment problem. In fact, Kuhn named the algorithm “Hungarian” in their honor. He published “The Hungarian method for the assignment problem” in Naval Research Logistics Quarterly in 1955. The key insight was to turn the primal assignment problem into a dual problem that could be solved iteratively.

2.2 The Dual Linear Program

Every assignment problem can be expressed as a linear program:

[ \begin{aligned} \text{Minimize} \quad & \sum*{i=1}^n \sum*{j=1}^n C*{ij} x*{ij} \ \text{subject to} \quad & \sum*{j=1}^n x*{ij} = 1, \quad i = 1,\dots,n \ & \sum*{i=1}^n x*{ij} = 1, \quad j = 1,\dots,n \ & x_{ij} \in {0,1}. \end{aligned} ]

Relaxing the integer constraint to (x_{ij} \geq 0) still gives an integer optimum because the constraint matrix is totally unimodular. The dual problem introduces variables (u_i) (for each row) and (v_j) (for each column) such that:

[ \max*{u,v} \quad \sum*{i=1}^n ui + \sum{j=1}^n vj ] subject to [ u_i + v_j \leq C{ij} \quad \forall i,j. ]

At optimality, the primal and dual objective values coincide (strong duality). The Hungarian algorithm systematically builds dual feasible (ui, v_j) while maintaining a complementary slackness condition: (x{ij} = 1) only if (ui + v_j = C{ij}). This is exactly the condition that a zero cost edge in the reduced cost matrix (\tilde{C}{ij} = C{ij} - u_i - v_j) can be part of an optimal assignment.

2.3 Munkres’ Refinement: The O(N^3) Version

James Munkres, in 1957, published a note showing that the Hungarian algorithm could be implemented in (O(n^3)) time by careful bookkeeping of “covered” rows and columns, and by avoiding repeated scanning of the entire matrix. The key was to treat the process of finding augmenting paths (which is essentially a shortest augmenting path in a bipartite graph) as a single Dijkstra-like pass, rather than repeated BFS from each uncovered zero. This is the version we will implement.

3. The Core Idea of the Hungarian Algorithm

3.1 Primal-Dual Relationship

The algorithm maintains dual variables (ui) (row potentials) and (v_j) (column potentials) such that (u_i + v_j \leq C{ij}). The reduced cost matrix (\tilde{C}{ij} = C{ij} - ui - v_j) has nonnegative entries. The idea is to start with some feasible dual variables (e.g., (u_i = \min_j C{ij}), (v_j = 0)) so that each row has at least one zero, then to find a complete set of assignments that use only zero entries of the reduced matrix.

If such a complete assignment exists, complementary slackness says it is optimal for both primal and dual. If not, we adjust the dual variables to create new zeros while preserving feasibility, and repeat.

3.2 The Zero Graph

Define a bipartite graph where left nodes are agents, right nodes are tasks, and an edge exists if (\tilde{C}_{ij} = 0). The algorithm seeks a perfect matching in this graph. If none exists, we add more edges by lowering some (u_i) and raising some (v_j) so that new zeros appear.

3.3 Covering Zeros with Lines

Kőnig’s theorem (a precursor to the Hungarian method) states that in a bipartite graph, the size of the maximum matching equals the minimum number of lines (rows or columns) needed to cover all zeros. The Hungarian algorithm uses this theorem to decide when the zero graph does not have a perfect matching: if the minimum cover (in terms of rows and columns) is less than (n), then we cannot yet have a perfect assignment. We then adjust dual variables to reduce the cover count.

4. The Hungarian Algorithm: Detailed Step-by-Step (O(N^3) Version)

We will describe the algorithm as a sequence of steps. This is the implementation you will find in efficient libraries like scipy.optimize.linear_sum_assignment (which uses the Jonker-Volgenant algorithm, a variant of Hungarian, also O(N^3)). We follow Munkres’ original description adapted for clarity.

4.1 Data Structures

We maintain:

• cost: original matrix (will be modified as reduced costs)
• u and v: dual variables for rows and columns (array of length n)
• p and way: arrays used during the augmenting path search (predecessors)
• minv: for each column, the minimum value encountered during path search (a standard Dijkstra-like approach)
• used: boolean array to mark columns already visited
• assignment: array of length n where assignment[j] = i if column j is assigned to row i; by default -1.

4.2 Step 0: Initialization

Compute initial row potentials: for each row (i), set (ui = \min_j C{ij}). That guarantees that every row has at least one zero after subtraction. Set all (vj = 0). Then compute the reduced cost matrix (\tilde{C}{ij} = C_{ij} - u_i - v_j).

At this point, start with an empty assignment.

4.3 Step 1: For each row, find an augmenting path (Hungarian expansion)

This is the heart of the O(N^3) algorithm. Instead of solving a matching from scratch every time, we assign rows one by one. For each unassigned row i0:

1. Initialize: Set minv[j] = INF for all columns. Mark all columns as not visited. Set prev_col to a sentinel (e.g., -1). The current row we are trying to match is i0.

2. Iterate:

   • For each unvisited column j, compute the gap delta = C[i][j] - u[i] - v[j]. Since the reduced cost is non-negative, this is the amount we would need to reduce the reduced cost to zero.
   • If delta < minv[j], update minv[j] = delta and set way[j] = prev_col (the column from which we arrived).
   • Find the column j0 with the smallest minv[j] among unvisited columns. This is the column that can be “cheapest” to add to the zero graph.
   • Mark column j0 as visited.
   • If column j0 is currently unassigned (i.e., assignment[j0] == -1), then we have found an augmenting path. Stop.
   • Otherwise, let i1 = assignment[j0] be the row that is currently matched to j0. Set the current row to i1 and repeat.

3. Augmenting the assignment:

   • Once we have found an unassigned column, we trace back through the way array to adjust the dual variables and update the assignment. This is similar to updating potentials in Dijkstra: we need to decrease u for the rows along the alternating path and increase v for the columns visited.

   The standard way to update duals after the path is found:

   • For each visited column j, set v[j] += delta_min, where delta_min is the smallest minv among columns (the final column’s minv). Actually the textbook update: u[i] for rows along the alternating path are decreased by delta_min, and v[j] for columns in the alternating tree are increased by delta_min (or similar). The precise update ensures that all previously zero entries in the alternating tree remain zero, and the new entry (the column we added) becomes zero.

   A cleaner way to implement, following the scipy.optimize source code’s approach: we already compute delta values for each column; after the path is found, we update u[i] = u[i] + delta_min for all rows visited (or something). I will present the standard algorithm from Kuhn and Munkres as given in “Algorithm 1” of many references. Let me define the exact steps used in most O(N^3) implementations.

   Standard pseudocode (from Munkres 1957, as implemented in munkres package):

   For each row i:

   j = 0  # column index
   minv = [INF]*n
   visited = [False]*n
   way = [-1]*n
   
   done = False
   while not done:
       visited[j] = True
       i0 = assignment[j]  # row currently matched to column j (may be -1)
       delta = INF
       for j1 in range(n):
           if not visited[j1]:
               cur = C[i0][j1] - u[i0] - v[j1]   # reduced cost
               if cur < minv[j1]:
                   minv[j1] = cur
                   way[j1] = j
               if minv[j1] < delta:
                   delta = minv[j1]
                   j = j1
       # Now delta is the minimum minv among unvisited columns
       for j1 in range(n):
           if visited[j1]:
               u[assignment[j1]] += delta
               v[j1] -= delta
           else:
               minv[j1] -= delta
       if assignment[j] == -1:
           # column j is free, we have an augmenting path
           done = True
       # else loop continues with new current column j
   

   I know this looks cryptic. I will explain the logic in detail in the implementation section.

4.4 Step 2: Repeating for all rows

After processing all rows, we have a perfect assignment. The algorithm finishes.

5. Implementation in Python: A Complete O(N^3) Solver

Let’s implement the exact algorithm described above, which is the classic Munkres version used in many production systems.

5.1 Code

import numpy as np

def hungarian(cost_matrix):
    """
    Solve the assignment problem using the Hungarian algorithm (O(n^3)).
    cost_matrix: 2D array-like, shape (n, n). Must be square.
    Returns: row_indices, col_indices such that cost_matrix[row_indices, col_indices] is minimal.
    """
    n = cost_matrix.shape[0]
    # Convert to float for INF and potential
    cost = cost_matrix.astype(float)

    # Dual variables
    u = np.zeros(n, dtype=float)
    v = np.zeros(n, dtype=float)
    # assignment: for each column j, the row assigned, or -1 if unassigned
    p = np.full(n, -1, dtype=int)   # p[j] = row assigned to column j
    # way array for backtracking
    way = np.empty(n, dtype=int)
    # minv array for each column
    minv = np.empty(n, dtype=float)

    for i in range(n):
        # For current row i, we will try to match it
        p[0] = i   # temporary: use column 0 as start point (sentinel)
        j0 = 0      # current column
        minv[:] = np.inf
        used = np.zeros(n, dtype=bool)   # visited columns
        # main loop for this row
        while True:
            used[j0] = True
            i0 = p[j0]   # row currently matched to column j0 (may be initial i)
            delta = np.inf
            j1 = 0
            for j in range(n):
                if not used[j]:
                    cur = cost[i0, j] - u[i0] - v[j]
                    if cur < minv[j]:
                        minv[j] = cur
                        way[j] = j0
                    if minv[j] < delta:
                        delta = minv[j]
                        j1 = j
            # Update dual variables
            for j in range(n):
                if used[j]:
                    u[p[j]] += delta
                    v[j] -= delta
                else:
                    minv[j] -= delta
            j0 = j1   # next column to explore
            if p[j0] == -1:
                break   # found free column

        # Augment: update assignment along the path using way array
        while True:
            j1 = way[j0]
            p[j0] = p[j1]
            j0 = j1
            if j0 == 0:
                break
        # Note: the path is from column j0 (the free column) back to column 0.
        # p[0] remains the start row i, but after augmentation the assignments shift.

    # Build final row->col assignment from p
    # p[j] gives row assigned to column j. To return row->col indices:
    row_indices = np.full(n, -1, dtype=int)
    col_indices = np.arange(n)
    for j in range(n):
        row = p[j]
        if row != -1:
            row_indices[row] = j
    # Return sorted order (optional)
    return row_indices[np.argsort(row_indices)], col_indices

5.2 Explanation of the Algorithm Step by Step

Let’s dissect the loop for a single row i.

Initialization:
We set p[0] = i (a sentinel). This means we consider a dummy column 0 that is assigned to row i. This dummy column is just a starting point. We set j0 = 0 as the current column.

Main loop:

1. Mark column j0 as visited.
2. Let i0 = p[j0] be the row currently occupying column j0. Initially, i0 = i.
3. For all unvisited columns j, compute the reduced cost cur = cost[i0][j] - u[i0] - v[j]. Keep track of the minimum cur (call it delta) and which column j1 gives that minimum. Also store for each column the predecessor column way[j] (the column from which we came) if we update its minv value.
4. After scanning all columns, we have found the minimum delta. Now we update dual variables:
   • For each visited column j, we add delta to u[p[j]] (the row assigned to that column) and subtract delta from v[j]. This ensures that all currently used zero edges remain zero.
   • For unvisited columns, we subtract delta from minv[j] (effectively shifting the potentials so that the new zero edge emerges).
5. Set j0 = j1 (the column with the smallest minv). If p[j0] == -1, that column is free, so we have found an augmenting path. Otherwise, we continue the loop with the new current column.

Augmentation:

After breaking out of the loop, we have a path from the current free column j0 back to the starting sentinel column 0 via the way array. We walk backwards: for every step from j0 to j1 = way[j0], we reassign p[j0] = p[j1]. This essentially rotates the assignments along the alternating path, matching the new row i with column j0 (which becomes assigned), and free up columns that were previously matched.

The pseudocode might seem magical. To understand why this works, consider the algorithm as a variant of the shortest augmenting path in a graph with node potentials. The minv array is like a distance label, and the dual updates correspond to decreasing the potentials of rows in the alternating tree. The entire process is analogous to running the Hungarian algorithm using Dijkstra’s algorithm on the residual graph.

5.3 Test with a Small Example

Let’s test with a 3x3 matrix:

cost = np.array([[1, 2, 3],
                 [2, 4, 6],
                 [3, 6, 9]])
row_ind, col_ind = hungarian(cost)
print(row_ind, col_ind)  # should be [0,1,2] with min cost 1+4+9=14

But the optimal assignment might be off-diagonal? Actually, all costs are linearly dependent, so any permutation gives 1+4+9, 1+6+6, 2+2+9, etc. The minimum is 1+4+9=14 (diagonal). The algorithm should output row_ind = [0,1,2].

Let’s also test a random square matrix against brute force for small n to verify correctness.

import itertools
import time

def brute_force(cost):
    n = cost.shape[0]
    best = np.inf
    best_perm = None
    for perm in itertools.permutations(range(n)):
        total = sum(cost[i, perm[i]] for i in range(n))
        if total < best:
            best = total
            best_perm = perm
    return best_perm, best

np.random.seed(0)
for n in [3,4,5]:
    cost = np.random.rand(n,n)*10
    hung_row, hung_col = hungarian(cost)
    hung_val = cost[hung_row, hung_col].sum()
    brute_perm, brute_val = brute_force(cost)
    print(f"n={n}: Hungarian cost = {hung_val:.2f}, Brute force cost = {brute_val:.2f}, Match={np.allclose(hung_val, brute_val)}")

The algorithm should pass.

6. Complexity Analysis: Why O(N^3) and Not O(N^4)

6.1 Naive Hungarian Algorithm: O(N^4)

The original Hungarian method in Kuhn’s 1955 paper computed the minimum cover of zeros after each dual adjustment by scanning the entire matrix to count zeros and find the maximum matching. This scanning step could take O(N^2) per iteration, and there could be O(N^2) iterations (since each iteration may add at least one zero, but the dual adjustments might be small). In the worst case, total O(N^4).

Specifically, the naive algorithm:

• For each row, find an augmenting path using BFS/DFS in the zero graph, which is O(N^2) per row.
• At each dual update, re-scan the matrix to find uncovered zeros.

That leads to O(N _ N _ N * N) in the worst case.

6.2 The Munkres O(N^3) Improvement

Munkres realized that the augmenting path search could be performed using a Dijkstra-like algorithm that runs in O(N^2) per row, and the dual updates can be merged into that search without an O(N^2) re-scan. The key observations:

• The inner loop for a given row visits each column at most once because of the used array.
• The update of minv and dual variables takes O(N) per iteration of the while loop, and the while loop runs at most O(N) times (each time it visits a new column).
• So per row, we have O(N * N) = O(N^2) time. Since there are N rows, total O(N^3).

Moreover, the assignment update (augmentation) is O(N) per row.

Thus the algorithm is O(N^3).

6.3 Detailed Breakdown

For each row i:

1. Initialization: O(N).
2. While loop:
   • Each iteration marks a new column visited (so at most N iterations).
   • The inner loop over all unvisited columns: O(N) per iteration (but we maintain minv and scan all columns each time? Actually we scan all columns each time; that gives O(N^2) per row).
   • The dual update loop over all columns: O(N) per iteration.
   • So the while loop body is O(N) per iteration, giving O(N^2) per row.
3. Path augmentation: O(N) to trace back.

Overall O(N * N) = O(N^2) per row, summing to O(N^3).

6.4 Can We Do Better?

The assignment problem is a special case of minimum cost flow in a bipartite graph. The best known theoretical bound for dense assignment is O(N^3) (Munkres 1957). For sparse matrices, we can use a variant of the Successive Shortest Path algorithm with Dijkstra and potentials, achieving O(N^2 log N + NM) where M is number of nonzero costs. For dense matrices, O(N^3) is optimal in the comparison model.

7. Practical Considerations

7.1 Floating Point and Large N

Our implementation uses float. For integer costs, it is safe. For floating point, tiny errors may cause the algorithm to behave incorrectly (e.g., not finding a perfect matching because of near-zero values). A common fix is to add a small epsilon when comparing to zero. For instance, we can modify the condition to treat cur < 0 as zero after subtracting a tiny tolerance. But the algorithm relies on the nonnegativity of reduced costs; if you use floating point, you must ensure that all operations preserve nonnegativity (e.g., by using np.nextafter to avoid negative zeros). In practice, using double precision and checking for abs(cur) < 1e-12 in the comparison can work.

7.2 Handling Non-Square Matrices (Unbalanced Assignment)

Often we have more tasks than agents or vice versa. The Hungarian algorithm assumes square matrices. To handle unbalanced problems:

• If (m < n) (more tasks than agents), we add dummy rows (agents) with zero cost for all tasks (or very large cost if the dummy assignment is not allowed).
• If (m > n), we add dummy columns (tasks) with zero cost.

But this increases the matrix size to max(m,n). For extremely unbalanced cases (e.g., 10 agents vs 1000 tasks), the O(N^3) time can be prohibitive. In such cases, we can use a different algorithm: the Jonker-Volgenant (LAPJV) algorithm, which specializes in large rectangular assignments and runs in O(N^3) in the worst case but often faster in practice.

7.3 Maximization Problem

To maximize profit rather than minimize cost, we can convert profit matrix P to cost by (C*{ij} = M - P*{ij}) where M is the maximum profit, or simply negate all entries and run the minimization algorithm (since the Hungarian algorithm works with negative costs as long as they are not too negative to cause floating point overflow). Alternatively, you can modify the dual variables to start with row maxima. The standard trick: subtract each row from its maximum, then proceed as usual.

7.4 Parallelization and GPU

Because the Hungarian algorithm is inherently sequential (it processes rows one by one, updating dual variables), it is not easily parallelizable in its standard form. However, for very large matrices (N>10,000), one can use the Auction algorithm (Bertsekas) which is highly parallelizable and runs in O(N^3) worst case but can be sped up with GPUs. The Hungarian algorithm remains the gold standard for moderate sizes.

8. Real-World Applications

8.1 Disaster Relief Routing (the opening scenario)

Using the Hungarian algorithm, we can assign vehicles to incident locations in O(500^3) = 125 million operations, which on a modern CPU takes a few hundred milliseconds. That is fast enough for real-time dispatch.

8.2 Ride-Sharing and Taxi Dispatch

Uber and Lyft have to match drivers to riders in real time. The assignment is often one-to-one (one driver per rider) but may involve multiple riders per vehicle (pooling) which is a more general vehicle routing problem. For the basic dispatch, Hungarian can be used to minimize total wait time or maximize revenue.

8.3 Multi-Target Tracking (Computer Vision)

In tracking objects across video frames, we need to match detections in consecutive frames. The Hungarian algorithm assigns each tracked object to the most likely new detection using a cost based on predicted positions (from Kalman filters). This is a core component of the SORT (Simple Online and Realtime Tracking) algorithm.

8.4 Job Scheduling and Resource Allocation

In cloud computing, assigning virtual machine instances to physical servers to minimize power usage or maximize resource utilization can be modeled as an assignment problem. Similarly, in manufacturing, assigning jobs to machines.

8.5 DNA Sequence Alignment

In bioinformatics, the Hungarian algorithm is used in the “assignment problem” for aligning reads to a reference genome when the cost is edit distance, though faster heuristics are often used.

9. Variations and Further Reading

• Jonker-Volgenant algorithm (LAPJV): A variant of the Hungarian algorithm that uses more aggressive dual updates and works well for rectangular problems. It is the default in scipy.optimize.linear_sum_assignment for non-square matrices.
• Auction algorithm: Developed by Dimitri Bertsekas, it uses bids and price updates. It’s parallelizable and can be faster for very large N.
• Shortest Augmenting Path algorithm with potentials (e.g., using Dijkstra): This is essentially what we implemented, but can be extended to handle sparse matrices efficiently.
• Assignment problem with side constraints: When assignments cannot be arbitrary (e.g., due to skill requirements, time windows), the problem becomes a more complex integer program. Lagrangian relaxation often yields subproblems that are assignment problems.

10. Conclusion

We have journeyed from the conceptual depths of the assignment problem—where brute force is laughably impossible—to the elegant, polynomial-time solution of the Hungarian algorithm. We saw how Kuhn and Munkres transformed combinatorial optimization with a primal-dual method that runs in O(N^3) time. We implemented the algorithm line by line, dissected its complexity, and explored where it fits in modern high-performance computing.

The Hungarian algorithm is a testament to the power of mathematical structure: what seems like an intractable combinatorial explosion yields to a clever blend of linear programming theory, graph theory, and careful bookkeeping. Next time you request a ride, a satellite is tasked to cover a region, or a robot maps its environment, remember that somewhere, an O(N^3) implementation of the Hungarian algorithm might be silently solving an assignment problem to optimize your world.

If you want to dive deeper, I encourage you to read the original papers by Kuhn (1955) and Munkres (1957), or check out the source code of scipy.optimize and the munkres Python package. And next time you need to assign a hundred tasks to a hundred workers, you won’t need a miracle—you’ll need the Hungarian algorithm.

Happy optimizing!

This blog post was written to provide a thorough, implementation-focused explanation of the Hungarian algorithm. All code examples are available in a Jupyter notebook linked below. Comments, corrections, and further insights are welcome.