Assignment problem (Hungarian / Jonker–Volgenant)

Deep dive

How it works

Parts:

• cost matrix C — n rows (workers) × m ≥ n columns (jobs), costs i64
• row potentials u, column potentials v — dual variables; reduced cost c_ij − u_i − v_j ≥ 0 always
• matching p — p[j] = row holding column j; the primal solution
• search arrays — minv[j] (cheapest reduced cost into column j), way[j] (predecessor column), used[j]

Operations, for each row i in turn:

1. Search from the new row. Each pass scans the cost row of the latest visited row and updates minv/way for unvisited columns.
2. Take the unvisited column with the smallest minv; call it δ. Add δ to u of every visited row and subtract it from v of every visited column; subtract δ from minv of the rest. Every visited edge stays ≥ 0 and the chosen edge becomes tight.
3. If that column is free, stop. If it is owned by a row, visit that row and go to step 1.
4. Walk way back from the free column and flip matched and unmatched edges along the alternating path. The matching grows by one and stays tight.

Steps 1-3 are Dijkstra on reduced costs (nonnegative, as Dijkstra needs), O(m·k) for k visited columns. Final cost = Σ over matched c_ij = Σu + Σv.

 reduced costs c_ij − u_i − v_j     [ ] = matched      inserting r1 when r2–c2 and r3–c3 are matched:
         c1    c2    c3                                 r1 ─tight→ c2 ═owned by═ r2 ─tight→ c1 (free)
 u=3 r1   1   [0]    0                                 flip: r1–c2 and r2–c1 become matched
 u=2 r2  [0]    0    3                                 r3 keeps c3
 u=2 r3   1    2    [0]
 v:       0   -2     0          Σu + Σv = 7 − 2 = 5 = cost of C = [[4,1,3],[2,0,5],[3,2,2]]

Code

const INF: i64 = i64::MAX / 4;

pub struct Assignment {
    pub col_of_row: Vec<usize>, // column chosen for each row
    pub total: i64,
    pub u: Vec<i64>, // row potentials, 1-based (index 0 unused)
    pub v: Vec<i64>, // column potentials, 1-based (index 0 is scratch)
}

/// Min-cost assignment of n rows to m >= n columns; `c` is row-major n*m.
/// Forbidden pairs: a large finite cost (e.g. 1 << 40), never i64::MAX.
pub fn hungarian(n: usize, m: usize, c: &[i64]) -> Assignment {
    assert!(n <= m && c.len() == n * m);
    let (mut u, mut v) = (vec![0i64; n + 1], vec![0i64; m + 1]);
    let mut p = vec![0usize; m + 1]; // p[j] = row matched to column j, 0 = free; p[0] = row being inserted
    let mut way = vec![0usize; m + 1]; // way[j] = previous column on the shortest alternating path to j
    for i in 1..=n {
        p[0] = i;
        let mut j0 = 0; // current column; column 0 is the virtual start holding row i
        let mut minv = vec![INF; m + 1]; // best reduced cost from the visited rows to each column
        let mut used = vec![false; m + 1];
        loop {
            used[j0] = true;
            let i0 = p[j0];
            let (mut delta, mut j1) = (INF, 0);
            for j in 1..=m {
                if used[j] { continue; }
                let cur = c[(i0 - 1) * m + (j - 1)] - u[i0] - v[j]; // reduced cost, >= 0
                if cur < minv[j] { minv[j] = cur; way[j] = j0; }
                if minv[j] < delta { delta = minv[j]; j1 = j; }
            }
            // Dijkstra step: shift potentials so the nearest column becomes tight (reduced cost 0)
            for j in 0..=m {
                if used[j] { u[p[j]] += delta; v[j] -= delta; } else { minv[j] -= delta; }
            }
            j0 = j1;
            if p[j0] == 0 { break; } // free column reached: augment
        }
        while j0 != 0 { let j1 = way[j0]; p[j0] = p[j1]; j0 = j1; } // flip the alternating path
    }
    let mut col_of_row = vec![0; n];
    for j in 1..=m { if p[j] != 0 { col_of_row[p[j] - 1] = j - 1; } }
    let total = (0..n).map(|i| c[i * m + col_of_row[i]]).sum();
    Assignment { col_of_row, total, u, v }
}

/// Maximum-weight assignment by negating the weights.
pub fn hungarian_max(n: usize, m: usize, w: &[i64]) -> Assignment {
    let neg: Vec<i64> = w.iter().map(|x| -x).collect();
    let mut a = hungarian(n, m, &neg);
    a.total = -a.total;
    a
}

/// Any shape: with more rows than columns, solve the transpose; returns (row, col) pairs and the total.
pub fn assign_any(n: usize, m: usize, c: &[i64]) -> (Vec<(usize, usize)>, i64) {
    if n <= m {
        let a = hungarian(n, m, c);
        return (a.col_of_row.iter().enumerate().map(|(i, &j)| (i, j)).collect(), a.total);
    }
    let t: Vec<i64> = (0..m * n).map(|k| c[(k % n) * m + k / n]).collect(); // m x n transpose
    let a = hungarian(m, n, &t);
    (a.col_of_row.iter().enumerate().map(|(j, &i)| (i, j)).collect(), a.total)
}

