3. Network Optimization

Many OR problems are about moving something through a network: goods through a supply chain, data through routers, vehicles through roads, workers onto jobs. Network optimization models these problems with graphs and solves them with specialized algorithms. These algorithms are much faster than a general LP solver, and their integer data always gives integer solutions. This chapter builds the graph vocabulary, then covers the classic problem families (shortest paths, spanning trees, maximum flow, minimum cost flow, transportation and assignment), and shows how each one is a special case of a linear program with special structure.

3.1 Graph Theory Fundamentals

This subchapter introduces the objects every network algorithm works on (graphs, paths, trees and their matrix representations) and explains total unimodularity, the property that makes network LPs return integer solutions.

Most examples in this chapter use the same small delivery network. A company ships goods from a depot $s$ to a customer $t$ through four hubs $a, b, c, d$. Each arc carries a unit shipping cost and a daily capacity:

  Delivery network, arc label = (cost per unit, capacity)

                        (7,6)
              a --------------------> c
            ^ |                       ^ \
    (2,10) /  |                       |  \ (1,8)
          /   | (1,5)          (2,4)  |   v
        s     |                       |     t
          \   |                       |   ^
    (4,8)  \  |                       |  / (6,7)
            v v                       | /
              b --------------------> d
                        (3,7)

Arc	$s \to a$	$s \to b$	$a \to b$	$a \to c$	$b \to d$	$c \to t$	$d \to c$	$d \to t$
Cost $c_{ij}$	2	4	1	7	3	1	2	6
Capacity $u_{ij}$	10	8	5	6	7	8	4	7

Graphs, Digraphs and Networks

Why it matters. Before you can optimize a network you need a precise language for it. The same few definitions describe road maps, computer networks, project schedules and job assignments, so learning them once lets you reuse every algorithm in this chapter across all of those domains.

Theory. A graph $G = (V, E)$ is a set of nodes (vertices) $V$ and a set of edges $E$, each edge joining two nodes. We write $n = |V|$ and $m = |E|$.

In an undirected graph, an edge ${i, j}$ can be traversed both ways (a two-way street, a cable).
In a directed graph (digraph), an arc $(i, j)$ goes from its tail $i$ to its head $j$ (a one-way street, a pipe with a valve).
A graph is simple if it has no loops $(i, i)$ and no parallel edges. Most algorithms assume this, and real models are easily converted (merge parallel arcs or split them with an extra node).
The degree of a node counts its incident edges. In a digraph, the out-degree and in-degree count arcs leaving and entering it. Every edge has two ends, so $\sum_i \deg(i) = 2m$ (handshake lemma). Every arc has one tail and one head, so $\sum_i \deg^+(i) = \sum_i \deg^-(i) = m$.
A simple digraph has at most $n(n-1)$ arcs. Graphs with $m$ close to $n^2$ are dense. Graphs with $m$ close to $n$ are sparse, and almost all real networks (roads, power grids) are sparse.

A network is a digraph with numbers attached. These are the data that turn a picture into an optimization model:

Symbol	Attached to	Meaning	Example
$c_{ij}$	arc	cost, length or time per unit of flow	shipping cost per pallet
$u_{ij}$	arc	capacity (upper bound on flow)	trucks per day on a route
$l_{ij}$	arc	lower bound on flow (often 0)	contractual minimum volume
$b_i$	node	supply ($b_i > 0$), demand ($b_i < 0$) or transshipment ($b_i = 0$)	factory output, customer order

Application	Nodes	Arcs	Typical question
Road network	intersections	road segments	fastest route (shortest path)
Telecom / Internet	routers	links	maximum throughput (max flow)
Supply chain	plants, warehouses, customers	shipping lanes	cheapest distribution (min cost flow)
Staffing	workers and jobs	"worker can do job"	best assignment (matching)
Project	events	activities	earliest completion (longest path in a DAG)

Representations. An algorithm needs fast access to the arcs leaving a node. The adjacency list stores, for each node, a list of its outgoing arcs. It uses $O(n + m)$ memory and is the default for sparse networks. The adjacency matrix and incidence matrix are covered later in this subchapter.

from collections import defaultdict

# Delivery network: (tail, head, cost per unit, capacity)
ARCS = [("s", "a", 2, 10), ("s", "b", 4, 8), ("a", "b", 1, 5), ("a", "c", 7, 6),
        ("b", "d", 3, 7), ("c", "t", 1, 8), ("d", "c", 2, 4), ("d", "t", 6, 7)]
NODES = ["s", "a", "b", "c", "d", "t"]

succ = defaultdict(list)                  # adjacency list: node -> outgoing arcs
out_deg, in_deg = defaultdict(int), defaultdict(int)
for i, j, cost, cap in ARCS:
    succ[i].append((j, cost, cap))
    out_deg[i] += 1
    in_deg[j] += 1

for v in NODES:
    print(f"{v}: out={out_deg[v]} in={in_deg[v]}  successors {[j for j, *_ in succ[v]]}")
n, m = len(NODES), len(ARCS)
print("sum of out-degrees =", sum(out_deg.values()), "= m")
print(f"density m / (n(n-1)) = {m / (n * (n - 1)):.2f}")
# s: out=2 in=0  successors ['a', 'b']      <- a pure source
# a: out=2 in=1  successors ['b', 'c']
# b: out=1 in=2  successors ['d']
# c: out=1 in=2  successors ['t']
# d: out=2 in=1  successors ['c', 't']
# t: out=0 in=2  successors []              <- a pure sink
# sum of out-degrees = 8 = m
# density m / (n(n-1)) = 0.27

Key takeaways

A graph is $G = (V, E)$. Edges are undirected, arcs are directed from tail to head.
A network is a digraph plus data: arc costs $c_{ij}$, capacities $u_{ij}$, lower bounds $l_{ij}$ and node supplies $b_i$.
Degree sums: $\sum \deg = 2m$ in undirected graphs and $\sum \deg^+ = \sum \deg^- = m$ in digraphs.
Real networks are sparse, so adjacency lists ($O(n + m)$ memory) are the standard representation.

Practice

(Beginner) A simple undirected graph has 7 nodes with degrees 3, 3, 2, 2, 2, 1, 1. How many edges does it have? Could a simple graph have degrees 3, 3, 3, 2, 2?

(Intermediate) Model a city with three one-way streets and two two-way streets as a digraph. How many arcs does it have, and how do you represent a two-way street with different travel times in each direction?

(Interview) "You need to store a social network with 1 billion users and about 200 friends per user. Adjacency matrix or adjacency list?" Hint: compare $n^2$ with $n + m$ in bytes, then think about which operations (neighbor iteration vs edge lookup) the application performs most.

Paths, Cycles and Connectivity

Why it matters. Almost every network question is about paths: is there a route, which route is cheapest, how many routes are disjoint? Connectivity tells you whether a problem has a solution at all (no path from $s$ to $t$ means no shipment), and cycles decide which algorithms apply (a network without directed cycles admits much faster methods).

Theory.

Term	Definition
Walk	Sequence of nodes $v_0, v_1, \dots, v_k$ with each consecutive pair joined by an edge (arc, in a digraph)
Path	A walk with no repeated nodes
Cycle	A walk with $v_0 = v_k$, $k \ge 1$ (at least 3 in a simple undirected graph), and no other repeated nodes
Directed path / cycle	Every arc is traversed from tail to head
Reachable	$j$ is reachable from $i$ if a directed path $i \to j$ exists
Connected (undirected)	Every pair of nodes is joined by a path
Strongly connected (digraph)	Every node is reachable from every other node
Weakly connected (digraph)	Connected after ignoring arc directions
DAG	Directed acyclic graph: a digraph with no directed cycle
Cut $(S, V \setminus S)$	Partition of the nodes. $\delta^+(S)$ = arcs leaving $S$, $\delta^-(S)$ = arcs entering $S$

Think of a digraph as a city of one-way streets. "Weakly connected" means the map is in one piece. "Strongly connected" means you can drive from anywhere to anywhere without breaking traffic rules.

Two search procedures answer most connectivity questions in $O(n + m)$ time:

Breadth-first search (BFS) explores nodes in layers of increasing arc count, so the first time it reaches a node it has found a path with the fewest arcs. Edmonds-Karp (Section 3.4) is built on this.
Depth-first search (DFS) goes as deep as possible before backtracking. A node that reaches a node still on the DFS stack (a back arc) reveals a directed cycle. DFS is also the basis of topological sorting and strongly connected components.

A cut separates the network. If $s \in S$ and $t \notin S$, every $s$-$t$ path must use at least one arc of $\delta^+(S)$. This simple fact is the seed of the max-flow min-cut theorem (Section 3.4).

Example: reachability and cycles in the delivery network.

from collections import deque

ARCS = [("s", "a"), ("s", "b"), ("a", "b"), ("a", "c"),
        ("b", "d"), ("c", "t"), ("d", "c"), ("d", "t")]

def successors(arcs):
    succ = {}
    for i, j in arcs:
        succ.setdefault(i, []).append(j)
        succ.setdefault(j, [])
    return succ

def bfs_path(succ, source, target):
    """Path with the fewest arcs (or None) and the set of reachable nodes."""
    pred = {source: None}
    queue = deque([source])
    while queue:
        v = queue.popleft()
        for w in succ[v]:
            if w not in pred:             # first visit = fewest arcs
                pred[w] = v
                queue.append(w)
    if target not in pred:
        return None, set(pred)
    path, v = [], target
    while v is not None:
        path.append(v)
        v = pred[v]
    return path[::-1], set(pred)

def find_cycle(succ):
    """Return a directed cycle as a node list, or None (DFS with 3 colors)."""
    WHITE, GREY, BLACK = 0, 1, 2          # unvisited / on the DFS stack / finished
    color = {v: WHITE for v in succ}
    stack = []
    def dfs(v):
        color[v] = GREY
        stack.append(v)
        for w in succ[v]:
            if color[w] == GREY:          # back arc -> cycle w ... v -> w
                return stack[stack.index(w):] + [w]
            if color[w] == WHITE and (cyc := dfs(w)):
                return cyc
        color[v] = BLACK
        stack.pop()
        return None
    for v in succ:
        if color[v] == WHITE and (cyc := dfs(v)):
            return cyc
    return None

succ = successors(ARCS)
print("fewest-arc path s->t:", bfs_path(succ, "s", "t")[0])
print("reachable from c:    ", sorted(bfs_path(succ, "c", "s")[1]))
print("directed cycle:      ", find_cycle(succ))
print("after adding c->a:   ", find_cycle(successors(ARCS + [("c", "a")])))
# fewest-arc path s->t: ['s', 'a', 'c', 't']
# reachable from c:     ['c', 't']          <- not strongly connected
# directed cycle:       None                <- the delivery network is a DAG
# after adding c->a:    ['a', 'b', 'd', 'c', 'a']

The fewest-arc path $s$-$a$-$c$-$t$ is not the cheapest one (cost 10 versus 9 for $s$-$a$-$b$-$d$-$c$-$t$). Counting arcs and adding costs are different problems, which is why Section 3.2 needs other algorithms.

Key takeaways

Paths have no repeated nodes, and cycles return to their start. In digraphs, both must respect arc directions.
Undirected graphs are connected or not. Digraphs can be weakly or strongly connected.
BFS finds fewest-arc paths and DFS finds cycles, both in $O(n + m)$.
Every $s$-$t$ path crosses every cut that separates $s$ from $t$.

Practice

(Beginner) List all directed $s$-$t$ paths in the delivery network and their costs. Which has the fewest arcs, and which is cheapest?

(Intermediate) Write a function that returns the strongly connected components of a digraph (Kosaraju: DFS finishing order, then DFS on the reversed graph). Test it on the delivery network plus arc $c \to a$.

(Interview) "A package manager must install dependencies before the packages that need them. How do you detect an impossible (circular) dependency set?" Hint: model packages as nodes and "needs" as arcs, then look for a back arc during DFS.

Trees and Spanning Trees

Why it matters. Trees are the minimal structures that keep a network connected. They appear everywhere in this chapter: the output of a shortest path algorithm is a tree of routes, a minimum spanning tree is the cheapest way to connect sites, and the bases of the network simplex method are spanning trees.

Theory. A tree is a connected undirected graph with no cycles. A forest is a graph with no cycles (a collection of trees). For a graph $T$ with $n$ nodes, the following are equivalent:

$T$ is connected and acyclic.
$T$ is connected and has exactly $n - 1$ edges.
$T$ is acyclic and has exactly $n - 1$ edges.
Every pair of nodes is joined by exactly one path.
$T$ is minimally connected: removing any edge disconnects it.
$T$ is maximally acyclic: adding any edge creates exactly one cycle.

Properties 5 and 6 are the basis of the exchange arguments used in Section 3.3: adding a non-tree edge creates a unique cycle, and removing any edge of that cycle gives another spanning tree.

A spanning tree of a connected graph $G$ is a tree that uses all nodes of $G$ and a subset of its edges. A rooted tree picks one node as root. Every other node then has a unique parent (the next node on its path to the root), so a whole tree fits in one pred array. This is how shortest path algorithms return their routes.

  A graph and one of its spanning trees (n = 6, n - 1 = 5 edges)

  A ---- B ---- D ---- F            A ---- B      D      F
  |    /        |    /              |             |    /
  |  /          |  /                |             |  /
  C ----------- E                   C ----------- E

How many spanning trees can a graph have? For the complete graph $K_n$, Cayley's formula gives $n^{n-2}$, which grows faster than exponentially. For any graph, Kirchhoff's matrix-tree theorem says the count equals any cofactor of the Laplacian $L = D - A$ (degree matrix minus adjacency matrix): delete one row and the matching column and take the determinant. The example checks both against brute force.

import itertools
import numpy as np

def is_spanning_tree(n, edges):
    """n-1 edges and no cycle  <=>  spanning tree (union-find cycle check)."""
    if len(edges) != n - 1:
        return False
    parent = list(range(n))
    def find(v):
        while parent[v] != v:
            v = parent[v]
        return v
    for u, v in edges:
        ru, rv = find(u), find(v)
        if ru == rv:                      # u and v already connected -> cycle
            return False
        parent[ru] = rv
    return True

def count_brute_force(n, edges):
    return sum(is_spanning_tree(n, list(T)) for T in itertools.combinations(edges, n - 1))

def count_matrix_tree(n, edges):
    """Kirchhoff: number of spanning trees = any cofactor of the Laplacian L = D - A."""
    L = np.zeros((n, n))
    for u, v in edges:
        L[u, u] += 1; L[v, v] += 1
        L[u, v] -= 1; L[v, u] -= 1
    return round(np.linalg.det(L[1:, 1:]))   # delete row 0 and column 0

K4 = [(u, v) for u in range(4) for v in range(u + 1, 4)]
# The graph drawn above (used again in Section 3.3): A=0, B=1, C=2, D=3, E=4, F=5
G = [(0, 1), (0, 2), (1, 2), (1, 3), (2, 4), (3, 4), (3, 5), (4, 5)]

for name, n, E in [("K4", 4, K4), ("G", 6, G)]:
    print(f"{name}: brute force = {count_brute_force(n, E)}, "
          f"matrix-tree = {count_matrix_tree(n, E)}, "
          f"subsets tried = {len(list(itertools.combinations(E, n - 1)))}")
# K4: brute force = 16, matrix-tree = 16, subsets tried = 20    <- Cayley: 4^2 = 16
# G: brute force = 30, matrix-tree = 30, subsets tried = 56

Enumerating trees is hopeless for real networks ($K_{20}$ has about $2.6 \times 10^{23}$ spanning trees), which is why Section 3.3 finds the best one greedily instead.

Key takeaways

A tree on $n$ nodes is connected, acyclic and has exactly $n - 1$ edges. Any two of these properties imply the third.
Adding an edge to a tree creates exactly one cycle, and removing an edge splits it into two components.
Rooted trees are stored as parent (predecessor) arrays. Shortest path and search algorithms return them.
The number of spanning trees is huge ($n^{n-2}$ for $K_n$), so optimization over trees needs clever algorithms, not enumeration.

Practice

(Beginner) Draw all 16 spanning trees of $K_4$. How many of them are paths, and how many are stars?

(Intermediate) Prove that a tree with $n \ge 2$ nodes has at least two leaves (nodes of degree 1). Hint: consider a longest path.

(Interview) "You are given $n$ nodes and $n - 1$ edges. How do you check in near-linear time whether they form a tree?" Hint: union-find detects a cycle as soon as an edge joins two nodes already in the same component.

Incidence and Adjacency Matrices

Why it matters. To connect networks with linear programming you need them in matrix form. The node-arc incidence matrix turns "flow in equals flow out" into the linear system $Nx = b$, and its structure explains why network LPs are so well behaved. The adjacency matrix connects graphs with linear algebra in another way: its powers count walks.

Theory.

The node-arc incidence matrix $N \in {-1, 0, 1}^{n \times m}$ has one row per node and one column per arc. Column $k$ for arc $(i, j)$ has $+1$ in row $i$ (the tail), $-1$ in row $j$ (the head) and zeros elsewhere. Some textbooks use the opposite signs, so check the convention.
Row $i$ of $Nx$ is (flow out of $i$) $-$ (flow into $i$). Flow conservation for all nodes is therefore $Nx = b$, with $b_i$ the supply of node $i$.
Every column sums to 0, so the rows are linearly dependent: $\text{rank}(N) \le n - 1$. If the graph is connected, $\text{rank}(N) = n - 1$ exactly, and in every network LP one conservation equation is redundant.
The columns of a cycle are linearly dependent (add them with $\pm 1$ signs along the cycle and you get 0). The columns of a forest are independent. So a set of $n - 1$ columns is a basis of the column space exactly when its arcs form a spanning tree. This is the link between the simplex method and trees (Section 3.5).
The adjacency matrix $A \in {0, 1}^{n \times n}$ has $A_{ij} = 1$ if arc $(i, j)$ exists. Entry $(A^k)_{ij}$ counts the walks from $i$ to $j$ with exactly $k$ arcs.

Representation	Size	Arc lookup $(i, j)$	List neighbors of $i$	Best for
Adjacency list	$O(n + m)$	$O(\deg(i))$	$O(\deg(i))$	algorithms on sparse graphs
Adjacency matrix	$O(n^2)$	$O(1)$	$O(n)$	dense graphs, algebraic tricks
Incidence matrix	$O(nm)$ (sparse: $O(m)$)	$O(m)$	$O(m)$	LP formulations

Example: the delivery network as matrices.

import numpy as np

NODES = ["s", "a", "b", "c", "d", "t"]
ARCS = [("s", "a"), ("s", "b"), ("a", "b"), ("a", "c"),
        ("b", "d"), ("c", "t"), ("d", "c"), ("d", "t")]
idx = {v: k for k, v in enumerate(NODES)}
n, m = len(NODES), len(ARCS)

N = np.zeros((n, m), dtype=int)           # node-arc incidence matrix
Adj = np.zeros((n, n), dtype=int)         # adjacency matrix
for k, (i, j) in enumerate(ARCS):
    N[idx[i], k] = 1                      # arc k leaves i
    N[idx[j], k] = -1                     # arc k enters j
    Adj[idx[i], idx[j]] = 1

print("    " + " ".join(f"{i}{j}" for i, j in ARCS))
for v in NODES:
    print(f"{v}: ", " ".join(f"{x:2d}" for x in N[idx[v]]))
print("column sums:", N.sum(axis=0), " rank:", np.linalg.matrix_rank(N), "= n - 1")

# A flow of 3 units along s-a-c-t: N x gives the net outflow of every node
x = np.array([3, 0, 0, 3, 0, 3, 0, 0])
print("N x =", N @ x, " (supply 3 at s, demand 3 at t, balance elsewhere)")

