Hungarian Algorithm, in depth: the LP dual, reduced costs, two hand traces, O(n^3) shortest augmenting paths and production use

The linear program and its dual

Write x(i,j) = 1 if row i is assigned to column j. The primal problem minimizes the sum of c(i,j) x(i,j) subject to every row summing to 1, every column summing to at most 1 (exactly 1 when the matrix is square), and x non-negative. The constraint matrix is the incidence matrix of a bipartite graph, which is totally unimodular, so the linear program has an integral optimal vertex and the relaxation loses nothing. That fact is what makes a polynomial exact algorithm possible.

The dual assigns a potential u(i) to each row and v(j) to each column and maximizes the sum of all potentials subject to u(i) + v(j) <= c(i,j) for every pair. Define the reduced cost c'(i,j) = c(i,j) - u(i) - v(j); the dual constraint says every reduced cost is non-negative. For any assignment, the cost equals the sum of potentials plus the sum of reduced costs on the chosen pairs, so the sum of potentials is a lower bound on every assignment. Complementary slackness then gives the optimality test: if a perfect assignment uses only pairs with reduced cost zero, called tight pairs, its cost equals the lower bound and nothing can beat it. The Hungarian algorithm is a procedure that keeps the potentials feasible and grows a matching made only of tight pairs until it is perfect.

Worked example: reductions as potentials

The textbook matrix version makes the potentials visible. Take the 4 by 4 matrix with rows (9, 2, 7, 8), (6, 4, 3, 7), (5, 8, 1, 8) and (7, 6, 9, 4). Subtracting each row's minimum, 2, 3, 1 and 4, is the same as setting u = (2, 3, 1, 4); the reduced rows become (7, 0, 5, 6), (3, 1, 0, 4), (4, 7, 0, 7) and (3, 2, 5, 0). Column 0 has no zero, so subtract its minimum, 3, which sets v(0) = 3 and leaves the other column potentials at 0. The reduced matrix is now (4, 0, 5, 6), (0, 1, 0, 4), (1, 7, 0, 7), (0, 2, 5, 0). The potentials sum to 13, so no assignment costs less than 13.

Now look for a perfect matching on the zeros. Row 0 has only column 1, row 2 has only column 2, so row 1 takes column 0 and row 3 takes column 3. The original costs are 2 + 6 + 1 + 4 = 13, equal to the bound, so this assignment is optimal, and the potentials are its proof. The diagram shows the tight edges and the matching chosen among them.

When zeros are not enough: the cover step

Reductions alone are not always enough. Take rows (1, 2, 3), (2, 4, 6) and (3, 6, 9). Row reduction gives (0, 1, 2), (0, 2, 4), (0, 3, 6); column reduction subtracts 0, 1 and 2, leaving (0, 0, 0), (0, 1, 2), (0, 2, 4) with a bound of 9. Rows 1 and 2 both have their only zero in column 0, so the largest tight matching has size 2. By Konig's theorem a maximum matching of size 2 has a vertex cover of size 2: row 0 and column 0 cover every zero.

The classical step takes the smallest uncovered value, here 1, subtracts it from every uncovered entry and adds it to entries covered twice. In dual terms that raises the potential of each uncovered row by 1 and lowers the potential of each covered column by 1, which keeps every reduced cost non-negative, creates at least one new zero, and raises the bound by delta times (n minus the cover size), here 1 times (3 - 2) = 1. The matrix becomes (1, 0, 0), (0, 0, 1), (0, 1, 3), the matching row 2 to column 0, row 1 to column 1, row 0 to column 2 is now tight, and its cost 3 + 4 + 3 = 10 equals the new bound. A brute-force check of all six permutations confirms 10 is the minimum.

Finding a minimum cover by hand at each step is where the matrix method gets slow and error-prone. The practical algorithm reorganizes the same dual adjustment around shortest paths.

The O(n^3) shortest augmenting path algorithm

The O(n^3) version adds rows one at a time. When row i arrives, it searches for an augmenting path from i to a free column, using reduced costs as edge lengths, in the style of Dijkstra's algorithm. It keeps minv[j], the smallest reduced cost seen to reach column j from the tree so far. Each iteration picks the unvisited column with the smallest minv, shifts all potentials by that amount so the chosen edge becomes tight while tree edges stay tight, and either reaches a free column, in which case it flips the path, or follows the column's current partner row and keeps growing. The code below follows the well-known compact formulation with a virtual column 0; it was checked against brute force on 600 random matrices, including negative costs and rectangular shapes.

INF = float("inf")

def hungarian(cost):
    """Min-cost assignment for an n x m matrix with n <= m.
    Returns (total, assign) where assign[i] is the column given to row i."""
    n, m = len(cost), len(cost[0])
    assert n <= m, "transpose the matrix first"
    u = [0.0] * (n + 1)          # row potentials (index 0 unused)
    v = [0.0] * (m + 1)          # column potentials; v[0] is the virtual column
    match = [0] * (m + 1)        # match[j] = row matched to column j (0 = free)
    way = [0] * (m + 1)          # predecessor column on the shortest path
    for i in range(1, n + 1):
        match[0] = i             # root the search tree at row i
        j0 = 0
        minv = [INF] * (m + 1)   # tentative reduced distance to each column
        used = [False] * (m + 1)
        while True:
            used[j0] = True
            i0, delta, j1 = match[j0], INF, 0
            for j in range(1, m + 1):
                if not used[j]:
                    cur = cost[i0 - 1][j - 1] - u[i0] - v[j]   # reduced cost
                    if cur < minv[j]:
                        minv[j], way[j] = cur, j0
                    if minv[j] < delta:
                        delta, j1 = minv[j], j
            for j in range(m + 1):                # dual adjustment
                if used[j]:
                    u[match[j]] += delta
                    v[j] -= delta
                else:
                    minv[j] -= delta
            j0 = j1
            if match[j0] == 0:                    # free column reached
                break
        while j0:                                 # augment along the path
            j1 = way[j0]
            match[j0] = match[j1]
            j0 = j1
    assign = [0] * n
    for j in range(1, m + 1):
        if match[j]:
            assign[match[j] - 1] = j - 1
    return sum(cost[i][assign[i]] for i in range(n)), assign

