The Algorithm Collection

Who This Vignette Is For

Audience: Researchers curious about optimization algorithms, developers choosing the right method for their problem, anyone wondering why modern algorithms beat the textbook Hungarian method by two orders of magnitude.

Prerequisites:

• Basic understanding of assignment problems (matching rows to columns)

• Familiarity with cost matrices (you want to minimize total cost)

• Comfort with big-O notation (helpful but not required)

What You’ll Learn:

• Why the textbook Hungarian method is slow and how later algorithms improved it

• How primal-dual, auction, and network flow approaches work

• Which algorithm to choose for your problem

• How couplr’s automatic selection works

Hungarian on a 30x30

set.seed(1)
cost_hg <- matrix(sample(1:100, 900, replace = TRUE), 30, 30)
lap_animate(cost_hg, method = "hungarian")

Each row’s processing grows a Dijkstra shortest-path tree on reduced costs from the free row. Orange edges are the scan frontier; a free column popping off the queue terminates the search and the red dashed augmenting path commits one new matched edge before the next free row is selected.

Auction on a 30x30

set.seed(2)
cost_a <- matrix(sample(1:100, 900, replace = TRUE), 30, 30)
lap_animate(cost_a, method = "auction")

Bidders compete: each unassigned row picks its most-valued column and bids (best - second) + eps. When a column already held by another row is grabbed, that row is freed and rejoins the bidding. Watch the matched count climb in fits as displacements ripple through.

Gabow-Tarjan bit-scaling on a 30x30

set.seed(3)
cost_gt <- matrix(sample(1:50, 900, replace = TRUE), 30, 30)
lap_animate(cost_gt, method = "gabow_tarjan")

The header here tracks Bit k / n_bits resolved rather than matched edges, because GT’s progress metric is bit-precision: each phase solves the assignment problem fully at the current bit-precision, then the next phase incorporates one more bit and rebuilds the matching from scratch (with the duals warm-starting from the previous phase). phase_start frames mark each bit being processed (MSB first). Within a bit, Hopcroft-Karp finds maximal augmenting paths in the equality graph; a single 1-feasibility Hungarian step closes any gap before the next bit is added.

Jonker-Volgenant pre-stages on a 30x30

set.seed(4)
cost_jv <- matrix(sample(1:100, 900, replace = TRUE), 30, 30)
lap_animate(cost_jv, method = "jv")

Column reduction (greedy back-to-front), reduction transfer (tighten v on singletons), and augmenting row reduction usually do most of the work; you can watch the matched counter race upward before the Dijkstra main loop runs for any rows still unmatched after the pre-stages.

The same call works for method = "munkres" (the matrix-form 1957 algorithm with star/prime zeros and cover lines) and method = "auction_scaled" (the epsilon-scaling variant that runs an entire auction at progressively tighter eps).

The LP Formulation

The assignment problem is a linear program. Let \(x_{ij} \in \{0,1\}\) indicate whether worker \(i\) is assigned to job \(j\):

\[ \begin{aligned} \min \quad & \sum_{i,j} c_{ij} x_{ij} \\ \text{s.t.} \quad & \sum_j x_{ij} = 1 \quad \forall i \quad \text{(each worker assigned once)} \\ & \sum_i x_{ij} = 1 \quad \forall j \quad \text{(each job filled once)} \\ & x_{ij} \geq 0 \end{aligned} \]

The constraint matrix is totally unimodular: every square submatrix has determinant 0, 1, or -1. This guarantees integer solutions even when we relax \(x_{ij} \in \{0,1\}\) to \(x_{ij} \geq 0\).

Duality: The Key to Efficiency

The dual LP introduces prices \(u_i\) for workers and \(v_j\) for jobs:

\[ \begin{aligned} \max \quad & \sum_i u_i + \sum_j v_j \\ \text{s.t.} \quad & u_i + v_j \leq c_{ij} \quad \forall i,j \end{aligned} \]

Complementary slackness links primal and dual: if \(x_{ij} = 1\) in an optimal assignment, then \(u_i + v_j = c_{ij}\) (the constraint is tight). This means optimal assignments use only tight edges where the dual constraint holds with equality.

The reduced cost of edge \((i,j)\) is \(\bar{c}_{ij} = c_{ij} - u_i - v_j\). An edge is tight when \(\bar{c}_{ij} = 0\).

Different algorithms exploit this duality in different ways.

Hungarian Algorithm (1955)

The algorithm everyone learns. Published by Harold Kuhn, based on work by Hungarian mathematicians Koenig and Egervary.