Tested here on 400 random cases (n ≤ 7, m ≤ n + 2, costs −30..69): total equals brute force; u_i + v_j ≤ c_ij everywhere, with equality on assigned pairs; zero-row padding to square gives the same total; hungarian_max and assign_any (n > m) match brute force and the transpose.

Complexity and constants

• O(n² m): each of the n rows runs at most min(i, m) scan passes of O(m). Square: O(n³).
• Memory: n·m·8 B for i64 costs (32 MB at n = 2000) plus O(m) per search. The scan reads one cost row sequentially per pass, so it is bandwidth-bound once the matrix exceeds cache.
• Cost drivers: how many passes each row needs (few for random costs, up to i for ties and structure), and the matrix size.

Efficiency

• Cost: measured on this Mac (Apple M4 Pro), one thread, i64, median of 5 runs of hungarian(n, n, ..).

Tied and structured inputs grow about 8× per doubling (n³); random inputs about 4-6× (roughly n^2.2 from 500 to 2000), presumably because most searches end after a few passes (not profiled). Bound: memory bandwidth over the cost matrix.

• Levers:
  • shortest augmenting paths with potentials instead of re-running matching on the tight graph — O(n⁴) → O(n³) (cp-algorithms) — always
  • Jonker–Volgenant initialisation phases before the search — JV 1987 (SciPy's variant omits them, per its docs; gain not measured here) — dense instances
  • sparse solver after gating away impossible pairs — avoids the n² matrix (SciPy min_weight_full_bipartite_matching, LAPJVsp) — most pairs are impossible
  • split into independent blocks (connected components of allowed pairs) — n³ becomes Σ n_k³ — gated, clustered problems
• Naive → standard → state of the art:
  • all permutations · n! orderings
  • Kuhn's matrix method · O(n⁴); potentials + shortest augmenting path (code above) · O(n³), 217 ms at n = 2000 random, 4.3 s tied (measured)
  • JV family (SciPy's modified JV, lap, lapjv) and sparse LAPJVsp · same dense O(n³) bound; sparse inputs handled; not benchmarked here

Implementation notes

• Rectangular: with n ≤ m the code leaves m − n columns free, same as padding with zero-cost dummy rows (tested). With n > m transpose (assign_any). cp-algorithms pads with zero or infinite rows and columns instead (page).
• Maximise: negate the weights (SORT negates IoU). Costs may be negative; potentials handle them.
• Forbidden pairs: use a large finite cost such as 2⁴⁰ and keep it far below INF; after solving, a total ≥ 2⁴⁰ means no feasible assignment. i64::MAX overflows in c − u − v.
• Duals: u, v are returned; they prove optimality (a checker needs only O(nm)) and give each pair's reduced cost.
• Floats: potentials accumulate rounding; scale to integers or compare with a tolerance.

Pitfalls

• Greedy "pick the cheapest pair, remove row and column" is not optimal: for [[1, 2], [2, 100]] it pays 101, optimal is 4.
• Mixing row-major and column-major indexing in c[i * m + j] silently solves the transposed problem.
• Assuming the optimum is unique: ties give several equal-cost assignments, and SciPy documents that output may vary between versions.

Variants and alternatives

In production

• SORT: association per frame; the reference code calls lap.lapjv, falling back to SciPy's linear_sum_assignment.
• DETR: HungarianMatcher runs linear_sum_assignment per image in the batch, on the CPU.
• SciPy: scipy.optimize.linear_sum_assignment is a modified Jonker–Volgenant (Crouse 2016) and min_weight_full_bipartite_matching is LAPJVsp.
• Google OR-Tools: a dedicated linear sum assignment solver next to MIP and CP-SAT models.

Implementations

Sources

• H. W. Kuhn, The Hungarian method for the assignment problem, Naval Research Logistics Quarterly 2(1-2), 1955. paper
• J. Munkres, Algorithms for the Assignment and Transportation Problems, J. SIAM 5(1), 1957. paper
• R. Jonker, A. Volgenant, A shortest augmenting path algorithm for dense and sparse linear assignment problems, Computing 38(4), 1987. paper
• D. F. Bertsekas, The auction algorithm: a distributed relaxation method for the assignment problem, Annals of Operations Research 14, 1988. paper
• D. F. Crouse, On implementing 2D rectangular assignment algorithms, IEEE Trans. Aerospace and Electronic Systems 52(4), 2016. paper
• A. Bewley, Z. Ge, L. Ott, F. Ramos, B. Upcroft, Simple Online and Realtime Tracking, 2016. paper · code
• N. Carion et al., End-to-End Object Detection with Transformers, ECCV 2020. paper · code
• cp-algorithms, Hungarian algorithm for solving the assignment problem. page
• SciPy documentation, linear_sum_assignment and OR-Tools linear sum assignment solver.