# (Adj^k)[s, t] counts the s->t walks with exactly k arcs
walks = [int(np.linalg.matrix_power(Adj, k)[idx["s"], idx["t"]]) for k in range(1, n)]
print("s->t walks with k = 1..5 arcs:", walks, " total:", sum(walks))
#     sa sb ab ac bd ct dc dt
# s:   1  1  0  0  0  0  0  0
# a:  -1  0  1  1  0  0  0  0
# b:   0 -1 -1  0  1  0  0  0
# c:   0  0  0 -1  0  1 -1  0
# d:   0  0  0  0 -1  0  1  1
# t:   0  0  0  0  0 -1  0 -1
# column sums: [0 0 0 0 0 0 0 0]  rank: 5 = n - 1
# N x = [ 3  0  0  0  0 -3]  (supply 3 at s, demand 3 at t, balance elsewhere)
# s->t walks with k = 1..5 arcs: [0, 0, 2, 2, 1]  total: 5

The network is a DAG, so every walk is a path, and the total of 5 matches the five $s$-$t$ paths listed in the previous topic's exercise.

Key takeaways

Incidence matrix: column of arc $(i, j)$ has $+1$ at $i$ and $-1$ at $j$. Flow conservation is $Nx = b$.
For a connected graph, $\text{rank}(N) = n - 1$, so one conservation constraint is always redundant.
Sets of $n - 1$ linearly independent incidence columns are exactly the spanning trees.
$(A^k)_{ij}$ counts $k$-arc walks from $i$ to $j$. Use adjacency lists for algorithms and matrices for formulations and algebra.

Practice

(Beginner) Write the incidence matrix of the directed cycle $1 \to 2 \to 3 \to 1$. Find a nonzero vector $x$ with $Nx = 0$ and interpret it as a flow.

(Intermediate) For an undirected graph, orient each edge arbitrarily and build $N$. Show numerically that $N N^\top$ equals the Laplacian $D - A$, whatever orientation you chose.

(Interview) "How would you count the number of distinct routes of length at most $k$ between two routers using matrix operations, and what is the catch if the network has cycles?" Hint: sum $A^1, \dots, A^k$, but remember that powers count walks, which may revisit routers.

Total Unimodularity

Why it matters. Integer programs are hard in general (Chapter 4), yet network problems with integer data are solved by plain LP and still return integer solutions: whole trucks, whole workers, whole units. Total unimodularity is the reason. It explains why the LP relaxation of a network problem is exact, and why the same trick fails for problems like general matching or the traveling salesman problem.

Theory. A matrix $A$ is totally unimodular (TU) if every square submatrix has determinant $0$, $+1$ or $-1$. In particular, every entry is in ${-1, 0, 1}$.

Theorem (Hoffman-Kruskal). If $A$ is TU and $b$ is integer, every vertex of ${x : Ax \le b,\ x \ge 0}$ (or ${x : Ax = b,\ x \ge 0}$) is integer. Hence the simplex method, which returns a vertex, returns an integer optimal solution.

Why. A vertex is a basic solution $x_B = B^{-1} b$. By Cramer's rule, $B^{-1} = \operatorname{adj}(B) / \det(B)$. The adjugate has integer entries, and $\det(B) = \pm 1$ because $B$ is a square submatrix of a TU matrix. So $x_B$ is integer. This is the same argument as in Section 2.2, with integrality added.

Recognizing TU matrices. Checking all submatrices is exponential, but a sufficient condition (Heller-Tompkins) covers network matrices. $A$ is TU if:

every entry is $0$, $+1$ or $-1$,
every column has at most two nonzero entries, and
the rows can be split into two groups $R_1, R_2$ such that a column with two nonzeros of the same sign has them in different groups, and a column with two nonzeros of opposite signs has them in the same group.

Matrix	Why the condition holds	Consequence
Incidence matrix of any digraph	each column has one $+1$ and one $-1$, so put all rows in $R_1$	shortest path, max flow, min cost flow LPs are integral
Incidence matrix of a bipartite undirected graph	each column has two $+1$, one in each side, so $R_1$ = left nodes, $R_2$ = right nodes	transportation, assignment, bipartite matching are integral
Incidence matrix of a graph with an odd cycle	fails: the triangle has determinant $\pm 2$	the matching LP can return halves

TU is preserved by transposing, multiplying a row or column by $-1$, duplicating rows, and appending identity columns (slack variables). So adding slacks or capacity bounds $x \le u$ keeps the network structure intact.

Example: when the LP is (and is not) integral.

import itertools
import numpy as np
from scipy.optimize import linprog

def is_totally_unimodular(A):
    """Brute force: every square submatrix has determinant in {-1, 0, 1}."""
    m, n = A.shape
    for k in range(1, min(m, n) + 1):
        for rows in itertools.combinations(range(m), k):
            for cols in itertools.combinations(range(n), k):
                d = round(np.linalg.det(A[np.ix_(rows, cols)]))
                if d not in (-1, 0, 1):
                    return False, (rows, cols, d)
    return True, None

def incidence(n, edges, directed):
    M = np.zeros((n, len(edges)))
    for k, (i, j) in enumerate(edges):
        M[i, k] = 1
        M[j, k] = -1 if directed else 1   # undirected: both endpoints get +1
    return M

triangle = [(0, 1), (1, 2), (0, 2)]
square = [(0, 1), (1, 2), (2, 3), (0, 3)]   # even cycle -> bipartite
print("directed triangle  :", is_totally_unimodular(incidence(3, triangle, True)))
print("undirected square  :", is_totally_unimodular(incidence(4, square, False)))
print("undirected triangle:", is_totally_unimodular(incidence(3, triangle, False)))

# Matching LP: max sum x_e  s.t. each vertex covered at most once, 0 <= x <= 1
for name, n, E in [("square", 4, square), ("triangle", 3, triangle)]:
    A = incidence(n, E, directed=False)
    res = linprog(-np.ones(len(E)), A_ub=A, b_ub=np.ones(n), bounds=(0, 1))
    print(f"matching LP on {name:8s}: x = {np.round(res.x, 2) + 0.0}, value = {-res.fun:.1f}")
# directed triangle  : (True, None)
# undirected square  : (True, None)
# undirected triangle: (False, ((0, 1, 2), (0, 1, 2), 2))   <- det = 2
# matching LP on square  : x = [1. 0. 1. 0.], value = 2.0   <- integral
# matching LP on triangle: x = [0.5 0.5 0.5], value = 1.5   <- fractional