The idea: Maintain dual prices \((u_i, v_j)\) such that \(u_i + v_j \leq c_{ij}\) for all pairs. Edges where equality holds are “tight”: the only edges that can appear in an optimal solution.

The algorithm:

1. Initialize prices. Find tight edges.

2. Find a maximum matching using only tight edges.

3. If the matching is complete, done.

4. Otherwise, update prices to create new tight edges. Repeat.

Complexity: \(O(n^3)\)

The problem: That \(O(n^3)\) hides a large constant. The price updates and augmenting path searches are expensive. For a 1000x1000 matrix, you might wait 10+ seconds.

cost <- matrix(c(10, 19, 8, 15, 10, 11, 9, 12, 14), nrow = 3, byrow = TRUE)
result <- lap_solve(cost, method = "hungarian")
print(result)
#> Assignment Result
#> =================
#> 
#> # A tibble: 3 × 3
#>   source target  cost
#>    <int>  <int> <dbl>
#> 1      1      3     8
#> 2      2      2    10
#> 3      3      1     9
#> 
#> Total cost: 27 
#> Method: hungarian

Hungarian works. It’s clean and easy to teach. But in 1987, two Dutch researchers found something faster.

Jonker-Volgenant Algorithm (1987)

Roy Jonker and Anton Volgenant asked: what if we start with a good guess and fix it?

The key insight: Column reduction. Before any sophisticated search, greedily assign each row to its cheapest available column. This often gets most of the matching right immediately.

The algorithm:

1. Column reduction: For each column, find the two smallest costs. The difference is the “advantage” of the best row.

2. Reduction transfer: Assign rows to columns, handling conflicts by dual variable updates.

3. Augmentation: For any remaining unmatched rows, use Dijkstra-style shortest path search.

Complexity: Still \(O(n^3)\), with a smaller constant: on the package’s benchmark it takes 0.52 of the shortest-augmenting-path Hungarian’s time at n = 1000.

set.seed(123)
n <- 100
cost <- matrix(runif(n * n, 0, 100), n, n)
result <- lap_solve(cost, method = "jv")
cat("Total cost:", round(get_total_cost(result), 2), "\n")
#> Total cost: 149.09

JV became the de facto standard. For dense problems up to a few thousand rows, it’s hard to beat.

But JV has a limitation: it’s fundamentally serial. Each augmenting path depends on the previous. For very large problems, we need a different approach.

Auction Algorithm (1988)

Dimitri Bertsekas asked: what if we thought of assignment as an economics problem?

The metaphor: Workers are buyers. Jobs are goods. Each job has a price. Workers bid for their favorite jobs. Prices rise when there’s competition. Equilibrium = optimal assignment.

The algorithm:

1. Each unassigned worker finds their best job (highest value minus price).

2. The worker bids: new price = old price + (best value - second-best value) + epsilon.

3. If someone else held that job, they become unassigned.

4. Repeat until everyone is assigned.

Why epsilon matters: Without epsilon, two workers could bid infinitely against each other, each raising the price by 0. The epsilon ensures progress.

It also costs exactness. Bidding stops at an assignment within \(n\epsilon\) of the optimum, which is optimal on integer costs once \(\epsilon < 1/n\) but not on real-valued ones, where two assignments can differ by less than any fixed \(n\epsilon\). couplr’s auctions therefore finish with a repair step: shortest paths over the final prices turn them into exact dual potentials, and any cheaper reassignment found on the way is applied, so the matching returned is optimal.

Complexity: \(O(n^2 \log(nC) / \epsilon)\) where \(C\) is the cost range.

couplr offers three Auction variants:

set.seed(123)
n <- 100
cost <- matrix(runif(n * n, 0, 100), n, n)
result <- lap_solve(cost, method = "auction")
cat("Total cost:", round(get_total_cost(result), 2), "\n")
#> Total cost: 149.09

Auction shines for large dense problems. But it’s sensitive to epsilon. Get it wrong and performance degrades, or the algorithm cycles forever.

The next algorithm makes epsilon-scaling systematic.

Cost-Scaling Algorithm / CSA (1995)

Andrew Goldberg and Robert Kennedy applied the push-relabel method for minimum-cost flow to assignment.

The idea: Start with \(\epsilon\) at the span of the costs. Each phase (a refine) divides epsilon by 10, clears the matching, and discharges the unmatched rows from a stack by double-push: the row takes its cheapest column, the row that held it becomes unmatched, and the column’s price drops to the row’s second-cheapest reduced cost less epsilon. Prices only fall, so a row keeps its three cheapest columns between scans and rescans its full row only when fewer than two of them are still below the fourth-cheapest value it saw. After \(O(\log C)\) phases the assignment is within \(n\epsilon\) of optimal, and a repair step closes the remaining gap.

