Algorithms

The linear assignment problem

Given an n × n cost matrix C with entries c_ij ∈ ℝ ∪ {+∞}, the linear assignment problem (LAP) asks for a permutation σ of {1, ..., n} that minimizes

sum_{i=1..n} c_{i, σ(i)}

over all permutations. Rectangular variants (n × m with n ≠ m) match every row to a distinct column when n ≤ m, or every column to a distinct row in the transposed case. torchmatch reduces rectangular instances to the square case by padding with a sentinel that any optimal assignment avoids (see Implementation notes).

The problem has two equivalent combinatorial faces. As minimum-weight perfect matching, think of rows as one set of nodes and columns as another; every (row, column) pair is a possible edge with weight c_ij. The LAP finds a set of n edges — one per row, one per column — with minimum total weight. As a linear program, relax the integrality constraint on the indicator variables x_ij ∈ {0, 1}:

min   sum_{ij} c_ij * x_ij
s.t.  sum_j x_ij = 1   for all i  (each row matched once)
      sum_i x_ij = 1   for all j  (each column matched once)
      x_ij >= 0

Birkhoff's theorem guarantees that, despite allowing x_ij to be any non-negative real number, every corner point of the feasible region is already a 0/1 assignment matrix. Solvers that find an optimal corner — as standard LP solvers do — therefore return an integer assignment directly, with no rounding needed.

The dual LP introduces variables u_i (one per row) and v_j (one per column):

max   sum_i u_i + sum_j v_j
s.t.  u_i + v_j <= c_ij   for all i, j

Define the reduced cost c̄_ij = c_ij − u_i − v_j. A feasible dual satisfies c̄_ij ≥ 0 everywhere. The equality subgraph G_= (the subgraph containing only edges whose reduced cost is zero — the candidates for the optimal matching) keeps exactly the edges with c̄_ij = 0. This is the key insight driving every Hungarian variant: an assignment is optimal if and only if every edge it uses has zero reduced cost, i.e., the current prices (u, v) exactly account for each matched cost. Equivalently, an assignment is optimal if and only if it is a perfect matching in G_= for some feasible dual (u, v).

The Hungarian method: history

The combinatorial roots are theorems on bipartite graphs by König [Konig1931] and their weighted extension by Egerváry [Egervary1931], both published in Hungarian in 1931. König proved that in a bipartite graph the maximum matching size equals the minimum vertex cover size; Egerváry extended the duality to weighted edges. Neither paper described an algorithm, but the structural results were what a primal-dual procedure would later exploit.

Kuhn formalized the algorithmic content in 1955 [Kuhn1955], coining the name "Hungarian method" in tribute to the two earlier results. His version solved the LAP in O(n^4) time. Munkres [Munkres1957] simplified the bookkeeping into a six-step procedure that tracks which matrix entries are candidates for the matching ("starred" zeros) and which are being explored in the current search step ("primed" zeros), using row and column flags ("covers") to mark already-matched rows and columns; the restatement made the algorithm easier to implement and remains the canonical textbook presentation. Independently, Tomizawa [Tomizawa1971] and Edmonds and Karp [EdmondsKarp1972] showed that choosing which alternating path to augment along — a path that alternates between matched and unmatched edges, flipping each to grow the matching by one — can change the total work; choosing shortest such paths lowers the bound to O(n^3), the modern asymptotic.

Subsequent work refined the inner loop without changing the asymptotic. Lawler [Lawler1976] rephrased the augmentation as a breadth-first search that discovers all vertex-disjoint augmenting paths simultaneously per outer iteration, trading more per-step work for fewer outer iterations and greater parallelism. Jonker and Volgenant [JonkerVolgenant1987] gave a Dijkstra-based successive-shortest-path variant whose practical constant factor is small enough that it became the default CPU solver for dense LAPs; their column-reduction warm start and reduction-transfer heuristic cut a noticeable fraction of the per-iteration work. Kuhn's own retrospective [Kuhn2010Variants] discusses the chain of attribution and is a useful entry point for the historiography.

Two algorithmic families address the LAP without augmenting paths. Bertsekas's auction algorithm [Bertsekas1979]; [Bertsekas1988] casts the problem as distributed bidding with an ε-scaling parameter that controls convergence precision; it parallelizes uniformly but the number of iterations depends on ε and can grow large for high-precision solutions. Push-relabel-style cost-scaling [GoldbergKennedy1995] reformulates the LAP as a network flow problem and applies flow-based optimization techniques. Both have full coverage in the secondary literature [BurkardDellAmicoMartello2012] and are absent from torchmatch.