On the triangle, the best real matching has 1 edge but the LP claims 1.5. Closing that gap needs extra constraints (Edmonds' blossom inequalities), which is the kind of work integer programming does in Chapter 4.

Key takeaways

TU: every square submatrix has determinant in ${-1, 0, 1}$.
With a TU constraint matrix and integer right-hand side, every vertex is integer, so LP solves the integer problem.
Digraph incidence matrices and bipartite graph incidence matrices are TU. That covers every problem in this chapter.
Odd cycles break TU. Non-bipartite matching, TSP and most combinatorial problems need integer programming.

Practice

(Beginner) Is $\begin{bmatrix} 1 & 1 & 0 \ 0 & 1 & 1 \ 1 & 0 & 1 \end{bmatrix}$ TU? Compute its determinant and relate the answer to the graph it represents.

(Intermediate) Show that a matrix with entries in ${0, 1}$ where the ones of each column are consecutive (an interval matrix) is TU, by checking several random examples with is_totally_unimodular. Which scheduling problems produce such matrices?

(Interview) "Why can you solve a max-flow problem with an LP solver and get integer flows, but you cannot do the same for a TSP?" Hint: compare the constraint matrices. One is a digraph incidence matrix, the other needs subtour elimination constraints that break total unimodularity.

3.2 Shortest Path Problems

This subchapter covers the problem of finding cheapest routes: from one node to all others, with or without negative costs, and between all pairs of nodes. It ends with negative cycles, where "shortest" stops making sense.

The length of a path is the sum of its arc costs $c_{ij}$ (distance, time, money, or even negative values such as rebates). The single-source shortest path problem asks for the shortest path from a source $s$ to every node. As an LP, send one unit of flow from $s$ to $t$ at minimum cost:

$$ \min \sum_{(i,j)} c_{ij} x_{ij} \quad \text{s.t.} \quad Nx = e_s - e_t, \quad x \ge 0, $$

where $e_s - e_t$ is $+1$ at $s$, $-1$ at $t$ and 0 elsewhere. The incidence matrix is TU, so an optimal vertex is a 0/1 vector that marks a path. The dual LP has one variable $d(i)$ per node:

$$ \max\ d(t) - d(s) \quad \text{s.t.} \quad d(j) - d(i) \le c_{ij} \ \ \text{for every arc } (i, j). $$

Picture each arc as a string of length $c_{ij}$ tied between its nodes, and pull $s$ and $t$ apart as far as the strings allow. The strings that become taut form the shortest path, and the node positions are the distance labels $d(i)$. Every algorithm below computes labels that satisfy Bellman's equations

$$ d(s) = 0, \qquad d(j) = \min_{(i, j) \in E} \big( d(i) + c_{ij} \big) \quad \text{for } j \ne s, $$

using the same basic step, relaxation: if $d(i) + c_{ij} < d(j)$, set $d(j) \leftarrow d(i) + c_{ij}$ and remember $i$ as the predecessor of $j$.

Algorithm	Computes	Negative arcs	Running time	Strategy
Dijkstra	one source to all nodes	not allowed	$O((n + m) \log n)$	label setting (greedy)
Bellman-Ford	one source to all nodes	allowed, detects negative cycles	$O(nm)$	label correcting
DAG shortest path	one source, acyclic graphs only	allowed	$O(n + m)$	topological order
Floyd-Warshall	all pairs	allowed, detects negative cycles	$O(n^3)$	dynamic programming

Dijkstra's Algorithm

Why it matters. Dijkstra's algorithm (1959) is the shortest path method used in practice: in GPS navigation, in network routing protocols such as OSPF, and as a subroutine inside min-cost flow algorithms. When all costs are nonnegative, it is the fastest general method.

Theory. Imagine water poured into the network at $s$, flowing along every arc at the same speed. Nodes get wet in order of their distance from $s$, and the first time water reaches a node it has arrived by a shortest path. Dijkstra simulates this by always settling the closest unsettled node.

Set $d(s) = 0$ and $d(v) = \infty$ for every other node. All nodes start unsettled.
Pick the unsettled node $v$ with the smallest label and mark it settled (permanent).
Relax every arc $(v, w)$ leaving $v$.
Repeat until every reachable node is settled.

Why the greedy choice is safe. When $v$ has the smallest tentative label, any other route to $v$ must leave the settled set through some unsettled node $w$ with $d(w) \ge d(v)$, and then add arcs of nonnegative cost. It cannot be shorter. With a negative arc this argument fails: a later detour can make an already settled node cheaper.

Priority queue	Running time	Use when
Plain array (scan for the minimum)	$O(n^2)$	dense graphs, $m \approx n^2$
Binary heap	$O((n + m) \log n)$	sparse graphs (the usual choice)
Fibonacci heap	$O(m + n \log n)$	theory, rarely faster in practice

Example: trace on the delivery network. Each row shows the labels after settling one node (settled labels in bold).

Settled node	$a$	$b$	$c$	$d$	$t$
$s$ (0)	2	4	$\infty$	$\infty$	$\infty$
$a$ (2)	2	3	9	$\infty$	$\infty$
$b$ (3)	2	3	9	6	$\infty$
$d$ (6)	2	3	8	6	12
$c$ (8)	2	3	8	6	9
$t$ (9)	2	3	8	6	9

import heapq

ARCS = [("s", "a", 2), ("s", "b", 4), ("a", "b", 1), ("a", "c", 7),
        ("b", "d", 3), ("c", "t", 1), ("d", "c", 2), ("d", "t", 6)]

def dijkstra(arcs, source):
    succ = {}
    for i, j, c in arcs:
        succ.setdefault(i, []).append((j, c))
        succ.setdefault(j, [])
    dist = {v: float("inf") for v in succ}
    pred = {v: None for v in succ}
    dist[source] = 0
    heap = [(0, source)]                  # (tentative distance, node)
    done = set()                          # settled (permanently labeled) nodes
    while heap:
        d, v = heapq.heappop(heap)
        if v in done:                     # stale heap entry, skip it
            continue
        done.add(v)
        print(f"settle {v} at {d}")
        for w, c in succ[v]:
            if d + c < dist[w]:           # relaxation
                dist[w], pred[w] = d + c, v
                heapq.heappush(heap, (dist[w], w))
    return dist, pred

def path_to(pred, v):
    path = []
    while v is not None:
        path.append(v)
        v = pred[v]
    return path[::-1]

dist, pred = dijkstra(ARCS, "s")
print("distances:", dist)
print("shortest s->t path:", path_to(pred, "t"), "cost", dist["t"])
# settle s at 0
# settle a at 2
# settle b at 3
# settle d at 6
# settle c at 8
# settle t at 9
# distances: {'s': 0, 'a': 2, 'b': 3, 'c': 8, 'd': 6, 't': 9}
# shortest s->t path: ['s', 'a', 'b', 'd', 'c', 't'] cost 9

# Negative arc: x is settled at 2 before the cheaper detour via y is seen
bad = [("s", "x", 2), ("x", "z", 1), ("s", "y", 4), ("y", "x", -3)]
print("with a negative arc:", dijkstra(bad, "s")[0], " (true distance to z is 2)")
# settle s at 0
# settle x at 2
# settle z at 3                           <- too early: z keeps the wrong label
# settle y at 4
# with a negative arc: {'s': 0, 'x': 1, 'z': 3, 'y': 4}  (true distance to z is 2)

The pred dictionary is a rooted tree (Section 3.1): the shortest path tree, which holds the best route from $s$ to every node at once.

Key takeaways

Dijkstra settles nodes in order of distance from the source, and a settled label is final.
It requires $c_{ij} \ge 0$. A single negative arc can give wrong answers silently.
With a binary heap it runs in $O((n + m) \log n)$, and it can stop as soon as the target is settled.
The predecessors form a shortest path tree from the source to every reachable node.

Practice

(Beginner) Run Dijkstra by hand from $a$ in the delivery network. Which nodes are unreachable, and what labels do they keep?

(Intermediate) Adding the same constant $K$ to every arc cost removes negative arcs. Show with the bad network that this does not preserve shortest paths, and explain why.

(Interview) "A navigation app needs the fastest route between two points on a continental road network with 50 million nodes. Plain Dijkstra is too slow. What would you do?" Hint: think about searching from both ends (bidirectional search), guiding the search towards the target (A* with a lower bound such as straight-line distance), and precomputation (contraction hierarchies).

Bellman-Ford Algorithm

Why it matters. Negative costs appear naturally: rebates, backhaul payments, profits modeled as negative costs, and the reduced costs inside min-cost flow algorithms. Bellman-Ford handles them, and it also detects negative cycles. Its distributed version is the "distance vector" routing used in early Internet protocols (RIP).

Theory. Bellman-Ford is a label-correcting method: labels are never final until the end, and any label may be lowered later.

Set $d(s) = 0$ and $d(v) = \infty$ for every other node.
Repeat $n - 1$ rounds: relax every arc once.
Stop early if a round changes nothing.

Why $n - 1$ rounds suffice. By induction, after round $k$ each label is at most the length of the shortest path that uses at most $k$ arcs. Without negative cycles, some shortest path is simple, and a simple path has at most $n - 1$ arcs. So after $n - 1$ rounds every label is exact. This is dynamic programming (Chapter 6) with the number of arcs as the stage.

Each round costs $O(m)$, so the total is $O(nm)$. That is slower than Dijkstra, but it works with any arc signs. A common speed-up processes only nodes whose label changed in the previous round (a FIFO queue of "dirty" nodes, sometimes called SPFA). The worst case is the same, but it is often much faster in practice.

Example: a backhaul payment. In the delivery network, suppose a partner pays 4 per unit for the leg $d \to c$, so $c_{dc} = -4$. The arcs are listed from $t$ backwards, the worst possible order, so that every round extends the known paths by only one arc.

def bellman_ford(nodes, arcs, source):
    dist = {v: float("inf") for v in nodes}
    pred = {v: None for v in nodes}
    dist[source] = 0
    for rnd in range(1, len(nodes)):      # at most n - 1 rounds
        changed = False
        for i, j, c in arcs:              # relax every arc once per round
            if dist[i] + c < dist[j]:
                dist[j], pred[j] = dist[i] + c, i
                changed = True
        print(f"round {rnd}:", {v: dist[v] for v in nodes})
        if not changed:                   # labels stable -> done early
            break
    return dist, pred

NODES = ["s", "a", "b", "c", "d", "t"]
# Delivery network where the d->c leg earns a backhaul payment: cost -4.
# Arcs listed in an unlucky order (from t backwards) to show several rounds.
ARCS = [("d", "t", 6), ("c", "t", 1), ("d", "c", -4), ("b", "d", 3),
        ("a", "c", 7), ("a", "b", 1), ("s", "b", 4), ("s", "a", 2)]

dist, pred = bellman_ford(NODES, ARCS, "s")
path, v = [], "t"
while v is not None:
    path.append(v)
    v = pred[v]
print("shortest s->t path:", path[::-1], "cost", dist["t"])
# round 1: {'s': 0, 'a': 2, 'b': 4, 'c': inf, 'd': inf, 't': inf}
# round 2: {'s': 0, 'a': 2, 'b': 3, 'c': 9, 'd': 7, 't': inf}
# round 3: {'s': 0, 'a': 2, 'b': 3, 'c': 3, 'd': 6, 't': 10}
# round 4: {'s': 0, 'a': 2, 'b': 3, 'c': 2, 'd': 6, 't': 4}
# round 5: {'s': 0, 'a': 2, 'b': 3, 'c': 2, 'd': 6, 't': 3}
# shortest s->t path: ['s', 'a', 'b', 'd', 'c', 't'] cost 3

The optimal path has 5 arcs, and with this arc order each round adds one arc to the known paths, so all $n - 1 = 5$ rounds are needed. If the arcs were listed in topological order, a single round would be enough (see the DAG topic below).

Key takeaways

Bellman-Ford relaxes all arcs in rounds. After round $k$, labels are correct for paths with at most $k$ arcs.
It handles negative arc costs and runs in $O(nm)$. Stop early when a round changes nothing.
The arc order changes the number of rounds needed, not the result.
If labels still change in round $n$, the graph has a negative cycle (see the last topic of this subchapter).

Practice

(Beginner) Re-run the example with the arcs in the order $s \to a, s \to b, a \to b, a \to c, b \to d, d \to c, d \to t, c \to t$. How many rounds are needed now?

(Intermediate) Implement the queue-based variant: keep a FIFO queue of nodes whose label decreased and relax only their outgoing arcs. Count relaxations on both arc orders and compare with the round-based version.

(Interview) "Why did distance-vector routing protocols suffer from the 'count to infinity' problem when a link failed?" Hint: each router only knows its neighbors' labels, and Bellman-Ford labels decrease quickly but increase only one step per exchange when a path disappears.

Floyd-Warshall Algorithm

Why it matters. Many applications need the distance between every pair of nodes: a distance matrix for vehicle routing (Chapter 4), the metric closure used by Steiner tree heuristics (Section 3.3), or a table of travel times between all depots. Floyd-Warshall computes all of them with three nested loops.

Theory. Number the nodes $1, \dots, n$ and let $D^k(i, j)$ be the length of the shortest $i$-$j$ path whose intermediate nodes all come from ${1, \dots, k}$. A shortest such path either avoids node $k$ or passes through it exactly once (no negative cycles), so

$$ D^k(i, j) = \min\big( D^{k-1}(i, j),\ D^{k-1}(i, k) + D^{k-1}(k, j) \big), $$

with $D^0(i, j) = c_{ij}$ (or $\infty$ if there is no arc, and $0$ on the diagonal). After $k = n$, $D^n$ holds all shortest distances. The update can be done in place in one $n \times n$ matrix.

Time $O(n^3)$, memory $O(n^2)$, and very simple code with good cache behavior.
Negative arcs are fine. A negative entry on the diagonal ($D(i, i) < 0$) signals a negative cycle through $i$.
To recover paths, store $\text{next}(i, j)$, the first hop on the best $i$-$j$ path.

All-pairs method	Time	Best when
Floyd-Warshall	$O(n^3)$	dense graphs, $n$ up to a few thousand
$n$ runs of Dijkstra	$O(n (n + m) \log n)$	sparse graphs with $c \ge 0$
Johnson (Bellman-Ford reweighting + $n$ Dijkstras)	$O(nm + n (n + m) \log n)$	sparse graphs with negative arcs

Example: four warehouses. Costs between warehouses $0..3$, with one arc ($2 \to 1$) earning a rebate of 2.

INF = float("inf")

def floyd_warshall(n, arcs):
    D = [[0 if i == j else INF for j in range(n)] for i in range(n)]
    nxt = [[None] * n for _ in range(n)]  # nxt[i][j] = first hop on the i->j path
    for i, j, c in arcs:
        D[i][j], nxt[i][j] = c, j
    for k in range(n):                    # allow k as an intermediate node
        for i in range(n):
            for j in range(n):
                if D[i][k] + D[k][j] < D[i][j]:
                    D[i][j] = D[i][k] + D[k][j]
                    nxt[i][j] = nxt[i][k] # go towards k first
    if any(D[v][v] < 0 for v in range(n)):
        raise ValueError("negative cycle")
    return D, nxt

def path(nxt, i, j):
    if nxt[i][j] is None:
        return None
    p = [i]
    while i != j:
        i = nxt[i][j]
        p.append(i)
    return p

ARCS = [(0, 1, 3), (0, 2, 8), (1, 2, 4), (1, 3, 7), (2, 3, 1), (3, 0, 2), (2, 1, -2)]
D, nxt = floyd_warshall(4, ARCS)
for i, row in enumerate(D):
    print(f"from {i}:", row)
print("path 0->3:", path(nxt, 0, 3), " path 3->1:", path(nxt, 3, 1))
# from 0: [0, 3, 7, 8]
# from 1: [7, 0, 4, 5]
# from 2: [3, -2, 0, 1]
# from 3: [2, 5, 9, 0]
# path 0->3: [0, 1, 2, 3]  path 3->1: [3, 0, 1]

The direct arc $0 \to 2$ costs 8, but the matrix reports 7 via node 1: intermediate nodes are exactly what the $k$ loop adds.

Key takeaways

Floyd-Warshall is dynamic programming over the set of allowed intermediate nodes. The $k$ loop must be the outermost loop.
It computes all-pairs distances in $O(n^3)$ time and $O(n^2)$ memory, with negative arcs allowed.
A negative diagonal entry reveals a negative cycle.
For large sparse graphs, repeated Dijkstra (with Johnson's reweighting if needed) is faster.

Practice

(Beginner) Write $D^0$ and $D^1$ (after allowing node 0 as an intermediate) for the example by hand.

(Intermediate) Change the rebate on $2 \to 1$ to $-5$. What happens, and which diagonal entry turns negative first?

(Interview) "Why must the loop over the intermediate node $k$ be the outermost loop? What goes wrong if you put it inside?" Hint: $D^k$ must be completely computed from $D^{k-1}$, and with $k$ inside you would combine distances that are not yet final.

Shortest Paths in DAGs

Why it matters. Many OR networks have no directed cycles: project schedules (activities point forward in time), multi-stage decisions (Chapter 6) and time-expanded networks. On a DAG, shortest and longest paths take a single pass in linear time, even with negative costs. The longest path in a project network is the critical path of CPM (Chapter 10).

Theory. A digraph is acyclic if and only if it has a topological order: a numbering of the nodes in which every arc goes from a lower number to a higher one. Kahn's algorithm finds one in $O(n + m)$ by repeatedly removing a node with no incoming arcs. If it gets stuck before all nodes are removed, the graph has a cycle.

Once the nodes are in topological order, process them one by one and relax their outgoing arcs. When node $v$ is processed, all arcs into $v$ come from nodes processed earlier, so $d(v)$ is already final. Each arc is relaxed exactly once: $O(n + m)$ in total.

Negative costs are fine: there is no settled-too-early problem because the order is fixed by the structure, not by the labels.
Longest paths: negate all costs, or replace $\min$ with $\max$. In general graphs the longest simple path is NP-hard (it contains the Hamiltonian path problem), but in a DAG it is easy.
Dynamic programming is shortest path on a DAG: states are nodes, decisions are arcs, and the stage order is a topological order.

Example: shortest and longest routes in the delivery network. In a project-scheduling reading of the same graph, the longest path gives the earliest completion time if each arc is an activity with that duration.

from collections import deque

ARCS = [("s", "a", 2), ("s", "b", 4), ("a", "b", 1), ("a", "c", 7),
        ("b", "d", 3), ("c", "t", 1), ("d", "c", 2), ("d", "t", 6)]

def topological_order(arcs):
    """Kahn's algorithm: repeatedly remove a node with no incoming arcs."""
    succ, indeg = {}, {}
    for i, j, _ in arcs:
        succ.setdefault(i, []).append(j)
        succ.setdefault(j, [])
        indeg[j] = indeg.get(j, 0) + 1
        indeg.setdefault(i, 0)
    queue = deque(v for v in succ if indeg[v] == 0)
    order = []
    while queue:
        v = queue.popleft()
        order.append(v)
        for w in succ[v]:
            indeg[w] -= 1
            if indeg[w] == 0:
                queue.append(w)
    if len(order) < len(succ):
        raise ValueError("graph has a cycle")
    return order

def dag_paths(arcs, source, longest=False):
    order = topological_order(arcs)
    sign = -1 if longest else 1           # longest path = shortest path on -c
    dist = {v: float("inf") for v in order}
    dist[source] = 0
    out = {}
    for i, j, c in arcs:
        out.setdefault(i, []).append((j, sign * c))
    for v in order:                       # each arc is relaxed exactly once
        for w, c in out.get(v, []):
            dist[w] = min(dist[w], dist[v] + c)
    return {v: sign * d for v, d in dist.items()}

print("topological order:", topological_order(ARCS))
print("shortest from s:  ", dag_paths(ARCS, "s"))
print("longest from s:   ", dag_paths(ARCS, "s", longest=True))
# topological order: ['s', 'a', 'b', 'd', 'c', 't']
# shortest from s:   {'s': 0, 'a': 2, 'b': 3, 'd': 6, 'c': 8, 't': 9}
# longest from s:    {'s': 0, 'a': 2, 'b': 4, 'd': 7, 'c': 9, 't': 13}   <- s-b-d-t

Key takeaways

A digraph is a DAG exactly when it has a topological order, which Kahn's algorithm finds in $O(n + m)$.
Relaxing arcs in topological order solves shortest paths in one pass, with negative costs allowed.
Longest paths in DAGs are equally easy. This is the core of critical path scheduling.
Every deterministic dynamic program can be seen as a shortest or longest path on a DAG of states.

Practice

(Beginner) Find a second topological order of the delivery network, and check that the DAG algorithm gives the same distances.

(Intermediate) Modify dag_paths to also count the number of distinct shortest paths to each node. Hint: when a label ties, add counts. When it improves, replace the count.

(Interview) "Longest path is NP-hard, yet project managers compute critical paths on networks with thousands of activities every day. Explain." Hint: time only moves forward, so project networks are acyclic.

Negative Cycle Detection

Why it matters. If a cycle has negative total cost, going around it again and again lowers the cost without limit, so "shortest walk" is $-\infty$. (Shortest simple path becomes NP-hard.) Usually this signals a modeling error, but sometimes the cycle is the answer: a currency arbitrage, an inconsistent set of scheduling constraints, or (Section 3.5) a way to improve a min-cost flow.

Theory.

Bellman-Ford test. Without negative cycles, labels stabilize within $n - 1$ rounds. If some label still decreases in round $n$, a negative cycle is reachable from the source.
Finding all of them. To detect a negative cycle anywhere (not only those reachable from $s$), add a virtual source with a 0-cost arc to every node. In code, start every label at 0.
Extracting the cycle. Take the node $v$ updated in round $n$ and follow predecessor links $n$ times. After $n$ steps you must be on the cycle, because there are only $n$ nodes. Then follow the predecessors until you return to the same node.
Floyd-Warshall test. $D(i, i) < 0$ for some $i$.

Application: arbitrage. If exchanging currency $i$ into $j$ gives rate $r_{ij}$, a round trip multiplies money by $r_{i_1 i_2} r_{i_2 i_3} \cdots r_{i_k i_1}$. Profit means the product is greater than 1. Take logarithms with costs $c_{ij} = -\log r_{ij}$:

$$ \prod r_{ij} > 1 \iff \sum \log r_{ij} > 0 \iff \sum c_{ij} < 0. $$

So an arbitrage opportunity is exactly a negative cycle.

Application: difference constraints. A system of constraints $x_j - x_i \le c_{ij}$ (for example "task $j$ starts at most 5 hours after task $i$") is feasible if and only if the graph with arc $(i, j)$ of cost $c_{ij}$ has no negative cycle. If it is feasible, the shortest path labels from a virtual source are a solution.

import math

def find_negative_cycle(nodes, arcs):
    """Bellman-Ford from a virtual source linked to every node with cost 0."""
    dist = {v: 0.0 for v in nodes}        # the virtual source gives every node label 0
    pred = {v: None for v in nodes}
    for _ in range(len(nodes)):           # n rounds: the n-th must change nothing
        last = None
        for i, j, c in arcs:
            if dist[i] + c < dist[j] - 1e-12:
                dist[j], pred[j] = dist[i] + c, i
                last = j
        if last is None:
            return None                   # converged: no negative cycle
    v = last
    for _ in range(len(nodes)):           # walk back n steps to land on the cycle
        v = pred[v]
    cycle, u = [v], pred[v]
    while u != v:
        cycle.append(u)
        u = pred[u]
    cycle.append(v)
    return cycle[::-1]                    # pred links point backwards

RATES = {("USD", "EUR"): 0.92, ("EUR", "USD"): 1.08, ("EUR", "GBP"): 0.86,
         ("GBP", "EUR"): 1.16, ("GBP", "USD"): 1.27, ("USD", "GBP"): 0.78,
         ("USD", "JPY"): 150.0, ("JPY", "USD"): 0.0066}
nodes = ["USD", "EUR", "GBP", "JPY"]
arcs = [(i, j, -math.log(r)) for (i, j), r in RATES.items()]  # product -> sum

cycle = find_negative_cycle(nodes, arcs)
gain = math.prod(RATES[(cycle[k], cycle[k + 1])] for k in range(len(cycle) - 1))
print("arbitrage cycle:", " -> ".join(cycle), f"  1 unit becomes {gain:.4f}")

RATES[("GBP", "USD")] = 1.25              # market corrects the GBP/USD quote
arcs = [(i, j, -math.log(r)) for (i, j), r in RATES.items()]
print("after correction:", find_negative_cycle(nodes, arcs))
# arbitrage cycle: USD -> EUR -> GBP -> USD   1 unit becomes 1.0048   <- 0.48% profit
# after correction: None

The small tolerance 1e-12 avoids false alarms from floating point rounding in the logarithms. In practice, transaction fees would be added to the costs as well.

Key takeaways

A negative cycle makes shortest walks unbounded. Bellman-Ford detects it when labels still change in round $n$.
Start all labels at 0 (virtual source) to search the whole graph, and follow predecessors $n$ steps to land on the cycle.
Multiplicative problems (exchange rates, reliabilities) become additive with $c_{ij} = -\log r_{ij}$.
Difference-constraint systems are feasible exactly when their graph has no negative cycle. Min-cost flow optimality uses the same test (Section 3.5).

Practice

(Beginner) Verify by hand that USD $\to$ EUR $\to$ GBP $\to$ USD multiplies money by $0.92 \times 0.86 \times 1.27$. What is the largest GBP/USD rate that removes all arbitrage?

(Intermediate) Tasks A, B, C must satisfy: B starts at least 3 hours after A, C starts at least 2 hours after B, and C starts at most 4 hours after A. Write these as difference constraints, build the graph, and decide feasibility. Then change "4 hours" to "6 hours".

(Interview) "How would you detect arbitrage among 200 currencies that update every few milliseconds?" Hint: think about the $O(nm)$ cost per run, incremental updates when only a few rates change, and the effect of fees and floating point tolerance.

3.3 Minimum Spanning Trees

This subchapter shows how to connect all nodes of a network as cheaply as possible, why simple greedy rules are optimal for this problem, and how the problem becomes hard when only some of the nodes need to be connected.

Given a connected undirected graph with edge weights $w_e$, a minimum spanning tree (MST) is a spanning tree with the smallest total weight. Typical uses are laying cable, pipelines or roads that must reach every site, clustering data points (Kruskal's algorithm stopped early), and building blocks for approximation algorithms such as the MST-based TSP heuristic (Chapter 4).

The running example is a campus with six buildings $A$-$F$. The edge weights are trenching costs for a fiber network (in thousands):

     A -----4----- B -----10----- D -----11----- F
     |           /                |            /
     2         5                  4          6
     |       /                    |        /
     C -----------3-------------- E ------

Cut and Cycle Properties

Why it matters. Greedy algorithms usually fail in optimization, but for spanning trees they are provably optimal. Two simple properties explain why. They also let you check by hand whether a single edge belongs to an MST, without computing the whole tree.

Theory.

Cut property. For any cut $(S, V \setminus S)$, a lightest edge crossing the cut belongs to some MST. If it is the unique lightest edge, it belongs to every MST.
Cycle property. For any cycle, a heaviest edge on the cycle belongs to no MST, provided it is strictly heavier than all other edges of the cycle.

Proof of the cut property (exchange argument). Let $e$ be a lightest edge crossing the cut, and let $T$ be an MST that does not contain $e$. Adding $e$ to $T$ creates a cycle (Section 3.1). That cycle crosses the cut once on $e$, so it must cross back on some other edge $f$. Since $e$ is lightest across the cut, $w_e \le w_f$. Then $T + e - f$ is a spanning tree with weight at most $w(T)$, so it is also an MST, and it contains $e$. The cycle property follows by a symmetric argument.

  Exchange argument

     S                 V \ S
   +-----+   e (light)   +-----+
   |  o--+---------------+--o  |
   |  |  |               |  |  |       T + e contains a cycle that
   |  o--+---------------+--o  |       crosses the cut twice: on e and f.
   +-----+   f (heavy)   +-----+       Swap f for e: the tree gets no heavier.

These two rules give a generic greedy algorithm: repeatedly take a lightest edge across a cut (Prim, Borůvka) or reject a heaviest edge on a cycle. Kruskal does both: it accepts an edge when it joins two components (lightest across their cut), and it rejects an edge when it closes a cycle (it is the heaviest edge of that cycle, since edges come in increasing order).

Two further consequences:

If all edge weights are distinct, the MST is unique.
The MST path between two nodes minimizes the largest edge weight on the path (the minimax or bottleneck path). This is useful for problems like "the route with the best worst link".

Example: checking the properties by brute force.

import itertools

EDGES = [("A", "B", 4), ("A", "C", 2), ("B", "C", 5), ("B", "D", 10),
         ("C", "E", 3), ("D", "E", 4), ("D", "F", 11), ("E", "F", 6)]
NODES = "ABCDEF"

def is_spanning_tree(tree):
    comp = {v: v for v in NODES}          # naive component labels
    for u, v, _ in tree:
        if comp[u] == comp[v]:
            return False                  # edge closes a cycle
        old, new = comp[u], comp[v]
        comp = {x: new if c == old else c for x, c in comp.items()}
    return True

trees = [T for T in itertools.combinations(EDGES, len(NODES) - 1) if is_spanning_tree(T)]
best = min(trees, key=lambda T: sum(w for *_, w in T))
print(f"{len(trees)} spanning trees, minimum weight {sum(w for *_, w in best)}:")
print("  ", [u + v for u, v, _ in best])

# Cut property: S = {A, B, C} vs {D, E, F}
S = set("ABC")
crossing = [(u + v, w) for u, v, w in EDGES if (u in S) != (v in S)]
print("edges crossing the cut:", crossing, "-> lightest is in the MST")

# Cycle property: cycle B-D-E-C-B
cycle = [("B", "D", 10), ("D", "E", 4), ("C", "E", 3), ("B", "C", 5)]
heaviest = max(cycle, key=lambda e: e[2])
print("heaviest edge on cycle B-D-E-C:", heaviest[0] + heaviest[1],
      "in MST?", heaviest in best)
# 30 spanning trees, minimum weight 19:
#    ['AB', 'AC', 'CE', 'DE', 'EF']
# edges crossing the cut: [('BD', 10), ('CE', 3)] -> lightest is in the MST
# heaviest edge on cycle B-D-E-C: BD in MST? False

Key takeaways

Cut property: the lightest edge across any cut is safe to add.
Cycle property: the heaviest edge on any cycle is safe to discard.
Both are proved by the same exchange argument: swapping edges along the unique cycle never makes the tree heavier.
Distinct weights give a unique MST, and MST paths are minimax (bottleneck) paths.

Practice

(Beginner) Use only the cut property, with cuts that isolate a single node, to show that edges $AC$ and $EF$ belong to every MST of the campus graph.

(Intermediate) Edges $AB$ and $DE$ both weigh 4. Is the campus MST unique anyway? Justify with the cycle property.

(Interview) "The weight of one edge of a graph decreases. How do you update the MST without recomputing it from scratch?" Hint: if the edge is not in the tree, add it and remove the heaviest edge on the cycle it creates.

Kruskal's Algorithm

Why it matters. Kruskal's algorithm (1956) is the simplest MST method to implement and explain: sort the edges and keep each one that does not close a cycle. It is the natural choice for sparse graphs and edge lists, and stopping it early gives single-linkage clustering.

Theory.

Sort the edges by increasing weight.
Start with every node in its own component (a forest of $n$ single-node trees).
For each edge ${u, v}$ in order: if $u$ and $v$ are in different components, add the edge and merge the components. Otherwise skip it, because it would close a cycle.
Stop after $n - 1$ edges.

The key data structure is union-find (disjoint-set union): find(v) returns a representative of $v$'s component and union(u, v) merges two components. With path compression and union by rank, a sequence of operations costs $O(\alpha(n))$ each, where $\alpha$ (inverse Ackermann) is at most 4 for any realistic $n$. Sorting dominates, so Kruskal runs in $O(m \log m)$.

Clustering. Stopping Kruskal when $k$ components remain (equivalently, removing the $k - 1$ heaviest MST edges) gives $k$ clusters that maximize the minimum distance between clusters. This is single-linkage clustering.

class UnionFind:
    """Disjoint sets with path compression and union by rank."""
    def __init__(self, items):
        self.parent = {v: v for v in items}
        self.rank = {v: 0 for v in items}

    def find(self, v):
        if self.parent[v] != v:
            self.parent[v] = self.find(self.parent[v])  # path compression
        return self.parent[v]

    def union(self, u, v):
        ru, rv = self.find(u), self.find(v)
        if ru == rv:
            return False                  # same component: edge would close a cycle
        if self.rank[ru] < self.rank[rv]:
            ru, rv = rv, ru
        self.parent[rv] = ru              # hang the shorter tree under the taller one
        self.rank[ru] += self.rank[ru] == self.rank[rv]
        return True

def kruskal(nodes, edges):
    uf, tree = UnionFind(nodes), []
    for u, v, w in sorted(edges, key=lambda e: e[2]):  # lightest first
        if uf.union(u, v):
            tree.append((u, v, w))
            print(f"take {u}{v} ({w})")
            if len(tree) == len(nodes) - 1:
                break                     # spanning tree complete
        else:
            print(f"skip {u}{v} ({w})  closes a cycle")
    return tree

EDGES = [("A", "B", 4), ("A", "C", 2), ("B", "C", 5), ("B", "D", 10),
         ("C", "E", 3), ("D", "E", 4), ("D", "F", 11), ("E", "F", 6)]
tree = kruskal("ABCDEF", EDGES)
print("MST weight:", sum(w for *_, w in tree))
# take AC (2)
# take CE (3)
# take AB (4)
# take DE (4)
# skip BC (5)  closes a cycle            <- A, B, C already connected
# take EF (6)                            <- 5 edges: stop, BD and DF never examined
# MST weight: 19

Key takeaways

Kruskal adds edges in increasing weight order and skips those that close a cycle.
Union-find makes the cycle test almost constant time, so sorting dominates: $O(m \log m)$.
The partial forests built along the way are optimal for clustering into $k$ groups.
It needs only an edge list, so it suits sparse graphs and data that arrives as a list of distances.

Practice

(Beginner) Run Kruskal by hand on the campus graph with $BC$ changed to weight 1. Which edge leaves the tree compared with the original MST?

(Intermediate) Use the MST to split the six buildings into 2 clusters, then 3 clusters. Which edges do you remove?

(Interview) "Implement union-find. Why do you need both path compression and union by rank, and what is the complexity without them?" Hint: without them, a chain of unions can create a tree of height $n$, so find becomes $O(n)$.

Prim's Algorithm

Why it matters. Prim's algorithm grows a single tree outward from a start node, like a crystal. It looks like Dijkstra and uses the same priority-queue machinery, and with an array implementation it is the best choice for dense graphs, where sorting all $m \approx n^2$ edges would be wasteful.

Theory.

Start with a tree containing one arbitrary root node.
Among all edges with exactly one endpoint in the tree, add a lightest one (cut property with $S$ = current tree).
Repeat until the tree spans all nodes.

The only difference from Dijkstra is the key used to pick the next node: Prim uses the weight of the single connecting edge, while Dijkstra uses the total distance from the source.

Algorithm	Strategy	Data structure	Time	Best for
Kruskal	global: lightest safe edge anywhere	sorting + union-find	$O(m \log m)$	sparse graphs, edge lists
Prim (binary heap)	local: grow one tree	priority queue	$O(m \log n)$	sparse to medium graphs
Prim (array)	local: grow one tree	array of best connections	$O(n^2)$	dense graphs, complete distance matrices
Borůvka	every component picks its lightest edge at once	union-find	$O(m \log n)$	parallel and distributed computing

Example: Prim from building $D$.

import heapq

def prim(edges, root):
    adj = {}
    for u, v, w in edges:
        adj.setdefault(u, []).append((w, v))
        adj.setdefault(v, []).append((w, u))
    in_tree, tree = {root}, []
    heap = [(w, root, v) for w, v in adj[root]]  # candidate edges leaving the tree
    heapq.heapify(heap)
    while heap and len(in_tree) < len(adj):
        w, u, v = heapq.heappop(heap)     # lightest edge crossing the cut
        if v in in_tree:
            continue                      # both ends already inside: stale entry
        in_tree.add(v)
        tree.append((u, v, w))
        print(f"add {u}{v} ({w})  tree = {''.join(sorted(in_tree))}")
        for w2, x in adj[v]:
            if x not in in_tree:
                heapq.heappush(heap, (w2, v, x))
    return tree

EDGES = [("A", "B", 4), ("A", "C", 2), ("B", "C", 5), ("B", "D", 10),
         ("C", "E", 3), ("D", "E", 4), ("D", "F", 11), ("E", "F", 6)]
tree = prim(EDGES, "D")
print("MST weight:", sum(w for *_, w in tree))
# add DE (4)  tree = DE
# add EC (3)  tree = CDE
# add CA (2)  tree = ACDE
# add AB (4)  tree = ABCDE
# add EF (6)  tree = ABCDEF
# MST weight: 19                         <- same tree as Kruskal, built in another order

Key takeaways

Prim grows one tree, always adding the lightest edge that leaves it.
It is Dijkstra with a different key: edge weight instead of path length.
Use a heap for sparse graphs ($O(m \log n)$) and a plain array for dense graphs ($O(n^2)$).
Different algorithms, start nodes and tie-breaks give the same weight, and with distinct weights the same tree.

Practice

(Beginner) Run Prim by hand from $A$. List the edges in the order they are added.

(Intermediate) Implement the $O(n^2)$ array version of Prim for a distance matrix, and test it on 200 random points in the plane with Euclidean distances. Compare its running time with the heap version on the complete graph.

(Interview) "Dijkstra and Prim look almost identical. Give a graph where the shortest path tree from a node differs from the MST." Hint: a triangle with weights 1, 1 and 1.5 is enough.

Steiner Tree Problem

Why it matters. Often only some nodes must be connected (the terminals: buildings, customers, chip pins), while other nodes (junctions, hubs) may be used if they help. This is the Steiner tree problem. A small change in the question turns the easy MST into an NP-hard problem, so it is a good example of how fragile tractability is.

Theory. Given an undirected graph with edge weights and a set of terminals $T \subseteq V$, find a minimum-weight tree that contains all terminals. The other nodes it uses are called Steiner nodes.

Case	Problem	Complexity
$\lvert T \rvert = 2$	shortest path	polynomial
$T = V$	minimum spanning tree	polynomial
general $T$	Steiner tree	NP-hard
fixed number of terminals	dynamic programming (Dreyfus-Wagner)	$O(3^{\lvert T \rvert} n + \dots)$, exponential only in $\lvert T \rvert$

MST heuristic (2-approximation).

Build the metric closure on the terminals: a complete graph whose edge weights are shortest path distances (Floyd-Warshall, Section 3.2).
Compute an MST of the closure.
Replace each closure edge by its shortest path in the original graph, and remove duplicate edges.

The result costs at most $2(1 - 1/\lvert T \rvert)$ times the optimum. Why: walk around the optimal Steiner tree, visiting every edge twice. This tour has cost $2 \cdot \text{OPT}$ and visits all terminals, and shortcutting it gives a spanning tree on the closure. Better ratios exist (about 1.39 with LP-based methods), and exact solutions use integer programming (Chapter 4).

Example: four buildings and an optional hub. Four buildings sit at the corners of a square. Cables along the sides cost 5 and cables from a central hub $H$ cost 3.

   A -------5------- B
   | \             / |
   |   3         3   |
   |     \     /     |
   5        H        5        terminals: A, B, C, D
   |     /     \     |        Steiner node: H
   |   3         3   |
   | /             \ |
   D -------5------- C

import itertools

def mst(nodes, edges):
    """Kruskal; returns (weight, edges) or None if the nodes are not connected."""
    comp = {v: v for v in nodes}
    def find(v):
        while comp[v] != v:
            v = comp[v]
        return v
    tree = []
    for u, v, w in sorted(edges, key=lambda e: e[2]):
        if u in comp and v in comp and find(u) != find(v):
            comp[find(u)] = find(v)
            tree.append((u, v, w))
    if len(tree) != len(nodes) - 1:
        return None
    return sum(w for *_, w in tree), tree

def steiner_exact(terminals, steiner_nodes, edges):
    """Try every subset of Steiner nodes: exponential, fine for tiny graphs."""
    best = None
    for k in range(len(steiner_nodes) + 1):
        for extra in itertools.combinations(steiner_nodes, k):
            res = mst(set(terminals) | set(extra), edges)
            if res and (best is None or res[0] < best[0]):
                best = res
    return best

def steiner_mst_heuristic(terminals, edges):
    """2-approximation: MST of the metric closure on the terminals, then expand paths."""
    nodes = {u for u, v, _ in edges} | {v for u, v, _ in edges}
    d = {(u, v): 0 if u == v else float("inf") for u in nodes for v in nodes}
    nxt = {}
    for u, v, w in edges:
        d[u, v] = d[v, u] = w
        nxt[u, v], nxt[v, u] = v, u
    for k, i, j in itertools.product(nodes, repeat=3):   # Floyd-Warshall, k outermost
        if d[i, k] + d[k, j] < d[i, j]:
            d[i, j], nxt[i, j] = d[i, k] + d[k, j], nxt[i, k]
    closure = [(a, b, d[a, b]) for a, b in itertools.combinations(terminals, 2)]
    used = set()
    for a, b, _ in mst(terminals, closure)[1]:          # expand each closure edge
        while a != b:
            used.add(frozenset((a, nxt[a, b])))
            a = nxt[a, b]
    weight = {frozenset((u, v)): w for u, v, w in edges}
    return sum(weight[e] for e in used), sorted("".join(sorted(e)) for e in used)

EDGES = [("A", "B", 5), ("B", "C", 5), ("C", "D", 5), ("D", "A", 5),
         ("A", "H", 3), ("B", "H", 3), ("C", "H", 3), ("D", "H", 3)]
T = ["A", "B", "C", "D"]
print("MST of terminals only:", mst(set(T), EDGES)[0])
print("MST heuristic        :", steiner_mst_heuristic(T, EDGES))
w, tree = steiner_exact(T, ["H"], EDGES)
print("exact Steiner tree   :", w, sorted(u + v for u, v, _ in tree))
# MST of terminals only: 15
# MST heuristic        : (15, ['AB', 'AD', 'BC'])   <- closure distances are 5 and 6
# exact Steiner tree   : 12 ['AH', 'BH', 'CH', 'DH'] <- the hub saves 20%

The heuristic misses the hub because each pair of terminals is cheaper to connect directly (5) than through $H$ (6). Only when all four terminals share $H$ does the hub pay off. Pairwise reasoning cannot see this, which is the source of the approximation gap.

Key takeaways

Steiner tree: connect a subset of terminals, optionally through extra nodes. It generalizes both shortest path and MST.
It is NP-hard in general, but polynomial for a fixed number of terminals.
The metric-closure MST heuristic is a 2-approximation that is simple and fast.
Shared junctions (hubs) are where Steiner trees beat terminal-only spanning trees.

Practice

(Beginner) In the example, for which spoke cost $h$ does the hub stop being useful? Compare $4h$ with 15.

(Intermediate) Three cities sit at the corners of an equilateral triangle with side 1. Compare the MST length with the length of the tree that meets at the center (the Fermat point). What ratio do you get?

(Interview) "You must connect 30 offices with leased lines and may also rent space in any of 200 data centers as junctions. How would you approach this?" Hint: start with the MST heuristic for a quick bound, then consider an integer programming formulation (Chapter 4) or local search that adds and removes Steiner nodes.

3.4 Maximum Flow Problems

This subchapter asks how much can be pushed through a capacitated network from a source to a sink. It develops the augmenting path and push-relabel algorithms, proves the max-flow min-cut theorem, and applies flows to matching.

Given a digraph with capacities $u_{ij} \ge 0$, a source $s$ and a sink $t$, the maximum flow problem is the LP

$$ \max\ v \quad \text{s.t.} \quad \sum_{j} x_{ij} - \sum_{j} x_{ji} = \begin{cases} v & i = s \ -v & i = t \ 0 & \text{otherwise} \end{cases}, \qquad 0 \le x_{ij} \le u_{ij}. $$

Applications range from pipeline and telecom throughput to evacuation planning, airline crew and baseball-elimination problems, image segmentation (via minimum cuts) and every kind of matching. In the delivery network (Section 3.1), the question is: how many units per day can reach $t$? The answer, found several ways below, is 13.

Flow Networks and Residual Graphs

Why it matters. A greedy "push flow along any open path" strategy can get stuck below the optimum, because an early routing decision can block better ones. The residual graph fixes this by keeping track of how flow can be added and how existing flow can be rerouted. Every max-flow and min-cost-flow algorithm in this chapter works on it.

Theory. A flow $x$ assigns a value $x_{ij}$ to every arc such that:

Capacity: $0 \le x_{ij} \le u_{ij}$ for every arc.
Conservation: inflow = outflow at every node except $s$ and $t$.
Value: $v(x)$ = net flow out of $s$ = net flow into $t$.

The residual graph $G_x$ describes the changes still possible. For each arc $(i, j)$ with flow $x_{ij}$:

Residual arc	Residual capacity	Exists when	Meaning
forward $(i, j)$	$u_{ij} - x_{ij}$	$x_{ij} < u_{ij}$	send more along the arc
backward $(j, i)$	$x_{ij}$	$x_{ij} > 0$	cancel (undo) flow already sent

An augmenting path is an $s$-$t$ path in the residual graph. Pushing $\delta$ = (smallest residual capacity on the path) along it increases the flow value by $\delta$ and keeps the flow feasible. A backward arc in the path means "reroute flow that previously went through $(i, j)$ somewhere else", like an "undo" button for earlier decisions.

Example: 6 units on $s$-$a$-$c$-$t$.

  Flow after pushing 6 units on s-a-c-t, label = flow/capacity

                   6/6
            a ----------------> c
          ^ |                   ^ \
    6/10 /  | 0/5          0/4  |  \ 6/8
        /   |                   |   v
      s     |                   |     t
        \   |                   |   ^
    0/8  \  |                   |  / 0/7
          v v                   | /
            b ----------------> d
                   0/7

CAP = {("s", "a"): 10, ("s", "b"): 8, ("a", "b"): 5, ("a", "c"): 6,
       ("b", "d"): 7, ("c", "t"): 8, ("d", "c"): 4, ("d", "t"): 7}

def check_flow(cap, flow, s, t):
    assert all(0 <= flow[a] <= cap[a] for a in cap), "capacity violated"
    nodes = {v for arc in cap for v in arc}
    for v in nodes - {s, t}:
        inflow = sum(f for (i, j), f in flow.items() if j == v)
        outflow = sum(f for (i, j), f in flow.items() if i == v)
        assert inflow == outflow, f"conservation violated at {v}"
    return sum(f for (i, j), f in flow.items() if i == s)   # flow value

def residual(cap, flow):
    """Residual arcs with positive residual capacity."""
    r = {}
    for (i, j), u in cap.items():
        if flow[i, j] < u:
            r[i, j] = u - flow[i, j]      # forward: room to push more
        if flow[i, j] > 0:
            r[j, i] = flow[i, j]          # backward: room to cancel flow
    return r

flow = {a: 0 for a in CAP}
for arc in [("s", "a"), ("a", "c"), ("c", "t")]:   # 6 units along s-a-c-t
    flow[arc] = 6
print("flow value:", check_flow(CAP, flow, "s", "t"))
for (i, j), r in sorted(residual(CAP, flow).items()):
    kind = "forward" if (i, j) in CAP else "backward"
    print(f"  {i}->{j}  residual {r:2d}  ({kind})")
# flow value: 6
#   a->b  residual  5  (forward)
#   a->s  residual  6  (backward)
#   b->d  residual  7  (forward)
#   c->a  residual  6  (backward)
#   c->t  residual  2  (forward)
#   d->c  residual  4  (forward)
#   d->t  residual  7  (forward)
#   s->a  residual  4  (forward)
#   s->b  residual  8  (forward)
#   t->c  residual  6  (backward)

Arc $a \to c$ is saturated, so it has no forward residual arc, only the backward arc $c \to a$.

Key takeaways

A flow respects capacities and conserves flow at every node except $s$ and $t$.
The residual graph has forward arcs (spare capacity) and backward arcs (flow that can be cancelled).
An augmenting path is an $s$-$t$ path in the residual graph. Pushing its bottleneck capacity increases the flow value.
Backward arcs let algorithms undo poor early choices, which is what makes augmenting path methods optimal.

Practice

(Beginner) Starting from the flow in the example, find an augmenting path that uses only forward arcs, and push as much as possible along it. What is the new flow value?

(Intermediate) In the delivery network, push 5 units along $s$-$a$-$b$-$d$-$t$ and 4 along $s$-$a$-$c$-$t$. Draw the residual graph and find an augmenting path that must use a backward arc.

(Interview) "Why is the greedy algorithm 'repeatedly push flow along any path with spare capacity in the original graph' not optimal?" Hint: build a small network where the first path chosen blocks two better paths, and show how a backward residual arc repairs it.

Ford-Fulkerson Method

Why it matters. Ford and Fulkerson (1956) gave the first general max-flow algorithm and, with it, the max-flow min-cut theorem. Their method is the template for every augmenting path algorithm. It also proves the integrality theorem: with integer capacities there is an integer maximum flow.

Theory.

Start with $x = 0$.
While the residual graph contains an $s$-$t$ path $P$:
- let $\delta = \min_{(i,j) \in P} r_{ij}$ (bottleneck residual capacity),
- push $\delta$ along $P$: increase flow on forward arcs, decrease it on the original arcs of backward arcs.
Return $x$. It is maximum (proved in the max-flow min-cut topic).

It is a method rather than an algorithm because it does not say which augmenting path to choose, and that choice matters a lot:

Integer capacities: each augmentation adds at least 1, so there are at most $v^$ augmentations of $O(m)$ each: $O(m v^)$. This is pseudo-polynomial, since $v^*$ can be exponential in the input size.
Bad choices are slow. In the network below, alternating between $s$-$a$-$b$-$t$ and $s$-$b$-$a$-$t$ (through the backward arc) pushes 1 unit at a time: 2000 augmentations instead of 2.
Irrational capacities: with unlucky path choices, Ford-Fulkerson may never terminate and can even converge to a non-optimal value.

  Pathological instance for arbitrary path choice

             a
    1000   / |  \   1000
          /  |   \
        s    | 1   t          good paths: s-a-t, s-b-t   ->  2 augmentations
          \  |   /            bad paths:  s-a-b-t, s-b-a-t, ... -> 2000 augmentations
    1000   \ v  /   1000
             b

Integrality theorem. If all capacities are integers, Ford-Fulkerson only ever adds and subtracts integers, so it ends with an integer maximum flow. This is the algorithmic version of total unimodularity (Section 3.1), and it is why flows can model indivisible things such as trucks, people and matches.

Example: DFS path choice on the delivery network.

def ford_fulkerson(cap, s, t):
    """Augment along any residual s-t path (found by DFS) until none is left."""
    res = dict(cap)                       # residual capacities
    for i, j in cap:
        res.setdefault((j, i), 0)         # reverse arcs start with 0
    adj = {}
    for i, j in res:
        adj.setdefault(i, []).append(j)

    def dfs(v, visited):
        if v == t:
            return [t]
        visited.add(v)
        for w in adj[v]:
            if res[v, w] > 0 and w not in visited:
                path = dfs(w, visited)
                if path:
                    return [v] + path
        return None

    value = 0
    while (path := dfs(s, set())):
        arcs = list(zip(path, path[1:]))
        delta = min(res[a] for a in arcs) # bottleneck residual capacity
        for i, j in arcs:
            res[i, j] -= delta
            res[j, i] += delta            # allow this flow to be undone later
        value += delta
        backward = [f"{i}{j}" for i, j in arcs if (i, j) not in cap]
        print(f"augment {delta} along {'-'.join(path)}"
              + (f"   uses backward arcs {backward}" if backward else ""))
    flow = {a: cap[a] - res[a] for a in cap}
    return value, flow

CAP = {("s", "a"): 10, ("s", "b"): 8, ("a", "b"): 5, ("a", "c"): 6,
       ("b", "d"): 7, ("c", "t"): 8, ("d", "c"): 4, ("d", "t"): 7}
value, flow = ford_fulkerson(CAP, "s", "t")
print("max flow value:", value)
print("arc flows:", {i + j: f for (i, j), f in flow.items()})
# augment 4 along s-a-b-d-c-t
# augment 1 along s-a-b-d-t
# augment 4 along s-a-c-t
# augment 1 along s-a-c-d-t   uses backward arcs ['cd']
# augment 2 along s-b-d-t
# augment 1 along s-b-a-c-d-t   uses backward arcs ['ba', 'cd']
# max flow value: 13
# arc flows: {'sa': 10, 'sb': 3, 'ab': 4, 'ac': 6, 'bd': 7, 'ct': 8, 'dc': 2, 'dt': 5}

The first path was the long cheap-looking route through all hubs. It used up capacity on $d \to c$ that later paths had to win back through the backward arc $c \to d$. Without backward arcs, the method would have stopped below 13.

Key takeaways

Ford-Fulkerson repeatedly augments along residual $s$-$t$ paths until none remains.
With integer capacities it terminates in at most $v^*$ augmentations, and the maximum flow is integral.
The running time depends on the path choice: arbitrary choices can be pseudo-polynomial, or non-terminating with irrational data.
Backward residual arcs are essential: they reroute flow committed by earlier augmentations.

Practice

(Beginner) Check flow conservation at every hub for the final flow printed by the example.

(Intermediate) Build the pathological network in code and write a path chooser that forces the bad alternating paths. Count the augmentations, then rerun with BFS.

(Interview) "Is Ford-Fulkerson a polynomial-time algorithm?" Hint: the bound $O(m v^*)$ depends on the numeric value of the answer, not on the number of bits needed to write the input.

Edmonds-Karp Algorithm

Why it matters. Edmonds and Karp (1972) showed that one simple rule, "always augment along a shortest path (fewest arcs)", makes Ford-Fulkerson strongly polynomial, with a running time independent of the capacities. It is the easiest max-flow algorithm with a guarantee, and the analysis idea (distances that only grow) appears again in more advanced algorithms.

Theory. Find each augmenting path with BFS in the residual graph.

Monotonicity lemma. The BFS distance $\text{dist}(s, v)$ in the residual graph never decreases from one augmentation to the next.
Critical arcs. Every augmentation saturates at least one arc (its bottleneck). Before the same arc can be the bottleneck again, flow must first go back along its reverse arc, and the distance of its tail must grow by at least 2. Distances are bounded by $n$, so each arc is critical $O(n)$ times.
Bound. There are $O(nm)$ augmentations of $O(m)$ BFS work each: $O(nm^2)$ in total, regardless of capacities.

On the pathological network above, BFS picks the 2-arc paths $s$-$a$-$t$ and $s$-$b$-$t$ first, so it finishes in 2 augmentations.

Algorithm	Time	Idea
Ford-Fulkerson (any path)	$O(m v^*)$	generic augmenting paths
Edmonds-Karp	$O(n m^2)$	shortest augmenting paths (BFS)
Capacity scaling	$O(m^2 \log U)$	augment only along paths with large residual capacity first
Dinic	$O(n^2 m)$, $O(m \sqrt{n})$ on unit networks	blocking flows in BFS level graphs
Push-relabel (FIFO)	$O(n^3)$	local pushes, no paths (next topics)
Orlin / recent almost-linear methods	$O(nm)$, $m^{1+o(1)}$	theoretical state of the art

Example: BFS on the delivery network.

from collections import deque

def edmonds_karp(cap, s, t):
    """Ford-Fulkerson with BFS: always augment along a fewest-arc residual path."""
    res = dict(cap)
    for i, j in cap:
        res.setdefault((j, i), 0)
    adj = {}
    for i, j in res:
        adj.setdefault(i, []).append(j)
    value = 0
    while True:
        pred = {s: None}
        queue = deque([s])
        while queue and t not in pred:    # BFS in the residual graph
            v = queue.popleft()
            for w in adj[v]:
                if res[v, w] > 0 and w not in pred:
                    pred[w] = v
                    queue.append(w)
        if t not in pred:
            break                         # no augmenting path: flow is maximum
        path, v = [], t
        while v is not None:
            path.append(v)
            v = pred[v]
        path.reverse()
        arcs = list(zip(path, path[1:]))
        delta = min(res[a] for a in arcs)
        for i, j in arcs:
            res[i, j] -= delta
            res[j, i] += delta
        value += delta
        print(f"augment {delta:2d} along {'-'.join(path):12s} ({len(arcs)} arcs)")
    return value, res

CAP = {("s", "a"): 10, ("s", "b"): 8, ("a", "b"): 5, ("a", "c"): 6,
       ("b", "d"): 7, ("c", "t"): 8, ("d", "c"): 4, ("d", "t"): 7}
value, res = edmonds_karp(CAP, "s", "t")
print("max flow value:", value)
# augment  6 along s-a-c-t      (3 arcs)
# augment  7 along s-b-d-t      (3 arcs)
# max flow value: 13                     <- 2 augmentations instead of 6 with DFS

Key takeaways

Edmonds-Karp is Ford-Fulkerson with BFS: always use a fewest-arc augmenting path.
Residual distances never decrease, and each arc is a bottleneck at most $O(n)$ times.
Running time $O(nm^2)$, independent of the capacity values (strongly polynomial).
Faster methods (Dinic, push-relabel) exist and are preferred for large instances.

Practice

(Beginner) Run Edmonds-Karp by hand on the pathological network and confirm it needs 2 augmentations.

(Intermediate) Add a counter of BFS node visits to edmonds_karp and to the DFS version of the previous topic. Generate random networks with 50 nodes and compare the total work.

(Interview) "Why does choosing the shortest augmenting path bound the number of iterations independently of the capacities?" Hint: follow how the BFS distance of an arc's tail must change between two times the arc is saturated.

Max-Flow Min-Cut Theorem

Why it matters. This theorem is one of the most useful results in combinatorial optimization. It says that the maximum throughput of a network equals the capacity of its narrowest bottleneck, and it gives a certificate: a cut you can show to a manager to prove that no plan can deliver more. It is LP duality (Section 2.4) in network form, and it implies matching and connectivity theorems such as König's and Menger's.

Theory. An $s$-$t$ cut is a partition $(S, \bar S)$ of the nodes with $s \in S$ and $t \in \bar S$. Its capacity counts only the arcs going forward from $S$ to $\bar S$:

$$ u(S) = \sum_{i \in S,\ j \in \bar S} u_{ij}. $$

Weak duality. For any flow $x$ and any cut $S$, adding the conservation equations over the nodes of $S$ gives

$$ v(x) = \sum_{i \in S, j \in \bar S} x_{ij} - \sum_{j \in \bar S, i \in S} x_{ji} \le u(S), $$

because forward arcs carry at most their capacity and backward arcs carry at least 0.

Theorem (Ford-Fulkerson, 1956). The following are equivalent:

$x$ is a maximum flow.
The residual graph $G_x$ has no augmenting path.
$v(x) = u(S)$ for some cut $S$.

So max flow value = min cut capacity.

Proof sketch of 2 $\Rightarrow$ 3. Let $S$ be the set of nodes reachable from $s$ in $G_x$. Since $t$ is not reachable, $S$ is an $s$-$t$ cut. Every arc from $S$ to $\bar S$ must be saturated (otherwise its head would be reachable), and every arc from $\bar S$ to $S$ must carry 0 (otherwise its backward residual arc would make its tail reachable). The inequality in weak duality therefore holds with equality.

LP view. The dual of the max-flow LP has a variable $p_i$ per node and $y_{ij}$ per arc. Because the matrix is TU, it has an optimal solution that is a 0/1 vector, with $p_i = 1$ for $i \in S$ and $y_{ij} = 1$ for the cut arcs. The min-cut problem is exactly the dual LP.

Example: certificate for the delivery network.

import itertools
# Reuses edmonds_karp and CAP from the previous topic.

value, res = edmonds_karp(CAP, "s", "t")

# Min cut from the final residual graph: S = nodes still reachable from s
S, stack = {"s"}, ["s"]
while stack:
    v = stack.pop()
    for (i, j), r in res.items():
        if i == v and r > 0 and j not in S:
            S.add(j)
            stack.append(j)
cut_arcs = [(i, j) for (i, j) in CAP if i in S and j not in S]
print("S =", sorted(S), " cut arcs:", cut_arcs,
      " capacity:", sum(CAP[a] for a in cut_arcs))

# Brute force over all 2^4 cuts (s always in S, t never)
cuts = []
for k in range(5):
    for extra in itertools.combinations("abcd", k):
        S2 = {"s", *extra}
        cuts.append((sum(u for (i, j), u in CAP.items() if i in S2 and j not in S2),
                     "s" + "".join(extra)))
cuts.sort()
print("three smallest cuts:", cuts[:3], " largest:", cuts[-1])
# augment  6 along s-a-c-t      (3 arcs)
# augment  7 along s-b-d-t      (3 arcs)
# S = ['a', 'b', 's']  cut arcs: [('a', 'c'), ('b', 'd')]  capacity: 13
# three smallest cuts: [(13, 'sab'), (15, 'sabc'), (15, 'sabcd')]  largest: (33, 'scd')

        S = {s, a, b}    |    S-bar = {c, d, t}
                         |
            a ---(6)-----+----> c
          ^ |            |      ^ \
         /  |            |      |  \
        s   |            |      |    t
         \  |            |      |  ^
          v v            |      | /
            b ---(7)-----+----> d
                         |
   cut arcs a->c (6) and b->d (7): capacity 13 = max flow
   arcs inside S (s->a, s->b, a->b) or inside S-bar (d->c, c->t, d->t) do not count

The certificate tells management exactly where to invest: extra capacity on $s \to a$ or $c \to t$ is useless until $a \to c$ or $b \to d$ is upgraded.

Key takeaways

Cut capacity counts only arcs from $S$ to $\bar S$. Every flow value is at most every cut capacity.
Max flow = min cut. The reachable set in the final residual graph is a minimum cut.
The theorem is LP strong duality for a TU problem, so the min cut is a combinatorial dual certificate.
Minimum cuts identify bottlenecks, and they solve segmentation, reliability and selection problems on their own.

Practice

(Beginner) Compute the capacity of the cut $S = {s, a, b, d}$ by hand. Why does arc $d \to c$ count but $c \to t$ not?

(Intermediate) Raise $u_{ac}$ from 6 to 10. Predict the new max flow and the new minimum cut, then verify with the code.

(Interview) "How many links must an attacker cut to disconnect two data centers, and how would you compute it?" Hint: give every link capacity 1. Max flow = min cut = the number of link-disjoint paths (Menger's theorem).

Push-Relabel Algorithm

Why it matters. Augmenting path algorithms move flow along whole $s$-$t$ paths, which is expensive on large graphs. Push-relabel (Goldberg and Tarjan, 1986) works locally, one arc at a time, and is the basis of the fastest practical max-flow codes. It is also used for image segmentation and other large min-cut problems.

Theory. Augmenting path methods keep a feasible flow and work towards optimality. Push-relabel does the opposite: it keeps a "super-optimal" preflow and works towards feasibility.

A preflow satisfies the capacities, but nodes may receive more than they send. The excess $e(v)$ = inflow - outflow is allowed to be positive. A node with $e(v) > 0$ (other than $s$, $t$) is active.
Each node has a height (label) $h(v)$. Valid heights satisfy $h(s) = n$, $h(t) = 0$, and $h(i) \le h(j) + 1$ for every residual arc $(i, j)$.
Push: if $v$ is active and the residual arc $(v, w)$ is admissible ($h(v) = h(w) + 1$), send $\min(e(v), r_{vw})$ to $w$.
Relabel: if $v$ is active but has no admissible arc, set $h(v) = 1 + \min{h(w) : r_{vw} > 0}$.

Think of water flowing downhill. Excess water at a node flows to lower neighbors through open pipes. A node holding water that cannot go anywhere is raised until water can flow out. Water that cannot reach $t$ is eventually raised above $h(s) = n$ and flows back to the source.

Saturate every arc leaving $s$, and set $h(s) = n$ and all other heights to 0.
While some node is active, push from it or relabel it.
The preflow is now a flow, and it is maximum.

Why the result is optimal. The height condition means a residual path can drop at most 1 height per arc. A residual $s$-$t$ path would drop from $n$ to 0 in at most $n - 1$ arcs, which is impossible. So there is never an augmenting path, and once all excesses are gone the flow is maximum by the max-flow min-cut theorem.

Variant	Time	Note
Generic	$O(n^2 m)$	any order of operations
FIFO queue of active nodes	$O(n^3)$	simple and robust
Highest-label first	$O(n^2 \sqrt{m})$	best in practice with heuristics

The practical speed comes from two heuristics: global relabeling (recompute exact heights by a backward BFS from $t$ now and then) and the gap heuristic (if no node has height $k$, every node above $k$ can be lifted to $n$ at once).

Example: FIFO push-relabel on the delivery network.

from collections import deque

def push_relabel(cap, s, t):
    """FIFO push-relabel. Heights guide excess 'downhill' towards t."""
    nodes = sorted({v for arc in cap for v in arc})
    res = dict(cap)
    for i, j in cap:
        res.setdefault((j, i), 0)
    adj = {v: [] for v in nodes}
    for i, j in res:
        adj[i].append(j)
    height = {v: 0 for v in nodes}
    excess = {v: 0 for v in nodes}
    height[s] = len(nodes)                # source starts n levels up
    for w in adj[s]:                      # saturate every arc out of s
        delta = res[s, w]
        res[s, w] -= delta
        res[w, s] += delta
        excess[w] += delta
    active = deque(v for v in nodes if excess[v] > 0 and v not in (s, t))
    pushes = relabels = 0
    while active:
        v = active.popleft()
        while excess[v] > 0:              # discharge v completely
            for w in adj[v]:
                if res[v, w] > 0 and height[v] == height[w] + 1:  # admissible arc
                    delta = min(excess[v], res[v, w])
                    res[v, w] -= delta
                    res[w, v] += delta
                    excess[v] -= delta
                    excess[w] += delta
                    pushes += 1
                    if w not in (s, t) and excess[w] == delta:   # w just became active
                        active.append(w)
                    if excess[v] == 0:
                        break
            else:                         # no admissible arc left: lift v
                height[v] = 1 + min(height[w] for w in adj[v] if res[v, w] > 0)
                relabels += 1
    return excess[t], height, pushes, relabels

CAP = {("s", "a"): 10, ("s", "b"): 8, ("a", "b"): 5, ("a", "c"): 6,
       ("b", "d"): 7, ("c", "t"): 8, ("d", "c"): 4, ("d", "t"): 7}
value, height, pushes, relabels = push_relabel(CAP, "s", "t")
print("max flow value:", value)
print("final heights:", height)
print(f"{pushes} pushes, {relabels} relabels")
# max flow value: 13
# final heights: {'a': 7, 'b': 7, 'c': 1, 'd': 1, 's': 6, 't': 0}
# 10 pushes, 8 relabels

The final heights show the minimum cut. Nodes $a$ and $b$ rose above $h(s) = 6$ to send their surplus (18 sent from $s$, only 13 deliverable) back to the source. They are the source side of the cut, together with $s$.

Key takeaways

Push-relabel keeps a preflow (excess allowed) and valid heights, and pushes excess downhill along admissible arcs.
Valid heights guarantee there is never an augmenting path, so the algorithm ends with a maximum flow.
FIFO selection gives $O(n^3)$. With highest-label selection plus gap and global relabeling heuristics, it is among the fastest practical methods.
Excess that cannot reach $t$ returns to $s$, and the heights reveal the minimum cut.

Practice

(Beginner) After the initialization, which nodes are active and what are their excesses?

(Intermediate) Change the active-node selection to "highest height first" (use a dictionary of buckets). Compare the number of pushes and relabels on the delivery network and on a random 100-node network.

(Interview) "Compare augmenting path and push-relabel algorithms: what invariant does each maintain, and what does each work towards?" Hint: one keeps feasibility (a flow) and works towards optimality, the other keeps an optimality certificate (no augmenting path) and works towards feasibility.

Bipartite Matching via Flows

Why it matters. Assigning workers to shifts, students to projects, ads to slots or organ donors to patients are all matching problems. When each pair is either allowed or not and the goal is to match as many as possible, a max-flow algorithm solves it directly. It is the first of several "reductions to flow" that make network algorithms so widely applicable.

Theory. A bipartite graph has nodes split into two sides $L$ (workers) and $R$ (jobs), with every edge between the sides. A matching is a set of edges with no shared endpoint.

Reduction. Add a source $s$ with an arc of capacity 1 to every worker, an arc of capacity 1 from every job to a sink $t$, and orient every allowed edge from worker to job (capacity 1, or $\infty$).

   s --1--> W1 --> Pack, Ship              Pack   --1--> t
   s --1--> W2 --> Pack                    Ship   --1--> t
   s --1--> W3 --> Pack, Ship              Audit  --1--> t
   s --1--> W4 --> Ship, Audit, Repair     Repair --1--> t
                   (worker -> job arcs, capacity 1)

With integer capacities there is an integer max flow (integrality theorem). Each worker sends at most 1 unit, so the arcs with flow 1 form a matching, and max flow = maximum matching size.
An augmenting path in this network is an alternating path: it starts at an unmatched worker, uses non-matching and matching edges alternately (matching edges are backward residual arcs), and ends at an unmatched job.
König's theorem (from max-flow min-cut): in a bipartite graph, the maximum matching size equals the minimum vertex cover size (fewest nodes touching every edge). The minimum cut gives the cover.
Hall's theorem: all of $L$ can be matched if and only if every subset $A \subseteq L$ has at least $|A|$ neighbors.
Simple augmenting paths take $O(nm)$. Hopcroft-Karp augments along many shortest paths at once, in $O(m \sqrt{n})$.

Example: four workers, four tasks. W1 and W3 can do Pack or Ship, W2 can only Pack, and W4 can Ship, Audit or Repair.

def max_bipartite_matching(can_do):
    """Augmenting paths on the unit-capacity flow network s -> workers -> jobs -> t."""
    match_job = {}                        # job -> worker currently assigned

    def try_assign(worker, seen):
        for job in can_do[worker]:
            if job in seen:
                continue
            seen.add(job)
            # job is free, or its current worker can be moved to another job
            if job not in match_job or try_assign(match_job[job], seen):
                match_job[job] = worker
                return True
        return False

    for worker in can_do:
        try_assign(worker, set())
    return {w: j for j, w in match_job.items()}

def koenig_cover(can_do, matching):
    """Minimum vertex cover from a maximum matching (Koenig's theorem)."""
    match_job = {j: w for w, j in matching.items()}
    Z_workers = {w for w in can_do if w not in matching}  # unmatched workers
    Z_jobs, frontier = set(), list(Z_workers)
    while frontier:                       # alternating paths: non-matching, then matching edge
        w = frontier.pop()
        for j in can_do[w]:
            if j not in Z_jobs:
                Z_jobs.add(j)
                if j in match_job and match_job[j] not in Z_workers:
                    Z_workers.add(match_job[j])
                    frontier.append(match_job[j])
    return sorted(set(can_do) - Z_workers) + sorted(Z_jobs)

CAN_DO = {"W1": ["Pack", "Ship"], "W2": ["Pack"], "W3": ["Pack", "Ship"],
          "W4": ["Ship", "Audit", "Repair"]}
matching = max_bipartite_matching(CAN_DO)
print("maximum matching:", matching, " size", len(matching))
print("minimum vertex cover:", koenig_cover(CAN_DO, matching))
# maximum matching: {'W2': 'Pack', 'W1': 'Ship', 'W4': 'Audit'}  size 3
# minimum vertex cover: ['W4', 'Pack', 'Ship']    <- size 3 = matching size (Koenig)

The cover is a proof that 3 is optimal. Every allowed pair involves W4, Pack or Ship, so no matching can use more than three edges. In Hall's terms, W1, W2 and W3 together have only two possible tasks.

Key takeaways

Bipartite matching is max flow on a unit-capacity network $s \to L \to R \to t$.
Integrality makes the flow a matching. Augmenting paths become alternating paths.
König: maximum matching = minimum vertex cover in bipartite graphs, a direct consequence of max-flow min-cut.
Weighted versions (best total value, not just most pairs) are the assignment problem (Section 3.6).

Practice

(Beginner) Use Hall's theorem to explain why W1, W2 and W3 cannot all be matched.

(Intermediate) Extend the model so that each worker can take up to 2 tasks and task Pack needs 2 workers. Which capacities change? Solve it with edmonds_karp.

(Interview) "Design an algorithm that assigns 10,000 students to 500 projects, each project having a capacity and each student listing 5 acceptable projects." Hint: this is a flow network with capacities on the project-to-sink arcs. Discuss Hopcroft-Karp or Dinic for speed, and min-cost flow if students rank their choices.

3.5 Minimum Cost Flow Problems

This subchapter presents the most general network model of the chapter: send given supplies to given demands through a capacitated network at minimum total cost. It covers the formulation and its optimality conditions, then three algorithms that reach the optimum from different directions.

Every problem seen so far is a special case:

Problem	Supplies $b_i$	Costs $c_{ij}$	Capacities $u_{ij}$
Shortest path $s \to t$	$b_s = 1$, $b_t = -1$	arc lengths	$\infty$
Maximum flow	all 0, plus an arc $t \to s$ with cost $-1$	0 on all other arcs	as given, $\infty$ on $t \to s$
Transportation (Section 3.6)	supplies at sources, demands at sinks	unit shipping costs	$\infty$
Assignment (Section 3.6)	$+1$ per worker, $-1$ per job	assignment costs	1 (or $\infty$)
Transshipment (Section 3.6)	any, with intermediate nodes at 0	unit costs	$\infty$ or given

The running example is the delivery network (Section 3.1): ship 10 units from $s$ to $t$ at minimum total cost. The optimal cost is 101.

Min-Cost Flow Formulation

Why it matters. One model covers distribution planning, production-inventory planning over time, crew and vehicle scheduling, and many problems in disguise. Because the constraint matrix is a network matrix, the LP returns integer solutions and specialized algorithms solve instances with millions of arcs. Recognizing that a problem is a min-cost flow is often the most valuable step in solving it.

Theory.

$$ \begin{aligned} \min \quad & \sum_{(i,j) \in E} c_{ij} x_{ij} \ \text{s.t.} \quad & \sum_{j : (i,j) \in E} x_{ij} - \sum_{j : (j,i) \in E} x_{ji} = b_i && \text{for every node } i \quad (Nx = b) \ & l_{ij} \le x_{ij} \le u_{ij} && \text{for every arc } (i, j) \end{aligned} $$

Modeling tricks that keep a problem in this form:

Situation	Transformation
$\sum_i b_i > 0$ (supply exceeds demand)	add a dummy sink with demand equal to the excess, linked from each source with cost 0
Lower bound $l_{ij} > 0$	substitute $x_{ij} = l_{ij} + x'{ij}$: set $b_i \mathrel{-}= l{ij}$, $b_j \mathrel{+}= l_{ij}$, $u'{ij} = u{ij} - l_{ij}$
Node capacity (a hub handles at most $k$ units)	split the node into $i_{in} \to i_{out}$ with capacity $k$
Inventory over time	copy each node once per period and link copies with "storage" arcs (cost = holding cost)
Convex piecewise linear cost	parallel arcs, one per linear piece, with increasing costs

Optimality conditions. The dual LP assigns a potential $p_i$ to each node (a "price" of the commodity at node $i$). Define the reduced cost $\bar c_{ij} = c_{ij} + p_i - p_j$: the cost of shipping on $(i, j)$ minus the price gain from $i$ to $j$. By complementary slackness (Section 2.4), a feasible flow $x$ is optimal if and only if some potentials $p$ satisfy

Reduced cost	Required flow	Interpretation
$\bar c_{ij} > 0$	$x_{ij} = l_{ij}$	arc too expensive: use the minimum
$\bar c_{ij} < 0$	$x_{ij} = u_{ij}$	arc profitable: use it fully
$l_{ij} < x_{ij} < u_{ij}$	$\bar c_{ij} = 0$	arc used partially: exactly break-even

An equivalent statement, used by the cycle-canceling algorithm: $x$ is optimal if and only if its residual graph has no negative-cost cycle. (Residual arcs cost $c_{ij}$ forward and $-c_{ij}$ backward.)

Example: the delivery network as an LP.

import numpy as np
from scipy.optimize import linprog

NODES = ["s", "a", "b", "c", "d", "t"]
ARCS = [("s", "a", 2, 10), ("s", "b", 4, 8), ("a", "b", 1, 5), ("a", "c", 7, 6),
        ("b", "d", 3, 7), ("c", "t", 1, 8), ("d", "c", 2, 4), ("d", "t", 6, 7)]
supply = np.array([10, 0, 0, 0, 0, -10])  # b_i: ship 10 units from s to t

idx = {v: k for k, v in enumerate(NODES)}
N = np.zeros((len(NODES), len(ARCS)))
for k, (i, j, _, _) in enumerate(ARCS):
    N[idx[i], k], N[idx[j], k] = 1, -1    # outflow - inflow = b_i
cost = [c for *_, c, _ in ARCS]
bounds = [(0, u) for *_, u in ARCS]

res = linprog(cost, A_eq=N, b_eq=supply, bounds=bounds, method="highs")
print("total cost:", res.fun)
for (i, j, c, u), x in zip(ARCS, res.x):
    print(f"  {i}->{j}: flow {x:4.1f} / {u:2d}   cost {c}")
duals = res.eqlin.marginals               # one dual per conservation row
p = {v: round(float(duals[0] - duals[idx[v]])) for v in NODES}   # shift so p_s = 0
print("potentials p:", p)
print("reduced costs:", {i + j: c + p[i] - p[j] for i, j, c, _ in ARCS})
# total cost: 101.0
#   s->a: flow 10.0 / 10   cost 2
#   s->b: flow  0.0 /  8   cost 4
#   a->b: flow  5.0 /  5   cost 1
#   a->c: flow  5.0 /  6   cost 7
#   b->d: flow  5.0 /  7   cost 3
#   c->t: flow  8.0 /  8   cost 1
#   d->c: flow  3.0 /  4   cost 2
#   d->t: flow  2.0 /  7   cost 6
# potentials p: {'s': 0, 'a': 2, 'b': 4, 'c': 9, 'd': 7, 't': 13}
# reduced costs: {'sa': 0, 'sb': 0, 'ab': -1, 'ac': 0, 'bd': 0, 'ct': -3, 'dc': 0, 'dt': 0}

Check the conditions: $a \to b$ and $c \to t$ have negative reduced cost and are full. The partially used arcs ($a \to c$, $b \to d$, $d \to c$, $d \to t$) have reduced cost 0. The potential $p_t = 13$ is the marginal cost of the last unit delivered to $t$, and $\bar c_{ct} = -3$ says that one more unit of capacity on $c \to t$ would save 3.

Key takeaways

Min-cost flow: minimize $c^\top x$ subject to $Nx = b$ and $l \le x \le u$. It contains shortest path, max flow, transportation and assignment.
Lower bounds, unbalanced supplies, node capacities, time periods and convex costs all fit the model with standard transformations.
Optimality: potentials whose reduced costs are $\ge 0$ on arcs at their lower bound, $\le 0$ at their upper bound and $= 0$ in between.
Equivalently, a flow is optimal exactly when its residual graph has no negative cycle.

Practice

(Beginner) Verify that the LP flow conserves flow at every hub and costs 101.

(Intermediate) Hub $d$ can handle at most 5 units per day. Split it into $d_{in} \to d_{out}$ and re-solve. What is the new cost?

(Interview) "A factory must decide production and inventory for 12 months with monthly demand, production capacity, production cost and holding cost. Show that this is a min-cost flow problem." Hint: one node per month, a production arc from a source into each month, an inventory arc from month $k$ to $k + 1$, and the demand as the node's requirement.

Cycle-Canceling Algorithm

Why it matters. Cycle canceling turns the optimality condition directly into an algorithm. It is the simplest correct min-cost flow method and the easiest one to prove. Its idea, "improve a feasible solution along an improving cycle", is also how local search and the simplex method improve solutions.

Theory. (Klein, 1967)

Find any feasible flow, for example with a max-flow algorithm from a super-source (linked to all supplies) to a super-sink (linked from all demands).
While the residual graph contains a negative-cost cycle $W$ (found with Bellman-Ford, Section 3.2):
- let $\delta$ be the smallest residual capacity on $W$,
- push $\delta$ units around $W$. Conservation is unaffected (it is a cycle), and the cost drops by $\delta \cdot |c(W)|$.
When no negative cycle remains, the flow is optimal.

Termination. With integer data, each cancellation lowers the cost by at least 1, so there are at most (initial cost - optimal cost) iterations. That bound is pseudo-polynomial. Always canceling a minimum mean cycle (smallest average cost per arc) makes the algorithm strongly polynomial (Goldberg and Tarjan, 1989).

A cycle in the residual graph mixes forward arcs (more flow) and backward arcs (less flow), so canceling it reroutes flow: some units move from an expensive route to a cheaper one.

Example: improving a poor starting plan. Start with 7 units on $s$-$b$-$d$-$t$ (cost 13 each) and 3 units on $s$-$a$-$c$-$t$ (cost 10 each).

def negative_cycle(nodes, arcs):
    """Bellman-Ford with a virtual source; returns a cycle of (i, j) arcs or None."""
    dist = {v: 0 for v in nodes}
    pred = {v: None for v in nodes}
    for _ in range(len(nodes)):
        last = None
        for (i, j), c in arcs.items():
            if dist[i] + c < dist[j]:
                dist[j], pred[j] = dist[i] + c, i
                last = j
        if last is None:
            return None
    for _ in range(len(nodes)):           # step back onto the cycle
        last = pred[last]
    cycle, v = [], last
    while True:
        cycle.append((pred[v], v))
        v = pred[v]
        if v == last:
            return cycle[::-1]

def cycle_canceling(arcs, flow):
    nodes = {v for i, j, *_ in arcs for v in (i, j)}
    cost = {(i, j): c for i, j, c, _ in arcs}
    cap = {(i, j): u for i, j, _, u in arcs}
    while True:
        residual = {}                     # residual arc -> unit cost
        for i, j in cap:
            if flow[i, j] < cap[i, j]:
                residual[i, j] = cost[i, j]       # push more forward
            if flow[i, j] > 0:
                residual[j, i] = -cost[i, j]      # cancel flow: refund its cost
        cycle = negative_cycle(nodes, residual)
        if cycle is None:
            return flow                   # no negative cycle: optimal
        delta = min(cap[a] - flow[a] if a in cap else flow[a[::-1]] for a in cycle)
        for a in cycle:
            if a in cap:
                flow[a] += delta
            else:
                flow[a[::-1]] -= delta
        gain = sum(residual[a] for a in cycle)
        route = "-".join([cycle[0][0]] + [j for _, j in cycle])
        print(f"cancel {route:15s} cost {gain:3d} per unit x {delta}"
              f" -> total {sum(cost[a] * flow[a] for a in cap)}")

ARCS = [("s", "a", 2, 10), ("s", "b", 4, 8), ("a", "b", 1, 5), ("a", "c", 7, 6),
        ("b", "d", 3, 7), ("c", "t", 1, 8), ("d", "c", 2, 4), ("d", "t", 6, 7)]
# Feasible but expensive start: 7 units on s-b-d-t and 3 units on s-a-c-t
flow = {(i, j): 0 for i, j, *_ in ARCS}
for path, units in [("sbdt", 7), ("sact", 3)]:
    for i, j in zip(path, path[1:]):
        flow[i, j] += units
print("start cost:", sum(c * flow[i, j] for i, j, c, _ in ARCS))
flow = cycle_canceling(ARCS, flow)
print("optimal flow:", {i + j: f for (i, j), f in flow.items()})
# start cost: 121
# cancel c-t-d-c         cost  -3 per unit x 4 -> total 109
# cancel d-b-s-a-c-t-d   cost  -3 per unit x 1 -> total 106
# cancel a-b-s-a         cost  -1 per unit x 5 -> total 101
# optimal flow: {'sa': 9, 'sb': 1, 'ab': 5, 'ac': 4, 'bd': 6, 'ct': 8, 'dc': 4, 'dt': 2}

Read the first cycle as a rerouting: send 4 more units on $c \to t$ (cost +1), take 4 units off $d \to t$ (refund 6), and send them on $d \to c$ instead (cost +2). The net change is $-3$ per unit. The final flow differs from the LP solution but has the same cost (101): this problem has several optimal flows.

Key takeaways

Cycle canceling starts from any feasible flow and pushes flow around negative residual cycles until none remain.
Each canceled cycle reroutes flow to a cheaper route while keeping conservation.
It terminates with integer data. Minimum-mean cycle selection makes it strongly polynomial.
It keeps feasibility and works towards optimality, like the primal simplex method.

Practice

(Beginner) Interpret the third cycle $a$-$b$-$s$-$a$ in words: which arcs gain flow and which lose it?

(Intermediate) Start cycle canceling from another feasible flow: 6 units on $s$-$a$-$c$-$t$ and 4 units on $s$-$b$-$d$-$t$. Does it need more or fewer cancellations, and does it end with the same flow?

(Interview) "You have a working delivery plan and a new, cheaper route opens. How can you update the plan without solving from scratch?" Hint: the old plan is still feasible. Look for negative cycles that use the new arc.

Successive Shortest Path Algorithm

Why it matters. Instead of fixing a bad plan, this algorithm builds a good one unit by unit: always ship the next units along the currently cheapest route. It is the min-cost version of augmenting paths, the method behind most assignment and matching codes, and it shows that the total cost is a convex function of the amount shipped.

Theory.

Start with $x = 0$ (or $x = l$) and potentials $p = 0$. This is optimal for shipping nothing.
While demand remains: find a shortest path from a supply node to a demand node in the residual graph, push as much as possible along it (limited by residual capacities and remaining supply and demand), and update the potentials.

Invariant. After each augmentation, the residual graph has no negative cycle, so the current flow is the cheapest way to ship its amount. The algorithm keeps optimality and works towards feasibility, the mirror image of cycle canceling.

Dijkstra with potentials. Residual graphs contain negative-cost backward arcs, so plain Dijkstra seems unusable. But with potentials $p$ equal to the previous shortest path distances, the reduced costs $c_{ij} + p_i - p_j$ are all nonnegative, and Dijkstra works on them (Johnson's reweighting). Path lengths change only by the constant $p_s - p_t$, so shortest paths are preserved.

Running time. Each augmentation sends at least one unit, so with total supply $B$ there are at most $B$ Dijkstra runs: $O(B (n + m) \log n)$. This is pseudo-polynomial. Capacity scaling makes it polynomial.

Example: shipping the 10 units.

import heapq

def successive_shortest_path(arcs, s, t, demand):
    """Send `demand` units s->t, always along a cheapest residual path."""
    cap, cost = {}, {}
    for i, j, c, u in arcs:
        cap[i, j], cost[i, j] = u, c
        cap[j, i], cost[j, i] = 0, -c     # reverse arc refunds the cost
    adj = {}
    for i, j in cap:
        adj.setdefault(i, []).append(j)
    pot = {v: 0 for v in adj}             # valid start: all costs are >= 0
    sent = total = 0
    while sent < demand:
        # Dijkstra on reduced costs c_ij + pot_i - pot_j (all >= 0)
        dist, pred = {s: 0}, {s: None}
        heap, done = [(0, s)], set()
        while heap:
            d, v = heapq.heappop(heap)
            if v in done:
                continue
            done.add(v)
            for w in adj[v]:
                rc = cost[v, w] + pot[v] - pot[w]
                if cap[v, w] > 0 and d + rc < dist.get(w, float("inf")):
                    dist[w], pred[w] = d + rc, v
                    heapq.heappush(heap, (d + rc, w))
        if t not in dist:
            raise ValueError("demand exceeds the maximum flow")
        for v in dist:                    # keep reduced costs nonnegative
            pot[v] += dist[v]
        path, v = [t], t
        while pred[v] is not None:
            v = pred[v]
            path.append(v)
        path.reverse()
        legs = list(zip(path, path[1:]))
        delta = min([demand - sent] + [cap[a] for a in legs])
        unit = sum(cost[a] for a in legs) # true cost of one unit on this path
        for i, j in legs:
            cap[i, j] -= delta
            cap[j, i] += delta
        sent += delta
        total += delta * unit
        print(f"send {delta} along {'-'.join(path):14s} unit cost {unit:2d}  total {total}")
    return total

ARCS = [("s", "a", 2, 10), ("s", "b", 4, 8), ("a", "b", 1, 5), ("a", "c", 7, 6),
        ("b", "d", 3, 7), ("c", "t", 1, 8), ("d", "c", 2, 4), ("d", "t", 6, 7)]
print("minimum cost:", successive_shortest_path(ARCS, "s", "t", demand=10))
# send 4 along s-a-b-d-c-t    unit cost  9  total 36
# send 4 along s-a-c-t        unit cost 10  total 76
# send 1 along s-a-b-d-t      unit cost 12  total 88
# send 1 along s-a-c-d-t      unit cost 13  total 101   <- uses backward arc c->d
# minimum cost: 101

The unit costs 9, 10, 12, 13 never decrease. The minimum cost of shipping $v$ units is therefore a convex, piecewise linear function of $v$, whose slopes are these path costs:

Units shipped $v$	0	4	8	9	10
Minimum cost	0	36	76	88	101
Cost of each extra unit up to the next column	9	10	12	13	-

Key takeaways

Successive shortest paths ship flow along the cheapest residual path, one augmentation at a time.
Each intermediate flow is optimal for its own amount. Feasibility is reached at the end.
Potentials make all reduced costs nonnegative, so Dijkstra can be used despite negative residual arcs.
Marginal path costs never decrease, so total cost is convex in the amount shipped. The algorithm is pseudo-polynomial in the total supply.

Practice

(Beginner) Using the output, what is the minimum cost of shipping 8 units? And 9?

(Intermediate) Replace Dijkstra-with-potentials by Bellman-Ford on the raw residual costs. Check that you get the same augmentations and compare the number of relaxations.

(Interview) "Why can you run Dijkstra on a residual graph that contains negative arc costs?" Hint: show that with potentials from the previous shortest path computation, every residual arc has nonnegative reduced cost, including the new backward arcs on the last augmenting path.

Network Simplex Method

Why it matters. The network simplex method is the simplex method of Chapter 2 specialized to network matrices. Bases become spanning trees, and every matrix operation becomes a tree operation with integer arithmetic. It is often the fastest method in practice for min-cost flow problems of moderate to large size, and commercial LP solvers include a dedicated network simplex optimizer for models with network structure.

Theory. The connection to the LP simplex method (Sections 2.2-2.3):

Simplex concept	Network simplex counterpart
Basis ($m$ columns of $A$)	spanning tree of $n - 1$ arcs (Section 3.1: independent incidence columns are trees)
Nonbasic variables at a bound	nontree arcs with $x_{ij} = l_{ij}$ or $x_{ij} = u_{ij}$
Basic solution $x_B = B^{-1} b$	tree flows computed leaf by leaf from conservation
Dual values $\pi = c_B^\top B^{-1}$	potentials from $\bar c_{ij} = 0$ on tree arcs, with the root at 0
Pricing (reduced costs)	$\bar c_{ij} = c_{ij} + p_i - p_j$ for nontree arcs
Ratio test	the entering arc plus the tree path between its ends form a unique cycle. Push flow around it until some arc hits a bound
Pivot	swap the entering and leaving arcs in the tree

One iteration.

Potentials: set $p_{root} = 0$ and walk the tree, using $p_j = p_i + c_{ij}$ on every tree arc $(i, j)$.
Pricing: find a nontree arc that violates optimality: $\bar c_{ij} < 0$ at lower bound, or $\bar c_{ij} > 0$ at upper bound. If none exists, stop: optimal.
Cycle: adding the entering arc to the tree creates exactly one cycle (Section 3.1). Orient it in the direction that lowers cost.
Ratio test: push $\delta$ = the smallest residual capacity around the cycle. An arc that reaches a bound leaves the tree.
Update the tree, the flows on the cycle, and the potentials of the subtree that was cut off.

Starting tree. As in the Big-M method (Section 2.3), add an artificial root connected to every node by an arc of large cost $M$ that carries that node's supply or demand. Pivots drive the artificial flow to zero.

Degeneracy. Many pivots have $\delta = 0$ (they change the tree but not the flow). Cunningham's strongly feasible trees rule out cycling and work well in practice.

Example: network simplex on the delivery network.

def network_simplex(nodes, arcs, supply, M=100):
    """Primal network simplex with a Big-M artificial root and Dantzig pricing."""
    arcs = list(arcs)
    for v in nodes:                       # artificial arcs link every node to root "r"
        if supply[v] >= 0:
            arcs.append((v, "r", M, float("inf")))
        else:
            arcs.append(("r", v, M, float("inf")))
    m_real = len(arcs) - len(nodes)
    flow = [0] * m_real + [abs(supply[v]) for v in nodes]
    tree = set(range(m_real, len(arcs)))  # initial basis: the artificial star
    while True:
        # 1. Potentials: reduced cost c_ij + p_i - p_j = 0 on every tree arc
        p, parent, stack = {"r": 0}, {"r": None}, ["r"]
        while stack:
            v = stack.pop()
            for k in tree:
                i, j, c, _ = arcs[k]
                if i == v and j not in p:
                    p[j], parent[j] = p[i] + c, (k, i)
                    stack.append(j)
                elif j == v and i not in p:
                    p[i], parent[i] = p[j] - c, (k, j)
                    stack.append(i)
        # 2. Pricing: pick the nontree arc that violates optimality the most
        best, enter = 0, None
        for k, (i, j, c, u) in enumerate(arcs):
            if k in tree:
                continue
            rc = c + p[i] - p[j]
            viol = -rc if flow[k] == 0 else rc if flow[k] == u else 0
            if viol > best:
                best, enter = viol, k
        if enter is None:                 # optimal: no violated arc
            cost = sum(arcs[k][2] * flow[k] for k in range(len(arcs)))
            return cost, flow[:m_real]
        # 3. Cycle = entering arc + tree path back. Push i->j if at 0, j->i if at u
        i, j, c, u = arcs[enter]
        tail, head = (i, j) if flow[enter] == 0 else (j, i)
        def to_root(v):                   # tree arcs from v up to the root
            path = []
            while parent[v] is not None:
                path.append((parent[v][0], v))
                v = parent[v][1]
            return path
        up_head, up_tail = to_root(head), to_root(tail)
        common = {k for k, _ in up_head} & {k for k, _ in up_tail}
        cycle = [(enter, +1 if tail == i else -1)]
        for k, v in up_head:              # flow goes up from head to the apex
            if k in common:
                break
            cycle.append((k, +1 if arcs[k][0] == v else -1))
        for k, v in up_tail:              # ... and down from the apex to tail
            if k in common:
                break
            cycle.append((k, +1 if arcs[k][1] == v else -1))
        # 4. Ratio test: the first arc to hit a bound leaves the tree
        room = [(arcs[k][3] - flow[k]) if d > 0 else flow[k] for k, d in cycle]
        delta = min(room)
        leave = cycle[room.index(delta)][0]
        for k, d in cycle:
            flow[k] += d * delta
        if leave != enter:
            tree = (tree - {leave}) | {enter}
        name = lambda k: arcs[k][0] + arcs[k][1]
        print(f"enter {name(enter)} (reduced cost {c + p[i] - p[j]:+4d}), "
              f"push {delta}, leave {name(leave)}")

NODES = ["s", "a", "b", "c", "d", "t"]
ARCS = [("s", "a", 2, 10), ("s", "b", 4, 8), ("a", "b", 1, 5), ("a", "c", 7, 6),
        ("b", "d", 3, 7), ("c", "t", 1, 8), ("d", "c", 2, 4), ("d", "t", 6, 7)]
supply = {"s": 10, "a": 0, "b": 0, "c": 0, "d": 0, "t": -10}
cost, flow = network_simplex(NODES, ARCS, supply)
print("minimum cost:", cost)
print("flows:", {i + j: f for (i, j, *_), f in zip(ARCS, flow)})
# enter ct (reduced cost -199), push 0, leave cr    <- degenerate pivots: real arcs
# enter dc (reduced cost -197), push 0, leave dr    <- replace the artificial arcs
# enter bd (reduced cost -194), push 0, leave br    <- (those to or from root r)
# enter ab (reduced cost -193), push 0, leave ar
# enter sa (reduced cost -191), push 4, leave dc
# enter ac (reduced cost -190), push 4, leave ct
# enter dt (reduced cost -188), push 1, leave ab
# enter sb (reduced cost -187), push 1, leave rt    <- last artificial flow removed
# minimum cost: 101
# flows: {'sa': 9, 'sb': 1, 'ab': 5, 'ac': 4, 'bd': 6, 'ct': 8, 'dc': 4, 'dt': 2}

This teaching version recomputes all potentials from scratch in each iteration. Production codes store the tree with parent, depth and "thread" pointers so that finding the cycle and updating potentials touches only the affected subtree, and they use partial ("candidate list") pricing instead of scanning every arc.

Key takeaways

Network simplex is the primal simplex method in which bases are spanning trees.
Potentials come from zero reduced costs on tree arcs. The entering arc closes a unique cycle, and the ratio test picks the leaving arc on it.
All arithmetic is additions and subtractions of integers: no matrix inversion, no rounding errors, integral solutions.
Big-M artificial arcs give a starting tree. Strongly feasible trees prevent cycling under degeneracy.

Practice

(Beginner) For the final flow, list the arcs strictly between their bounds. Check that they contain no cycle, and explain why this holds for every basic (tree) solution.

(Intermediate) Change the pricing rule to "first violated arc found" (Bland-like). Count pivots and compare with Dantzig pricing.

(Interview) "Why is the network simplex method much faster per iteration than the general simplex method on the same LP?" Hint: a basis is a tree, so $B^{-1} b$ and $c_B^\top B^{-1}$ are computed by walking the tree in $O(n)$ instead of by factorizing a matrix.

3.6 Transportation and Assignment Problems

This subchapter studies the two most famous special cases of min-cost flow: shipping from sources to destinations (transportation) and pairing agents with tasks (assignment). Their special structure supports simple tableau methods, which are still the standard way these problems are taught and solved by hand, plus the Hungarian algorithm and the transshipment extension.

The running example has three warehouses supplying four customers. Costs are per unit shipped:

                  C1     C2     C3     C4   | Supply
          W1       8      6     10      9   |   30
          W2       9     12     13      7   |   40
          W3      14      9     16      5   |   30
          ----------------------------------+------
          Demand  20     25     35     20   |  100

Total supply equals total demand (100), so the problem is balanced. The optimal cost is 870.

Transportation Problem Formulation

Why it matters. The transportation problem (Hitchcock, 1941) models any situation where a homogeneous good moves directly from $m$ sources to $n$ destinations: plants to warehouses, warehouses to stores, or even production periods to demand periods. It was one of the first LP models ever solved, and its structure (each variable appears in exactly one supply row and one demand column) makes it much easier than a general LP.

Theory. Let $x_{ij}$ be the amount shipped from source $i$ (supply $s_i$) to destination $j$ (demand $d_j$) at unit cost $c_{ij}$:

$$ \min \sum_{i=1}^{m} \sum_{j=1}^{n} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_{j} x_{ij} = s_i \ \ (i = 1..m), \quad \sum_{i} x_{ij} = d_j \ \ (j = 1..n), \quad x_{ij} \ge 0. $$

Balance. A feasible solution needs $\sum_i s_i = \sum_j d_j$. If supply exceeds demand, add a dummy destination that absorbs the excess at cost 0 (unused capacity). If demand exceeds supply, add a dummy source whose cost is the penalty for unmet demand.
Rank. The $m + n$ equations sum to the same total in two ways, so one is redundant: the rank is $m + n - 1$. A basic feasible solution has at most $m + n - 1$ positive cells, and if it has fewer it is degenerate.
Basis = tree. Draw sources and destinations as nodes and basic cells as edges. Every basis is a spanning tree of this bipartite graph, so the problem is a min-cost flow on a bipartite network (Section 3.5).
Integrality. The constraint matrix is the incidence matrix of a bipartite graph, which is TU (Section 3.1). With integer supplies and demands, every basic solution is integer.
Forbidden routes get a very large cost $M$, or the variable is simply omitted.

import numpy as np
from scipy.optimize import linprog

COST = np.array([[8, 6, 10, 9],           # rows: warehouses W1..W3
                 [9, 12, 13, 7],          # columns: customers C1..C4
                 [14, 9, 16, 5]])
SUPPLY = np.array([30, 40, 30])
DEMAND = np.array([20, 25, 35, 20])
m, n = COST.shape

A_eq, b_eq = [], []
for i in range(m):                        # row i ships exactly its supply
    row = np.zeros((m, n)); row[i, :] = 1
    A_eq.append(row.ravel()); b_eq.append(SUPPLY[i])
for j in range(n):                        # column j receives exactly its demand
    col = np.zeros((m, n)); col[:, j] = 1
    A_eq.append(col.ravel()); b_eq.append(DEMAND[j])

res = linprog(COST.ravel(), A_eq=A_eq, b_eq=b_eq, bounds=(0, None), method="highs")
X = res.x.reshape(m, n)
print("optimal cost:", res.fun)
print(X.astype(int))
print("rank of the constraint matrix:", np.linalg.matrix_rank(np.array(A_eq)), "= m + n - 1")
print("positive shipments:", int((X > 0).sum()))
# optimal cost: 870.0
# [[ 0 15 15  0]
#  [20  0 20  0]
#  [ 0 10  0 20]]
# rank of the constraint matrix: 6 = m + n - 1
# positive shipments: 6                  <- a nondegenerate basic solution

Key takeaways

Transportation: sources with supplies, destinations with demands, one variable per route, equality constraints.
Balance the problem with a dummy source or destination. Block forbidden routes with a large cost.
The rank is $m + n - 1$, so a basic solution has at most $m + n - 1$ positive shipments, and they form a spanning tree.
Integer data gives integer optimal shipments because the matrix is TU.

Practice

(Beginner) Customer C4's demand rises to 30. Balance the problem with a dummy row, and explain what the dummy's shipments mean.

(Intermediate) Draw the bipartite graph of the optimal basis and check that it is a spanning tree on the 7 nodes.

(Interview) "A product is made in 4 factories with different costs per unit and shipped to 30 regional warehouses. Can the production decision be included in the same transportation model?" Hint: add the production cost of factory $i$ to every $c_{ij}$ in its row, and use the factory capacity as its supply (with a dummy destination for spare capacity).

Northwest Corner Rule

Why it matters. Tableau methods for the transportation problem need a starting basic feasible solution. The northwest corner rule produces one mechanically, with no arithmetic beyond subtraction. It ignores costs, so it is usually a poor start, but it shows clearly how a basic solution of the transportation problem is shaped.

Theory.

Start at the top-left (northwest) cell $(1, 1)$.
Ship $x_{ij} = \min(\text{remaining supply}_i, \text{remaining demand}_j)$ and subtract it from both.
If row $i$ is exhausted, move down. Otherwise (column $j$ is satisfied), move right.
Repeat until the bottom-right cell is reached.

Each step exhausts one row or column, and the last step exhausts both, so there are exactly $m + n - 1$ allocated cells. The cells form a staircase, which never contains a cycle, so it is a basis. Degeneracy: if a row and a column are exhausted at the same time (not at the end), move only one way and record the next cell with a zero allocation, so that the basis still has $m + n - 1$ cells.

Starting method	Uses costs?	Effort	Typical quality
Northwest corner	no	lowest	poor
Least-cost (cheapest cell first)	yes, greedily	low	fair
Vogel's approximation	yes, regret-based	moderate	often optimal or near-optimal

def northwest_corner(supply, demand):
    supply, demand = list(supply), list(demand)   # work on copies
    i = j = 0
    X = [[0] * len(demand) for _ in supply]
    basis = []
    while i < len(supply) and j < len(demand):
        q = min(supply[i], demand[j])     # ship as much as possible in cell (i, j)
        X[i][j] = q
        basis.append((i, j))
        supply[i] -= q
        demand[j] -= q
        if supply[i] == 0 and i < len(supply) - 1:
            i += 1                        # row exhausted: move down
        else:
            j += 1                        # column satisfied: move right
    return X, basis

COST = [[8, 6, 10, 9], [9, 12, 13, 7], [14, 9, 16, 5]]
SUPPLY, DEMAND = [30, 40, 30], [20, 25, 35, 20]
X, basis = northwest_corner(SUPPLY, DEMAND)
for row in X:
    print(row)
print("basic cells:", basis)
print("cost:", sum(COST[i][j] * X[i][j] for i, j in basis))
# [20, 10, 0, 0]
# [0, 15, 25, 0]
# [0, 0, 10, 20]
# basic cells: [(0, 0), (0, 1), (1, 1), (1, 2), (2, 2), (2, 3)]   <- the staircase
# cost: 985                               <- 13% above the optimum of 870

The plan sends 15 units on the expensive route W2-C2 (cost 12) and 10 on W3-C3 (cost 16) only because those cells happen to lie on the diagonal.

Key takeaways

The northwest corner rule fills cells from the top-left, moving down or right as rows and columns are exhausted.
It always gives a basic feasible solution with $m + n - 1$ cells (zero allocations included when degenerate).
It ignores costs, so it is fast but usually far from optimal.
Its staircase shape shows the tree structure of transportation bases.

Practice

(Beginner) Apply the northwest corner rule to supplies $(10, 20)$ and demands $(10, 5, 15)$. Is the result degenerate?

(Intermediate) Implement the least-cost method (always allocate to the cheapest remaining cell) and compare its cost with the northwest corner on the example.

(Interview) "If the northwest corner rule ignores costs, why is it still taught?" Hint: it shows the structure of a basis (exactly $m + n - 1$ cells, no cycles), and any starting basis is enough for the MODI method to reach the optimum.

Vogel's Approximation Method

Why it matters. A good starting solution saves many improvement iterations, and on small problems it is often already optimal. Vogel's approximation method (VAM) is the best of the classic tableau heuristics, and its key idea, regret, appears in many other heuristics such as regret-based insertion for vehicle routing.

Theory. For each row and column, the penalty (regret) is the difference between its two cheapest remaining costs. It measures how much you lose if that line does not get its cheapest cell.

Compute the penalty of every remaining row and column.
Pick the line with the largest penalty, the one where a wrong choice hurts most (break ties arbitrarily).
Allocate as much as possible to the cheapest remaining cell of that line.
Cross out the exhausted row or column (only one, if both are exhausted) and repeat until everything is allocated.

def vogel(cost, supply, demand):
    supply, demand = list(supply), list(demand)
    rows, cols = set(range(len(supply))), set(range(len(demand)))
    X = [[0] * len(demand) for _ in supply]
    basis = []

    def penalty(values):
        v = sorted(values)                # gap between the two cheapest options
        return v[1] - v[0] if len(v) > 1 else v[0]   # last cell: allocation is forced

    while rows and cols:
        cand = [(penalty([cost[i][j] for j in cols]), "row", i) for i in rows]
        cand += [(penalty([cost[i][j] for i in rows]), "col", j) for j in cols]
        p, kind, k = max(cand, key=lambda c: c[0])   # largest regret first
        if kind == "row":
            i, j = k, min(cols, key=lambda j: cost[k][j])
        else:
            i, j = min(rows, key=lambda i: cost[i][k]), k
        q = min(supply[i], demand[j])
        X[i][j] = q
        basis.append((i, j))
        print(f"penalty {p} on {kind} {k + 1}: ship {q} in cell (W{i + 1}, C{j + 1})")
        supply[i] -= q
        demand[j] -= q
        if supply[i] == 0 and len(rows) > 1:
            rows.discard(i)               # cross out exactly one line
        else:
            cols.discard(j)
    return X, basis

COST = [[8, 6, 10, 9], [9, 12, 13, 7], [14, 9, 16, 5]]
SUPPLY, DEMAND = [30, 40, 30], [20, 25, 35, 20]
X, basis = vogel(COST, SUPPLY, DEMAND)
for row in X:
    print(row)
print("cost:", sum(COST[i][j] * X[i][j] for i, j in basis))
# penalty 4 on row 3: ship 20 in cell (W3, C4)    <- W3: 9 - 5 = 4 is the largest regret
# penalty 5 on row 3: ship 10 in cell (W3, C2)
# penalty 6 on col 2: ship 15 in cell (W1, C2)
# penalty 4 on row 2: ship 20 in cell (W2, C1)
# penalty 13 on row 2: ship 20 in cell (W2, C3)
# penalty 10 on row 1: ship 15 in cell (W1, C3)
# [0, 15, 15, 0]
# [20, 0, 20, 0]
# [0, 10, 0, 20]
# cost: 870                               <- already optimal

Key takeaways

Vogel's method allocates first where the regret (gap between the two cheapest options) is largest.
It uses cost information globally and usually lands at or near the optimum.
It needs more work per step than the northwest corner rule, but saves improvement iterations later.
The result is still only a starting solution: check optimality with MODI.

Practice

(Beginner) Compute the initial row and column penalties of the example by hand and confirm the first allocation.

(Intermediate) Generate 100 random $5 \times 6$ balanced instances. Compare the average gap to the optimum for the northwest corner, least-cost and Vogel methods.

(Interview) "Explain the idea of 'regret' in heuristics and give an example outside transportation." Hint: in vehicle routing, insert first the customer whose best and second-best insertion positions differ the most.

MODI Method

Why it matters. The modified distribution (MODI, or u-v) method is the simplex method in transportation-tableau form. It decides whether a solution is optimal using only additions and subtractions, and improves it otherwise. It shows LP duality at work on a small table that you can solve by hand.

Theory. Give each row a dual value $u_i$ and each column a dual value $v_j$ (these are the LP duals of the supply and demand constraints).

Duals: for every basic cell, solve $u_i + v_j = c_{ij}$. There are $m + n$ unknowns and $m + n - 1$ equations, so set $u_1 = 0$.
Reduced costs: for every nonbasic cell, $d_{ij} = c_{ij} - u_i - v_j$. This is how much the total cost changes per unit if the cell is used.
Optimality test: if all $d_{ij} \ge 0$, the solution is optimal. Otherwise, the cell with the most negative $d_{ij}$ enters.
Loop (stepping stone): the entering cell and some basic cells form a unique closed loop that alternates between rows and columns, turning only at basic cells. Mark the cells $+, -, +, -, \dots$ starting with $+$ at the entering cell.
Pivot: $\theta$ = the smallest allocation among the $-$ cells. Add $\theta$ to the $+$ cells and subtract it from the $-$ cells. A $-$ cell that drops to 0 leaves the basis.

  First MODI iteration from the northwest corner solution (allocations shown)

              C1     C2      C3      C4
     W1       20     10       .       .
     W2        .   (-)15  (+)25       .
     W3        .   (+) .  (-)10      20      entering W3C2 with d = 9 - 9 - 6 = -6
                                             theta = min(15, 10) = 10, W3C3 leaves

This is network simplex (Section 3.5) on the bipartite network: the basis is a tree, $u$ and $v$ are node potentials, and the loop is the cycle the entering arc closes.

from collections import deque
# Reuses northwest_corner from the Northwest Corner Rule topic.

def modi(cost, X, basis):
    m, n = len(cost), len(cost[0])
    basis = list(basis)
    name = lambda c: f"W{c[0] + 1}C{c[1] + 1}"
    while True:
        # 1. Dual values: u_i + v_j = c_ij on every basic cell (u_1 = 0)
        u, v = [None] * m, [None] * n
        u[0] = 0
        while None in u or None in v:
            for i, j in basis:
                if u[i] is not None and v[j] is None:
                    v[j] = cost[i][j] - u[i]
                elif v[j] is not None and u[i] is None:
                    u[i] = cost[i][j] - v[j]
        # 2. Reduced costs d_ij = c_ij - u_i - v_j of the nonbasic cells
        d = {(i, j): cost[i][j] - u[i] - v[j]
             for i in range(m) for j in range(n) if (i, j) not in basis}
        enter = min(d, key=d.get)
        z = sum(cost[i][j] * X[i][j] for i, j in basis)
        print(f"cost {z}: u = {u}, v = {v}, smallest d = {d[enter]} at {name(enter)}")
        if d[enter] >= 0:
            return X, z                   # every d_ij >= 0: optimal
        # 3. Loop: path from row enter[0] to column enter[1] through basic cells
        start, goal = ("r", enter[0]), ("c", enter[1])
        pred, queue = {start: None}, deque([start])
        while queue:
            node = queue.popleft()
            kind, k = node
            for i, j in basis:
                if (kind == "r" and i == k) or (kind == "c" and j == k):
                    nxt = ("c", j) if kind == "r" else ("r", i)
                    if nxt not in pred:
                        pred[nxt] = (node, (i, j))
                        queue.append(nxt)
        loop, node = [enter], goal
        while pred[node] is not None:
            node, cell = pred[node]
            loop.append(cell)
        # 4. Signs alternate +, -, +, - around the loop. theta = min over '-' cells
        minus = loop[1::2]
        theta = min(X[i][j] for i, j in minus)
        leave = next(c for c in minus if X[c[0]][c[1]] == theta)
        for k, (i, j) in enumerate(loop):
            X[i][j] += theta if k % 2 == 0 else -theta
        basis.remove(leave)
        basis.append(enter)
        print(f"   loop {' '.join(map(name, loop))}, theta = {theta}, {name(leave)} leaves")

COST = [[8, 6, 10, 9], [9, 12, 13, 7], [14, 9, 16, 5]]
SUPPLY, DEMAND = [30, 40, 30], [20, 25, 35, 20]
X, z = modi(COST, *northwest_corner(SUPPLY, DEMAND))
for row in X:
    print(row)
# cost 985: u = [0, 6, 9], v = [8, 6, 7, -4], smallest d = -6 at W3C2
#    loop W3C2 W2C2 W2C3 W3C3, theta = 10, W3C3 leaves
# cost 925: u = [0, 6, 3], v = [8, 6, 7, 2], smallest d = -5 at W2C1
#    loop W2C1 W1C1 W1C2 W2C2, theta = 5, W2C2 leaves
# cost 900: u = [0, 1, 3], v = [8, 6, 12, 2], smallest d = -2 at W1C3
#    loop W1C3 W2C3 W2C1 W1C1, theta = 15, W1C1 leaves
# cost 870: u = [0, 3, 3], v = [6, 6, 10, 2], smallest d = 2 at W1C1   <- all d >= 0
# [0, 15, 15, 0]
# [20, 0, 20, 0]
# [0, 10, 0, 20]

Each pivot lowers the cost by $|d| \cdot \theta$: $6 \cdot 10 = 60$, then $5 \cdot 5 = 25$, then $2 \cdot 15 = 30$. The final $u, v$ are an optimal dual solution, and the dual objective $\sum_i s_i u_i + \sum_j d_j v_j = 0 + 120 + 90 + 120 + 150 + 350 + 40 = 870$ matches the primal cost (strong duality).

Key takeaways

MODI computes $u_i + v_j = c_{ij}$ on basic cells and reduced costs $d_{ij} = c_{ij} - u_i - v_j$ elsewhere.
All $d_{ij} \ge 0$ means optimal. Otherwise, the most negative cell enters along its unique stepping-stone loop.
$\theta$ is the smallest allocation on the minus cells of the loop, and the cell that reaches 0 leaves.
MODI is the simplex method (network simplex) written on the transportation tableau. $u, v$ are the dual prices.

Practice

(Beginner) In the final solution, compute $d_{ij}$ for all six empty cells by hand and confirm they are nonnegative.

(Intermediate) The route W2-C4 becomes cheaper: 4 instead of 7. Use the final $u, v$ to decide whether the plan is still optimal. If not, perform one more MODI iteration by hand and check the result with the code.

(Interview) "How do you handle degeneracy in MODI, when the starting solution has fewer than $m + n - 1$ positive cells?" Hint: add zero allocations ($\varepsilon$ cells) in cells that do not create a loop, so that the $u, v$ system can still be solved.

Assignment Problem

Why it matters. Assigning $n$ workers to $n$ jobs, machines to tasks, drivers to ride requests or radar readings to tracked objects at minimum total cost is one of the most frequently solved optimization problems. There are $n!$ possible assignments ($3.6$ million for $n = 10$), but polynomial algorithms solve it exactly.

Theory. With $x_{ij} = 1$ if worker $i$ does job $j$:

$$ \min \sum_{i,j} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_j x_{ij} = 1 \ \ \forall i, \quad \sum_i x_{ij} = 1 \ \ \forall j, \quad x_{ij} \in {0, 1}. $$

It is a transportation problem with all supplies and demands equal to 1, so the matrix is TU and the LP relaxation $0 \le x \le 1$ already has integer optimal vertices.
It is highly degenerate: a basis has $2n - 1$ cells but only $n$ of them are positive. MODI wastes many pivots on $\theta = 0$, which is why special algorithms (Hungarian, auction, shortest augmenting paths) are used.
Variants:

Situation	Modeling trick
More jobs than workers (or the reverse)	add dummy workers (or jobs) with cost 0
Worker $i$ cannot do job $j$	$c_{ij} = M$ (large) or remove the cell
Maximize total profit	minimize $-p_{ij}$, or $\max p - p_{ij}$
Worker can take $k$ jobs	duplicate the worker $k$ times (or use transportation)
Minimize the worst assignment cost	bottleneck assignment: binary search on a threshold plus bipartite matching

Example: technicians and jobs (hours needed).

import itertools
import numpy as np
from scipy.optimize import linprog, linear_sum_assignment

C = np.array([[10, 19, 8, 15],            # rows: technicians T1..T4
              [10, 18, 7, 17],            # columns: jobs J1..J4 (hours)
              [13, 16, 9, 14],
              [12, 19, 8, 18]])
n = len(C)

# 1. Brute force: n! permutations (24 here, 3.6 million for n = 10)
best = min(itertools.permutations(range(n)), key=lambda p: sum(C[i, p[i]] for i in range(n)))
print("brute force :", [f"T{i + 1}->J{j + 1}" for i, j in enumerate(best)],
      sum(int(C[i, best[i]]) for i in range(n)))

# 2. LP relaxation with 0 <= x <= 1: integral because the matrix is TU
A_eq = [np.kron(np.eye(n)[i], np.ones(n)) for i in range(n)]    # each row once
A_eq += [np.kron(np.ones(n), np.eye(n)[j]) for j in range(n)]   # each column once
res = linprog(C.ravel(), A_eq=A_eq, b_eq=np.ones(2 * n), bounds=(0, 1), method="highs")
print("LP solution :", np.round(res.x).reshape(n, n).astype(int).tolist(), res.fun)

# 3. Library routine (a shortest augmenting path variant of the Hungarian method)
rows, cols = linear_sum_assignment(C)
print("scipy       :", [(int(r) + 1, int(c) + 1) for r, c in zip(rows, cols)],
      int(C[rows, cols].sum()))
# brute force : ['T1->J4', 'T2->J1', 'T3->J2', 'T4->J3'] 49
# LP solution : [[0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]] 49.0
# scipy       : [(1, 4), (2, 1), (3, 2), (4, 3)] 49

Note that job J3 is the cheapest job for every technician. The greedy idea "give everyone their favorite job" fails immediately, because only one technician can take J3.

Key takeaways

The assignment problem pairs $n$ agents with $n$ tasks, each used exactly once, at minimum total cost.
It is a transportation problem with unit supplies and demands, so its LP relaxation is integral.
Its extreme degeneracy makes general simplex methods inefficient, so specialized $O(n^3)$ algorithms are used.
Dummies handle unequal sizes, large costs block pairs, and negation turns profit maximization into minimization.

Practice

(Beginner) Add a fifth job J5 with times $(11, 12, 13, 14)$ for T1-T4. How do you make the problem square, and what does the dummy's assignment mean?

(Intermediate) Solve the bottleneck version (minimize the longest single job) of the example by trying thresholds $t$ in increasing order and checking whether a perfect matching exists using only cells with $c_{ij} \le t$.

(Interview) "A ride-hailing app must match 2,000 drivers to 1,500 ride requests every few seconds. How would you model and solve it?" Hint: rectangular assignment with dummy requests, forbidden pairs beyond a pickup radius, and a fast solver (shortest augmenting path or auction), possibly warm-started from the previous round.

Hungarian Algorithm

Why it matters. The Hungarian algorithm (Kuhn, 1955, named after the Hungarian mathematicians König and Egerváry) was one of the first polynomial algorithms for a combinatorial problem and is still a standard method for assignment. It is a primal-dual algorithm: it adjusts dual prices until an assignment that uses only zero reduced-cost cells exists.

Theory. Two facts drive the method:

Reduction does not change the optimal assignment. Subtracting a constant from a whole row (or column) lowers the cost of every assignment by that same constant, because every assignment uses exactly one cell per row. The subtracted constants add up to a lower bound on the optimal cost (they are the dual variables $u_i$, $v_j$).
König-Egerváry. The maximum number of independent zeros (no two in the same row or column) equals the minimum number of lines (rows or columns) that cover all zeros.

Steps.

Subtract each row's minimum from that row.
Subtract each column's minimum from that column.
Cover all zeros with the minimum number of lines. If $n$ lines are needed, an assignment using only zeros exists: read it off and stop.
Otherwise, let $\delta$ = the smallest uncovered entry. Subtract $\delta$ from every uncovered entry and add $\delta$ to every entry covered twice. Go to step 3.

Step 4 creates a new zero without making any entry negative, and it raises the lower bound by $\delta \times (\text{number of uncovered rows} - \text{number of covered columns}) > 0$. The method runs in $O(n^3)$ with careful bookkeeping.

Trace on the technician example.

  Original            Rows reduced         Columns reduced       (lines: T1, T3, J3)
      J1 J2 J3 J4     (-8,-7,-9,-8)        (-2,-7,-0,-5)
  T1  10 19  8 15      2 11  0  7           0  4  0  2   <- line
  T2  10 18  7 17      3 11  0 10           1  4  0  5
  T3  13 16  9 14      4  7  0  5           2  0  0  0   <- line
  T4  12 19  8 18      4 11  0 10           2  4  0  5
                                                  ^ line             3 lines < 4

  Adjust by 1          (lines: T3, J1, J3)  Adjust by 2           4 independent zeros
   0  4  1  2                                 0  2  1 [0]          T1 -> J4   15
   0  3  0  4                                [0] 1  0  2           T2 -> J1   10
   2  0  1  0   <- line                       4 [0] 3  0           T3 -> J2   16
   1  3  0  4                                 1  1 [0] 2           T4 -> J3    8
   ^     ^ lines        3 lines < 4                                total      49

Lower bound bookkeeping: rows $32$ + columns $14$ + first adjustment $1 \times (2 - 1)$ + second adjustment $2 \times (3 - 2)$ $= 49$, equal to the cost of the assignment, which proves it optimal.

def hungarian(C):
    """Textbook Hungarian method: reduce, cover zeros with fewest lines, adjust."""
    n = len(C)
    M = [row[:] for row in C]
    for row in M:                         # 1. row reduction
        rmin = min(row)
        row[:] = [x - rmin for x in row]
    for j in range(n):                    # 2. column reduction
        cmin = min(M[i][j] for i in range(n))
        for i in range(n):
            M[i][j] -= cmin
    while True:
        # 3. Maximum set of independent zeros = max matching on zero cells
        match_col = {}                    # column -> row
        def augment(i, seen):
            for j in range(n):
                if M[i][j] == 0 and j not in seen:
                    seen.add(j)
                    if j not in match_col or augment(match_col[j], seen):
                        match_col[j] = i
                        return True
            return False
        for i in range(n):
            augment(i, set())
        if len(match_col) == n:
            return {i: j for j, i in match_col.items()}, M
        # 4. Fewest lines covering all zeros (Koenig): alternating paths from free rows
        matched_rows = set(match_col.values())
        Zr, Zc = {i for i in range(n) if i not in matched_rows}, set()
        frontier = list(Zr)
        while frontier:
            i = frontier.pop()
            for j in range(n):
                if M[i][j] == 0 and j not in Zc:
                    Zc.add(j)
                    if j in match_col and match_col[j] not in Zr:
                        Zr.add(match_col[j])
                        frontier.append(match_col[j])
        cov_rows, cov_cols = set(range(n)) - Zr, Zc
        # 5. Adjust by the smallest uncovered value
        delta = min(M[i][j] for i in range(n) for j in range(n)
                    if i not in cov_rows and j not in cov_cols)
        for i in range(n):
            for j in range(n):
                if i not in cov_rows and j not in cov_cols:
                    M[i][j] -= delta      # uncovered: subtract
                elif i in cov_rows and j in cov_cols:
                    M[i][j] += delta      # covered twice: add
        print(f"{len(match_col)} lines (rows {sorted(r + 1 for r in cov_rows)}, "
              f"cols {sorted(c + 1 for c in cov_cols)}), adjust by {delta}")

C = [[10, 19, 8, 15], [10, 18, 7, 17], [13, 16, 9, 14], [12, 19, 8, 18]]
assign, M = hungarian(C)
print("final reduced matrix:", M)
print("assignment:", {f"T{i + 1}": f"J{j + 1}" for i, j in sorted(assign.items())},
      "cost", sum(C[i][j] for i, j in assign.items()))
# 3 lines (rows [1, 3], cols [3]), adjust by 1
# 3 lines (rows [3], cols [1, 3]), adjust by 2
# final reduced matrix: [[0, 2, 1, 0], [0, 1, 0, 2], [4, 0, 3, 0], [1, 1, 0, 2]]
# assignment: {'T1': 'J4', 'T2': 'J1', 'T3': 'J2', 'T4': 'J3'} cost 49

The code finds the minimum line cover the same way König's theorem is proved: from a maximum matching of zeros (Section 3.4), lines go on the matched rows that cannot be reached by alternating paths and on the columns that can.

Key takeaways

Row and column reductions keep the optimal assignment and build a lower bound (the dual solution).
If $n$ lines are needed to cover all zeros, an optimal assignment on zero cells exists.
Otherwise, subtract the smallest uncovered value from uncovered cells and add it to doubly covered cells, which creates new zeros and raises the bound.
The method is primal-dual and runs in $O(n^3)$. Library versions (for example scipy.optimize.linear_sum_assignment) use the shortest augmenting path form of the same idea.

Practice

(Beginner) Solve by hand: $\begin{bmatrix} 4 & 1 & 3 \ 2 & 0 & 5 \ 3 & 2 & 2 \end{bmatrix}$. How many adjustment steps are needed?

(Intermediate) Use hungarian to maximize total profit for the profit matrix $\begin{bmatrix} 7 & 5 & 8 \ 6 & 9 & 4 \ 8 & 6 & 7 \end{bmatrix}$ by converting it into a cost matrix. Check the answer by brute force.

(Interview) "Why does subtracting the smallest uncovered value and adding it at the intersections never destroy the zeros you need, and why does the algorithm terminate?" Hint: zeros covered once are unchanged, and each adjustment either increases the matching size or the set of reachable rows, while the lower bound strictly increases.

Transshipment Problem

Why it matters. Real distribution networks rarely ship directly from factories to customers. Goods go through hubs, cross-docks and ports, and sometimes from one hub to another. The transshipment problem adds these intermediate nodes. It can be solved as a min-cost flow, or converted back into a transportation problem so that tableau methods still apply.

Theory. Nodes are sources ($b_i > 0$), sinks ($b_i < 0$) and transshipment nodes ($b_i = 0$, or any node that can also relay goods). Arcs can connect any pair, including hub to hub. Without arc capacities, this is the uncapacitated min-cost flow problem:

$$ \min \sum c_{ij} x_{ij} \quad \text{s.t.} \quad Nx = b, \quad x \ge 0. $$

Conversion to a transportation problem. Let $B$ = total supply (a buffer large enough for everything to pass through any node).

Node type	Appears as	Supply / demand in the tableau
Pure source	row	its supply
Pure sink	column	its demand
Transshipment node	both a row and a column	row supply $B$, column demand $B$ (plus its own supply or demand, if any)

Cell (node $k$ row, node $k$ column) costs 0: units "shipped to itself" are buffer that was not used. Missing arcs get cost $M$.

Example: two factories, two hubs, three customers. Factory F1 makes 40 units and F2 makes 30. Customers D1, D2 and D3 need 25, 25 and 20. Goods can go through hubs H1 and H2, from H1 to H2, or directly from F1 to D1:

Arc	F1 $\to$ H1	F1 $\to$ H2	F2 $\to$ H1	F2 $\to$ H2	H1 $\to$ H2	H1 $\to$ D1	H1 $\to$ D2	H2 $\to$ D2	H2 $\to$ D3	F1 $\to$ D1
Unit cost	3	5	6	2	1	4	6	3	4	9

import numpy as np
from scipy.optimize import linprog

# Supply > 0, demand < 0, pure transshipment hubs = 0
SUPPLY = {"F1": 40, "F2": 30, "H1": 0, "H2": 0, "D1": -25, "D2": -25, "D3": -20}
ARCS = {("F1", "H1"): 3, ("F1", "H2"): 5, ("F2", "H1"): 6, ("F2", "H2"): 2,
        ("H1", "H2"): 1, ("H1", "D1"): 4, ("H1", "D2"): 6, ("H2", "D2"): 3,
        ("H2", "D3"): 4, ("F1", "D1"): 9}
nodes, arcs = list(SUPPLY), list(ARCS)

N = np.zeros((len(nodes), len(arcs)))
for k, (i, j) in enumerate(arcs):
    N[nodes.index(i), k], N[nodes.index(j), k] = 1, -1
res = linprog([ARCS[a] for a in arcs], A_eq=N, b_eq=[SUPPLY[v] for v in nodes],
              bounds=(0, None), method="highs")
print("minimum cost:", res.fun)
for (i, j), x in zip(arcs, res.x):
    if x > 1e-9:
        print(f"  {i} -> {j}: {x:4.0f} units at {ARCS[i, j]}")

# Same problem as a transportation problem: every hub is both a source and a
# destination with a buffer B = total supply (unused buffer = "ship to itself" at 0)
B = sum(s for s in SUPPLY.values() if s > 0)
sources = ["F1", "F2", "H1", "H2"]
dests = ["H1", "H2", "D1", "D2", "D3"]
sup = [SUPPLY[v] + (B if v.startswith("H") else 0) for v in sources]
dem = [-SUPPLY[v] + (B if v.startswith("H") else 0) for v in dests]
BIG = 10**6                               # forbidden cell (no arc)
cost = np.array([[0 if i == j else ARCS.get((i, j), BIG) for j in dests] for i in sources])
m, n = cost.shape
A = [np.kron(np.eye(m)[i], np.ones(n)) for i in range(m)]
A += [np.kron(np.ones(m), np.eye(n)[j]) for j in range(n)]
res2 = linprog(cost.ravel(), A_eq=A, b_eq=sup + dem, bounds=(0, None), method="highs")
print(f"transportation form: supplies {sup}, demands {dem}, cost {res2.fun}")
# minimum cost: 450.0
#   F1 -> H1:   40 units at 3
#   F2 -> H2:   30 units at 2
#   H1 -> H2:   15 units at 1             <- hub-to-hub transfer
#   H1 -> D1:   25 units at 4
#   H2 -> D2:   25 units at 3
#   H2 -> D3:   20 units at 4
# transportation form: supplies [40, 30, 70, 70], demands [70, 70, 25, 25, 20], cost 450.0

The direct route F1 $\to$ D1 (9) is never used, because F1 $\to$ H1 $\to$ D1 costs only $3 + 4 = 7$. Some of F1's goods even travel F1 $\to$ H1 $\to$ H2 $\to$ D2/D3. A model that only allowed factory-to-hub-to-customer routes would miss this.

Key takeaways

Transshipment adds intermediate nodes that relay goods. It is the uncapacitated min-cost flow problem.
Solve it directly as a network flow, or convert it to a transportation problem by giving each relay node a buffer $B$ as both supply and demand.
Self-shipments of a relay node at cost 0 represent buffer capacity that was not used.
Allowing hub-to-hub and multi-stop routes can find plans that a two-stage model would miss.

Practice

(Beginner) Compute the cost of the plan "F1 ships 25 directly to D1 and 15 via H1 to D2, F2 ships everything via H2" and compare it with 450.

(Intermediate) Hub H2 can handle at most 50 units. Add the restriction with node splitting (Section 3.5) and re-solve. How much does the limit cost?

(Interview) "When would you prefer the transportation-tableau conversion over a direct network flow solver, and vice versa?" Hint: think about hand calculation and teaching versus model size, sparsity (most cells of the converted tableau are forbidden) and available software.

← Operative Research