Why it’s fast: Each phase is cheap because the previous phase’s solution is a good starting point. The algorithm exploits its own progress.

Complexity: \(O(n^3)\) amortized, often faster in practice.

set.seed(456)
n <- 100
cost <- matrix(runif(n * n, 0, 100), n, n)
result <- lap_solve(cost, method = "csa")
cat("Total cost:", round(get_total_cost(result), 2), "\n")
#> Total cost: 192.48

CSA often wins benchmarks for medium-large dense problems. It’s the workhorse.

But there’s an even stranger approach: what if instead of scaling costs, you scaled bits?

Gabow-Tarjan Algorithm (1989)

Harold Gabow and Robert Tarjan developed a clever algorithm based on binary representations. It’s also one of the most complex to implement.

The insight: Integer costs have a natural scale: binary digits. Process costs from most significant to least significant bit. At each scale, solve a simpler problem. Use that solution to warm-start the next scale.

The algorithm (simplified):

1. Initialize at the coarsest scale (most significant bit only).

2. Double the scale: multiply all costs by 2. This “doubles” the current solution’s slack.

3. Restore 1-feasibility (see below).

4. Use Hungarian-style search for augmenting paths.

5. Repeat until all bits are processed.

What is 1-feasibility? Standard dual feasibility requires \(u_i + v_j \leq c_{ij}\) for all edges. 1-feasibility relaxes this: we allow \(u_i + v_j \leq c_{ij} + 1\). The “1” comes from the current bit position. At each scaling phase, we only need reduced costs to be within 1 of optimal. When we refine to the next bit, we tighten the bound. After processing all \(\log C\) bits, the slack is less than 1, which for integers means exactly 0: true optimality.

Complexity: \(O(\sqrt{n}\, m \log(nC))\) where \(C\) is the maximum cost and \(m\) is the number of finite edges. For dense graphs with \(m = n^2\), this is \(O(n^{5/2} \log(nC))\) — sub-cubic in \(n\), which is why Gabow and Tarjan’s result was a breakthrough. It improves on the earlier \(O(n^{3/4} m \log C)\) scaling bound of Gabow (1985). The improvement over the \(O(n^3 \log C)\) bound for plain cost-scaling comes from finding a maximal vertex-disjoint set of augmenting paths at each bit, rather than one path at a time.

Rarely seen outside academic papers. The bookkeeping across scaling phases is complex enough that most implementations skip it.

set.seed(42)
n <- 50
# Use integer costs with large range - Gabow-Tarjan's strength
cost <- matrix(sample(1:100000, n * n, replace = TRUE), n, n)
result <- lap_solve(cost, method = "gabow_tarjan")
cat("Total cost:", get_total_cost(result), "\n")
#> Total cost: 151632

Gabow-Tarjan is primarily of theoretical interest. It provides the best known worst-case bounds for integer costs. That closes the scaling family, and there’s a completely different way to think about the problem entirely.

Network Simplex

The simplex method, specialized for networks. The key insight: for network flow problems, the simplex basis corresponds to a spanning tree of the graph. This makes basis operations (pivoting) much faster than general LP.

The algorithm:

1. Initialize: Find a spanning tree \(T\) that supports a feasible flow (e.g., all flow through a single hub).

2. Compute potentials: For tree edges, set \(\pi_i - \pi_j = c_{ij}\). This determines all node potentials uniquely (up to a constant).

3. Price non-tree edges: For each edge \((i,j) \notin T\), compute reduced cost \(\bar{c}_{ij} = c_{ij} - \pi_i + \pi_j\).

4. Pivot: If any non-tree edge has \(\bar{c}_{ij} < 0\), adding it to \(T\) creates a cycle. Push flow around the cycle until some tree edge saturates. Remove that edge; the non-tree edge enters the basis.

5. Repeat until all reduced costs are non-negative (optimality).

Why trees?: In a tree, there’s exactly one path between any two nodes. This means: - Flow is uniquely determined by supplies/demands

• Potentials can be computed in \(O(n)\) by tree traversal

• Each pivot changes \(O(n)\) potentials (along a tree path), not \(O(n^2)\)

Complexity: \(O(n^3)\) typical, strongly polynomial with anti-cycling rules.

