Min-Cost Bipartite Matching, in depth: epsilon-complementary slackness, the auction algorithm, cost scaling and price wars

The assignment problem as prices

Take an n by n integer cost matrix c. A perfect matching picks one column for each row with no column repeated. Minimising total cost is a linear program whose relaxation already has integral optimal vertices, so duality applies cleanly. It is easiest to state the auction as maximisation, so set the benefit a[i][j] = -c[i][j] and give every column a price p[j]. The net value of column j to row i is a[i][j] - p[j].

Exact optimality says every row holds a column of maximum net value; such prices are an optimal dual. The Hungarian method keeps that condition exact while growing the matching. Cost scaling lets the condition be slightly loose instead.

Epsilon-complementary slackness

Epsilon-complementary slackness (eps-CS). A row i assigned to column j satisfies eps-CS if a[i][j] - p[j] >= max_k (a[i][k] - p[k]) - eps. Every row is within eps of its best deal at current prices.

The bound. If a perfect matching satisfies eps-CS for every row, its total benefit is within n times eps of the optimum. Sum the inequality over all rows: the matched benefit plus the sum of all prices is at least the dual objective (the sum of each row's best net value plus all prices) minus n eps, and the dual objective is an upper bound on any matching's benefit by weak duality.

Exactness by integer scaling. Multiply every cost by n + 1 before starting. Now any two different matchings differ in total by a multiple of n + 1, so a matching within n eps of the optimum with eps = 1 is within n, which is less than n + 1, and must be optimal. That is the whole trick that lets the algorithm run on integers with eps never smaller than 1. Floating-point eps with a tolerance is where most homemade auctions go wrong.

The auction: bids and prices

A phase starts with every row free. A free row i scans its columns, finds the best net value v1 at column j and the second best v2, raises p[j] by v1 - v2 + eps, and takes column j. If someone held j, that row is evicted and rejoins the queue. Raising by exactly that amount makes j worth v2 - eps to row i, which keeps eps-CS for the new owner, and every price rises by at least eps per bid, which is what forces termination.

This is the code that produced every number on this page:

from collections import deque

def auction_phase(benefit, price, eps, stats):
    """One epsilon phase of the Gauss-Seidel auction. benefit: n x n ints."""
    n = len(benefit)
    owner = [-1] * n            # owner[j] = row holding column j
    assigned = [-1] * n         # assigned[i] = column held by row i
    free = deque(range(n))
    while free:
        i = free.popleft()
        row = benefit[i]
        best_j, v1, v2 = -1, None, None
        for j in range(n):
            v = row[j] - price[j]
            if v1 is None or v > v1:
                v2, v1, best_j = v1, v, j
            elif v2 is None or v > v2:
                v2 = v
        if v2 is None:
            v2 = v1 - eps        # n == 1
        price[best_j] += v1 - v2 + eps
        stats["bids"] += 1
        prev = owner[best_j]
        owner[best_j] = i
        assigned[i] = best_j
        if prev != -1:
            assigned[prev] = -1
            free.append(prev)
    return assigned

def min_cost_assignment(cost, theta=5):
    """Minimise sum cost[i][assign[i]] over perfect matchings; integer costs."""
    n = len(cost)
    scale = n + 1
    benefit = [[-c * scale for c in row] for row in cost]   # maximise -cost
    C = max(abs(b) for row in benefit for b in row) or 1
    price = [0] * n
    eps = max(1, C // theta)
    stats = {"bids": 0, "phases": 0}
    while True:
        assigned = auction_phase(benefit, price, eps, stats)
        stats["phases"] += 1
        if eps == 1:
            break
        eps = max(1, eps // theta)
    total = sum(cost[i][assigned[i]] for i in range(n))
    return assigned, total, stats

Each new phase drops the matching but keeps the prices, and the loop always ends with a phase at eps = 1.

Worked example: a 3 by 3 trace

Take the cost matrix with rows [9, 2, 7], [6, 4, 3] and [5, 8, 1]. By inspection the optimum is 2 + 6 + 1 = 9: row 0 to column 1, row 1 to column 0, row 2 to column 2. Scaling by n + 1 = 4 gives benefits from -36 to -4, so C = 36 and the first eps is 36 // 5 = 7.

The coarse phase ends at cost 2 + 3 + 5 = 10, one unit above optimal and inside its n eps guarantee of 21 scaled units (5.25 cost units). The eps = 1 phase reuses the prices, needs four bids including one eviction, and lands on the optimum of 9, matching scipy.

Why scale epsilon: price wars, measured

An unscaled auction at eps = 1 is still correct, so why scale? Price wars. When several rows want the same few columns with nearly tied values, each bid raises a price by about eps, the evicted row bids on a near-identical column, and prices creep up one eps at a time. Bids grow with C divided by eps, which is pseudo-polynomial.

At large eps the same prices move in big jumps and settle roughly where they belong; each later phase only fixes errors the size of the previous eps. The theory gives O(n m log(nC)) for the scaled auction on m edges. The measurement is more nuanced. Bids counted with the code above, totals verified against scipy:

On random dense matrices the unscaled auction is often cheaper: the costs are spread out, so wars are short, and scaling pays for its extra phases. At n = 400 with costs 0..100 there are many near-ties and the unscaled run needed five times as many bids. With identical rows, where every row ranks the columns the same way, it is a rout: 2.86 million bids against 1,670 at n = 80, roughly 85 seconds against 0.06 in Python. Treat scaling as insurance against the worst case, not as a speedup on every input. For theta, 4 gave the fewest bids on the identical-row tests (1,565 at n = 80, against 1,989 for 2 and 2,772 for 16), so moderate factors win and the curve is flat enough that tuning it is rarely worth much.

The cost-scaling family

The auction has relatives you will meet in papers and solver docs:

• Push-relabel with cost scaling. Goldberg and Tarjan generalised the idea to min-cost flow. A bid is a push of one unit of flow plus a relabel, which is a price rise. The live push-relabel article covers the max-flow version that this builds on.
• CSA (Goldberg and Kennedy, 1995). A cost-scaling assignment code with heuristics such as price updates and arc fixing. It is the reference point for large sparse assignment in benchmarks.
• Gabow and Tarjan (1989). Scaling plus Hopcroft-Karp style phases of shortest augmenting paths, giving O(sqrt(n) m log(nC)). It is mainly a theoretical bound; practical codes rarely implement it directly.
• Jacobi (parallel) bidding. All free rows compute bids at once and each column accepts the highest. This is what GPU implementations exploit, because the bid computation is a row-wise max and second-max.

Sparse, infeasible and rectangular inputs

Real instances are usually sparse: a detection can only match tracks within a gate. A bid then scans only the row's own edges. Two traps come with sparsity.

First, infeasibility: if no perfect matching exists, some rows will bid forever while prices climb without bound. Check feasibility first with Hopcroft-Karp from the live bipartite matching article, or add a dummy column per row at a large finite cost so unmatched rows have somewhere to go and you can read them back as unassigned. Second, single-edge rows: a row with one edge has no second-best value. Treat the missing v2 as minus infinity capped at a finite floor, or the bid increment overflows.

Rectangular problems need care too: pad with zero-cost dummy rows, or use the forward-reverse auction, in which columns also bid for rows.

Failure modes

• Floating-point eps. With real costs and a tiny eps, rounding can make the increment zero or negative and the loop never ends. Scale to integers (for example multiply by 1,000 and round) before running.
• Forgetting the n + 1 scale. Without it, eps = 1 only guarantees a result within n of optimal. The trace above would have stopped at cost 10.
• Overflow. Costs times n + 1, plus prices that rise past the cost range, can exceed 32 bits on large instances. Use 64-bit integers and check that C times (n + 1) times a small factor fits.
• No termination on infeasible input. Put a bid budget on each phase and raise an error instead of spinning.
• Trusting the answer. Confirm eps-CS at eps = 1 and compare samples against a reference solver.

Trade-offs and when to use which

The auction's practical edge is the warm start. In multi-object tracking, frame t + 1 has nearly the same costs as frame t, so reusing the previous prices often makes the next solve a handful of bids. Most Hungarian implementations, scipy's included, accept no warm start. If your problem has extra constraints, model it as flow instead with the live min-cost flow article.

What to do next

1. Default to scipy.optimize.linear_sum_assignment for dense matrices; it is exact and its docstring documents the algorithm.
2. Switch to an auction or CSA code when the graph is sparse and large, when you solve a stream of similar instances, or when you need GPU parallelism.
3. Convert costs to integers, multiply by n + 1, and finish with a phase at eps = 1.
4. Start with theta between 4 and 8 and measure bids, not seconds, on your own data.
5. Check feasibility first or add dummy edges, and set a bid budget per phase.
6. Test against a reference solver on random and adversarial inputs, including identical rows, which is where unscaled auctions collapse.
7. Keep prices between solves when consecutive instances are similar.