The primal-dual template

The dual variables (u, v) act as prices or potentials assigned to rows and columns: they encode how much each row and column "should cost" in an optimal assignment and guide the search toward edges that can be in the optimum. The primal-dual approach works by maintaining feasible prices at all times and improving the matching only along edges where the price exactly accounts for the cost — guaranteeing that each improvement step moves toward optimality without backtracking.

Every Hungarian variant in torchmatch instantiates the same primal-dual template:

1. Initialize duals (u, v) so that c̄_ij ≥ 0 everywhere. The standard choice is u_i = min_j c_ij (row reduction) followed by v_j = min_i (c_ij − u_i) (column reduction). This places at least one zero in every row and column.
2. Match in the equality subgraph. Mark (star) a maximal set of zero-reduced-cost edges with no two sharing a row or column — this gives the initial partial matching M in G_=.
3. Augment. While M is not perfect, search G_= for an alternating path from an unmatched row to an unmatched column; if found, augment M along it. Each augmentation grows |M| by one.
4. Dual update. If no augmenting path exists, find the minimum reduced cost δ over edges that connect the current search frontier to columns not yet reachable, then subtract δ from the row potentials and add it to the column potentials ("shift duals") so that at least one previously unreachable edge now has zero reduced cost and can be explored. Return to step 3.

Total work is bounded by n outer augmentations, each costing O(n^2) for the search and the dual update, giving O(n^3). The three torchmatch implementations differ in how step 3 (the path search) is organized and in what data structures step 4 maintains; the correctness argument is the same.

Munkres' classical (single-path)

The implementation in cuda/munkres.cu, exposed as munkres, follows the Munkres [Munkres1957] six-step state machine. The state consists of three pieces of bookkeeping:

• A 0/1 cover bit per row and per column.
• A starred mark on at most one zero per row and per column, representing the current matching.
• A primed mark used to construct an augmenting path.

Each outer iteration finds at most one new starred zero, so the algorithm terminates in at most n outer iterations. Inner work per iteration is dominated either by the augmenting-path search (linear in the number of currently uncovered zeros) or by the dual update (O(n^2) over the uncovered submatrix). [Tomizawa1971]; [EdmondsKarp1972] tightened the original O(n^4) bound to O(n^3) through a careful choice of augmenting paths.

On GPU, per-iteration work decomposes as follows. Column-min and matrix-min reductions use CUB's BlockReduce and DeviceReduce primitives. The compressed-zeros structure — a compact list of near-zero-cost edges — avoids rescanning the full matrix at every iteration. Augmentation runs serially within an iteration; the surrounding reductions parallelize cleanly. The implementation matches Munkres' state machine in shape, with kernel fusion for the dual update plus the next iteration's zero recompression; see cuda/munkres.cu for the details.

When it wins. Munkres' classical is fastest on sparse or tied LAP instances: integer-valued costs from a small support, gated edge sets after spatial or appearance filtering, or anything that produces many equality-subgraph edges per iteration. The single-path inner loop has the lowest per-iteration overhead of the three CUDA solvers, which pays off when the augmenting paths are short.

Lawler's tree augmentation

The implementation in cuda/lawler.cu, exposed as lawler, follows Lawler's [Lawler1976] reformulation. Instead of finding one augmenting path per outer iteration, a breadth-first search builds all vertex-disjoint augmenting paths reachable from unmatched rows simultaneously, then augments along each. The asymptotic complexity remains O(n^3); the work per outer iteration is larger and the number of outer iterations is correspondingly smaller.

The BFS structure suits GPU execution. Each BFS level is expanded in parallel across GPU threads; a standard parallel-prefix operation (Thrust exclusive_scan) compactly gathers the newly discovered columns into a list for the next level, in O(n) work per level. The implementation issues a single cooperative-launch kernel that fuses the entire expansion loop, avoiding host syncs between BFS levels; an offline fallback path uses a host-driven multi-kernel sequence for cases where the cooperative launch cannot fit (very large n).

When it wins. Lawler's variant is fastest on dense LAP instances at large n (≥ 512 in the current benchmark sweep). The single-path classical procedure runs largely serially within an iteration; tree augmentation exposes enough parallelism to saturate the SMs.