Best for: Problems where you need dual variables for sensitivity analysis. Problems already formulated as network flows. Cases where you want guaranteed finite convergence.

set.seed(789)
n <- 100
cost <- matrix(runif(n * n, 0, 100), n, n)
result <- lap_solve(cost, method = "network_simplex")
cat("Total cost:", round(get_total_cost(result), 2), "\n")
#> Total cost: 156.88

Network Simplex is a standard tool in operations research. Not always the fastest, but reliable and provides rich dual information.

Push-Relabel Algorithm

Goldberg and Tarjan’s push-relabel algorithm (1988), adapted for minimum-cost flow.

The key idea: Traditional algorithms maintain a valid flow and search for augmenting paths. Push-relabel inverts this: maintain a preflow (flow that may violate conservation at intermediate nodes) and work locally to eliminate violations.

Key concepts:

• Excess: \(e(v) = \text{flow in} - \text{flow out}\). Preflow allows \(e(v) \geq 0\) at non-sink nodes.

• Height function: \(h: V \to \mathbb{Z}_{\geq 0}\) with \(h(t) = 0\) and \(h(u) \leq h(v) + 1\) for residual edges \((u,v)\).

• Admissible edge: Residual edge \((u,v)\) where \(h(u) = h(v) + 1\) (flow goes “downhill”).

The algorithm:

1. Initialize: Saturate all edges from source. Set \(h(s) = n\), \(h(v) = 0\) for \(v \neq s\).

2. While any node \(v\) has excess \(e(v) > 0\):

   • Push: If admissible edge \((v, w)\) exists, push \(\min(e(v), \text{capacity})\) flow along it.

   • Relabel: If no admissible edge, set \(h(v) = 1 + \min\{h(w) : (v,w) \text{ is residual}\}\).

3. When no excess remains at intermediate nodes, we have a valid maximum flow.

For minimum-cost flow: Modify to only push along edges with zero reduced cost, and relabel using potential updates. This gives the cost-scaling push-relabel variant.

Complexity: \(O(n^2 m)\) for max-flow, \(O(n^3 \log(nC))\) for min-cost flow.

Strengths: Highly parallelizable since pushes are local operations. Excellent cache behavior. Dominates in practice for max-flow; competitive for min-cost flow on dense graphs.

set.seed(222)
n <- 100
cost <- matrix(runif(n * n, 0, 100), n, n)
result <- lap_solve(cost, method = "push_relabel")
cat("Total cost:", round(get_total_cost(result), 2), "\n")
#> Total cost: 169.7

Two network perspectives. Same problem. Different algorithmic approaches.

But all these algorithms assume dense, square matrices. Real problems are messier.

HK01: Binary Costs

When costs are only 0 or 1, we don’t need the full machinery.

The algorithm: Hopcroft-Karp for maximum cardinality matching, run on zero-cost edges first. Then add 1-cost edges as needed.

Complexity: \(O(n^{2.5})\) for binary costs.

set.seed(101)
n <- 100
cost <- matrix(sample(0:1, n^2, replace = TRUE, prob = c(0.3, 0.7)), n, n)
result <- lap_solve(cost, method = "hk01")
cat("Total cost:", get_total_cost(result), "\n")
#> Total cost: 0

When you have binary costs and large \(n\), HK01 is dramatically faster.

SAP and LAPMOD: Sparse Problems

When 80% of entries are forbidden (Inf or NA), why store them?

SAP (Shortest Augmenting Path) and LAPMOD use sparse representations: adjacency lists instead of dense matrices.

Complexity: \(O(n^2 + nm)\) where \(m\) is the number of allowed edges.

set.seed(789)
n <- 100
cost <- matrix(Inf, n, n)
edges <- sample(1:(n^2), floor(0.2 * n^2))  # Only 20% allowed
cost[edges] <- runif(length(edges), 0, 100)
result <- lap_solve(cost, method = "sap")
cat("Total cost:", round(get_total_cost(result), 2), "\n")
#> Total cost: 794.94

Sparse storage is what these two are written for, but on the regime grid behind the dispatcher neither "sap" nor "lapmod" was often quicker than "jv" on sparse matrices, so method = "auto" sends sparse input to "jv" and both stay available by name.

"sap_dense" runs the same successive-shortest-path search with the opposite storage assumption. Its Dijkstra step picks the next column by scanning every unfinalized one instead of pulling from a heap, which costs \(O(m)\) per extraction and \(O(n m^2)\) over a full solve. On a dense cost matrix, where every column is a candidate at every step, the scan does the work a heap would do anyway without paying for the heap.