C = [[9, 2, 7, 8], [6, 4, 3, 7], [5, 8, 1, 8], [7, 6, 9, 4]]
print(hungarian(C))   # (13, [1, 0, 2, 3])

Complexity and certificates

Each phase adds one row to the matching. Inside a phase, every iteration marks one new column as used, and every column it can reach before a free one is already matched, so phase i runs at most i iterations, each scanning m columns. The whole run is O(n^2 m), which is O(n^3) for a square matrix. Memory is O(n m) for the matrix plus O(m) of working arrays. The naive matrix method, which recomputes a maximum matching and a cover after each adjustment, is O(n^4) and worse in practice. For a 1,000 by 1,000 matrix the cubic bound is about a billion elementary steps: fine in C, slow in pure Python.

The final potentials are a certificate you can check independently in O(n m): verify that every reduced cost is non-negative, within a floating point tolerance, and that every assigned pair has reduced cost zero. Doing this in tests catches implementation bugs that a handful of hand-picked examples will not.

Variants: rectangular, maximization, forbidden pairs

Rectangular matrices. With more columns than rows the code above works directly, leaving some columns unused. With more rows than columns, transpose, solve, and invert the result. Padding with dummy zero-cost columns also works but makes the matrix larger.

Maximization. To maximize total score, negate the matrix, or subtract every entry from the maximum entry. Both give the same assignment; negation keeps the total easier to read back.

Forbidden pairs. SciPy accepts infinite entries and raises ValueError when no feasible assignment exists. In a hand-rolled solver like the one above, use a large finite cost M instead, because infinite entries poison the potential arithmetic. Choose M larger than n times the largest real cost, then check the answer: if any assigned pair costs M, no feasible perfect assignment exists and the caller must decide what partial result means.

Thresholded matching. Trackers often want 'assign only if cost is below tau'. Solve with forbidden pairs above tau, then drop any assignment at cost M.

In practice: libraries and alternatives

In Python, use scipy.optimize.linear_sum_assignment. Its documentation describes the implementation as a modified Jonker-Volgenant algorithm with no initialization, after Crouse's 2016 paper on rectangular assignment. Jonker-Volgenant is a shortest-augmenting-path method in the same family as the code above, with tuned search and bookkeeping. It accepts rectangular matrices and a maximize flag and returns row and column index arrays.

import numpy as np
from scipy.optimize import linear_sum_assignment

C = np.array([[9, 2, 7, 8], [6, 4, 3, 7], [5, 8, 1, 8], [7, 6, 9, 4]])
rows, cols = linear_sum_assignment(C)
print(cols, C[rows, cols].sum())          # [1 0 2 3] 13

# Tracking-style gating: forbid pairs above a distance threshold.
TAU, M = 50.0, 1e6
D = np.random.rand(30, 40) * 100           # detections x tracks
G = np.where(D <= TAU, D, M)
r, c = linear_sum_assignment(G)
keep = G[r, c] < M                         # drop forced, forbidden pairs
matches = list(zip(r[keep], c[keep]))

Dense Hungarian is the right tool up to a few thousand rows with a full cost matrix. Beyond that, or when most pairs are forbidden, model the problem as min-cost flow on the sparse graph, or use an auction algorithm, which parallelizes well. When the costs are all equal and only the size of the matching matters, an unweighted algorithm is far faster. Greedy assignment is tempting but carries no optimality guarantee.

Failure modes

• Infinity in a hand-rolled solver. inf - inf produces NaN and the potentials silently go wrong. Use a large finite M and check for it afterwards.
• M too small. If M is below the cost of a genuinely cheaper reshuffle, the solver will use a forbidden pair to save money elsewhere. Size M from the data.
• Wrong orientation. Passing a tall matrix to code that assumes n <= m fails or misreports; transpose and map indices back explicitly.
• Maximizing by accident. Similarity matrices fed to a minimizer produce the worst assignment, which still looks plausible. Test on a case with an obvious answer.
• Unstable ties. Equal costs admit several optimal assignments; frame-to-frame trackers then flicker. Add a tiny deterministic tie-breaker, such as a penalty for changing the previous assignment.
• Cubic surprise. A batch that grows from 500 to 5,000 rows becomes a thousand times slower. Put a size guard and a fallback in the code path.

Trade-offs

What to do next

1. Reproduce the two hand traces above, then run the code and confirm the totals 13 and 10.
2. Add a brute-force cross-check over random small matrices to your own implementation's tests.
3. Write a certificate checker that verifies non-negative reduced costs and tight assigned pairs.
4. Read bipartite matching in depth for the unweighted problem and the flow reduction behind it.
5. Work through Kuhn's algorithm to see the augmenting path idea without weights.
6. Study Konig's theorem and vertex cover to understand the cover step in the matrix method.
7. Model a sparse version of your problem with min-cost flow and compare its running time with the dense solver.
8. In production, call linear_sum_assignment, gate forbidden pairs with a finite M, and add a size guard.