Jonker-Volgenant (successive shortest path)

The CPU implementations in cpu/jonker_*.{h,cpp}, exposed as jonker_scalar, jonker_dense, jonker_compact, follow the Jonker-Volgenant procedure [JonkerVolgenant1987]. The mechanical structure differs from the primed/starred-zeros family but the algebraic substrate (dual variables, equality subgraph, alternating paths) is identical.

For each unmatched row i, run a Dijkstra-like search over the reduced cost graph: maintain tentative shortest distances d_j from i to each column, expand the column with smallest d_j whose row predecessor lies on the current alternating path, and stop when an unmatched column is reached. The search costs O(n^2) per row in the dense case; the total is O(n^3). The dual update folds into the search: at the end, set u_i ← u_i + d_J − d_{J'} for each row on the augmenting path (here J is the reached unmatched column and J' is the column it was matched to before the update). This preserves dual feasibility and produces the new equality-subgraph edges required for the next iteration.

The JV speedup over Munkres' classical on dense problems comes from three places. First, the column-reduction warm start finds a good initial matching cheaply, often matching O(n) rows before the SSP loop begins. Second, the per-row Dijkstra terminates as soon as an unmatched column is reached; on smooth cost distributions this settles after a small fraction of n expansions. Third, the reduction-transfer step in JV's column-reduction phase reuses dual updates from one iteration to seed the next, cutting per-iteration arithmetic.

Crouse [Crouse2016] extended JV to non-square cost matrices without the pad-and-solve overhead. The rectangular formulation tracks dual variables only on the min(n, m) side and reformulates the inner loop for the asymmetric structure. torchmatch's jonker_dense uses this rectangular variant; jonker_compact is square-only internally and the wrapper pads when needed.

Within the JV family, torchmatch ships three CPU variants:

The batched variants (jonker_dense_batch, jonker_compact_batch, and the _unpacked flavours) run at::parallel_for across independent problems on CPU. jonker_dense_batch also has a CUDA backend, a shared-memory tiled kernel that solves one problem per CUDA block, constrained to (B, K, K) square with K ≤ MAX_TILE. The tiled kernel keeps all working data in fast on-chip shared memory for the duration of a single problem, avoiding repeated reads from the slower global GPU memory that a straightforward per-problem kernel launch would require.

When it wins. Almost everywhere except integer-tied costs. On dense LAPs at n ≤ 1024, the AVX2 inner loop beats both CUDA Hungarian implementations by 10× to 100× on a current consumer GPU, since the algorithm needs neither host syncs nor cross-block communication.

Comparing the three

All three implementations solve the same problem, return optimal assignments, and share the same O(n^3) worst-case complexity. They differ in the constant factor and in the parallelism profile.

Practical picking guidance lives in Choosing the right op; the empirical numbers behind those rules live in Benchmarks.

Implementation notes

A few invariants are shared across all torchmatch solvers.

+∞ rewriting. Forbidden edges are encoded as +∞ on input. The wrapper rewrites +∞ cells to a per-call sentinel (max_finite + 1) · (K + 1) where K = max(n, m). Any match that routes through a sentinel is provably worse than any all-finite assignment, so optimal solvers avoid such cells; the output repacker maps any sentinel-routed match to −1 ("row unmatched"). NaN inputs are rejected at the boundary, since NaN carries no defensible forbidden-edge meaning.

Padding to square. Rectangular (n, m) inputs are padded with the sentinel to (K, K). The padding is invisible at the API boundary: outputs have length n, and matches into padded columns become −1.

Workspace caching. Each backend owns a process-level WorkspaceCache<(K, device) → Layout> that amortizes the per-call cudaMalloc cost. The first call at a given (K, device) constructs the workspace; subsequent calls reuse it. On the CUDA side, host-readable iteration flags coalesce into a single managed ControlBlock per workspace slot so per-iteration host reads avoid a cudaMemcpy round trip.

CUDA graph compatibility. The three CUDA Hungarian ops (munkres, hybrid, lawler) carry the cudagraph_unsafe tag: their inner loops do host-side cudaStreamSynchronize to read managed-memory iteration flags. Under torch.compile(mode="reduce-overhead") they trigger a graph break. The CUDA backend of jonker_dense_batch has no host syncs and is fully capturable.