set.seed(789)
n <- 60
cost <- matrix(runif(n * n, 0, 100), n, n)
result <- lap_solve(cost, method = "sap_dense")
cat("Total cost:", round(get_total_cost(result), 2), "
")
#> Total cost: 153.98

Ramshaw-Tarjan: Rectangular Problems (2012)

Most algorithms assume square matrices. When you have \(n\) workers and \(m > n\) jobs, standard approaches pad with \(m - n\) dummy workers at zero cost. This works but wastes effort: the algorithm processes \(m \times m\) entries when only \(n \times m\) matter.

Ramshaw and Tarjan (2012) developed an algorithm that handles rectangularity natively by exploiting the structure of unbalanced bipartite graphs.

The key insight: In a rectangular assignment, we match all \(n\) rows but only \(n\) of the \(m\) columns. The dual problem has different structure: row duals \(u_i\) are unconstrained, but column duals \(v_j\) must satisfy \(v_j \leq 0\) for unmatched columns.

The algorithm:

1. Maintain dual variables \((u, v)\) with \(u_i + v_j \leq c_{ij}\) and \(v_j \leq 0\) for free columns.

2. Use a modified Dijkstra search that respects the asymmetric dual constraints.

3. When augmenting, update duals to preserve feasibility without padding.

Complexity: \(O(nm \log n)\) using Fibonacci heaps, or \(O(nm + n^2 \log n)\) with simpler structures.

For highly rectangular problems (e.g., matching 100 treatments to 10,000 controls), this avoids the \(O(m^3)\) cost of padding to square.

set.seed(333)
n_rows <- 30
n_cols <- 100  # Highly rectangular: 30 x 100
cost <- matrix(runif(n_rows * n_cols, 0, 100), n_rows, n_cols)
result <- lap_solve(cost, method = "ramshaw_tarjan")
cat("Matched", sum(result$assignment > 0), "of", n_rows, "rows\n")
#> Warning: Unknown or uninitialised column: `assignment`.
#> Matched 0 of 30 rows

When you have significantly more columns than rows (or vice versa), Ramshaw-Tarjan avoids the wasted work of padding. The newest algorithm in couplr’s collection, and essential for large-scale matching problems where treatment and control pools have very different sizes.

K-Best Solutions (Murty’s Algorithm)

What if you want the 2nd best assignment? The 3rd best? The k-th best?

cost <- matrix(c(10, 19, 8, 15, 10, 18, 7, 17, 13, 16, 9, 14, 12, 19, 8, 18),
               nrow = 4, byrow = TRUE)
kbest <- lap_solve_kbest(cost, k = 5)
summary(kbest)
#> # A tibble: 5 × 4
#>    rank solution_id total_cost n_assignments
#>   <int>       <int>      <dbl>         <int>
#> 1     1           1         49             4
#> 2     2           2         50             4
#> 3     3           3         50             4
#> 4     4           4         50             4
#> 5     5           5         51             4

Use cases: Sensitivity analysis. Alternative plans when the optimal is infeasible. Understanding how costs affect solutions.

Bottleneck Assignment

Minimize the maximum edge cost instead of the sum.

cost <- matrix(c(5, 9, 2, 10, 3, 7, 8, 4, 6), nrow = 3, byrow = TRUE)
result <- bottleneck_assignment(cost)
cat("Bottleneck (max edge):", result$bottleneck, "\n")
#> Bottleneck (max edge): 6

Use cases: Load balancing. Fairness constraints. Worst-case optimization.

Sinkhorn: Soft Assignment

Entropy-regularized optimal transport. Instead of hard 0/1 assignment, produce a doubly-stochastic transport plan.

cost <- matrix(c(1, 2, 3, 4), nrow = 2)
result <- sinkhorn(cost, lambda = 10)
print(round(result$transport_plan, 3))
#>      [,1] [,2]
#> [1,] 0.25 0.25
#> [2,] 0.25 0.25

Use cases: Probabilistic matching. Domain adaptation. Wasserstein distances.

Dual Variables

Extract dual prices for sensitivity analysis.

cost <- matrix(c(10, 19, 8, 15, 10, 18, 7, 17, 13), nrow = 3, byrow = TRUE)
result <- assignment_duals(cost)
cat("Row duals (u):", result$u, "\n")
#> Row duals (u): 3 8 8
cat("Col duals (v):", result$v, "\n")
#> Col duals (v): -1 2 5

Use cases: Shadow prices. Identifying critical assignments. Marginal cost analysis.