Ford-Fulkerson Deep Dive: the augmenting-path method, why path choice decides running time, non-termination and measured augmentation counts

The method and why it is correct

A flow network is a directed graph with a source s, a sink t and a non-negative capacity c(u,v) on each edge. A flow assigns f(u,v) to each edge so that two invariants hold: capacity (0 ≤ f ≤ c) and conservation (every vertex except s and t has equal inflow and outflow). The value of the flow is the net amount leaving s.

The residual graph records what you could still change. A forward residual edge u→v has capacity c(u,v) − f(u,v): room to push more. A backward residual edge v→u has capacity f(u,v): flow you could cancel. An augmenting path is any s-t path in the residual graph, and its bottleneck is the smallest residual capacity along it. Pushing the bottleneck along the path preserves both invariants and raises the flow value by exactly that bottleneck.

ford_fulkerson(G, s, t, choose_path):
    f = 0 on every edge
    while (P = choose_path(residual(G, f), s, t)) exists:
        b = min residual capacity on P
        for each edge (u, v) in P:
            if (u, v) is a forward edge: f(u, v) += b
            else:                        f(v, u) -= b    # cancel earlier flow
    S = vertices reachable from s in residual(G, f)
    return f, cut(S, V - S)                              # value(f) == capacity(cut)

When the loop stops, the set S of vertices still reachable from s defines a cut. Every edge leaving S is saturated (otherwise its head would be reachable) and every edge entering S carries zero flow (otherwise its backward residual edge would extend S). So the flow equals the cut's capacity, and since no flow can exceed any cut, both are optimal. That one argument is the method's correctness proof, and it holds no matter which rule picked the paths.

An implementation with paired residual edges

The representation decides how cheap each step is. Store every edge next to its residual twin, so edge e and e ^ 1 are partners and pushing flow is two array updates, with no dictionary lookups and no special cases for backward edges. Parallel edges and antiparallel pairs also work without merging.

from collections import deque

class Network:
    def __init__(self, n):
        self.adj = [[] for _ in range(n)]
        self.to, self.cap, self.orig = [], [], []

    def add_edge(self, u, v, c):
        # edge e and its residual twin e ^ 1 are stored side by side
        for a, b, k in ((u, v, c), (v, u, 0)):
            self.adj[a].append(len(self.to))
            self.to.append(b); self.cap.append(k); self.orig.append(k)

def augment(net, path):
    b = min(net.cap[e] for e in path)
    for e in path:
        net.cap[e] -= b
        net.cap[e ^ 1] += b
    return b

def ford_fulkerson(net, s, t, find):
    flow = rounds = 0
    while (path := find(net, s, t)) is not None:
        flow += augment(net, path)
        rounds += 1
    return flow, rounds

def bfs_path(net, s, t, delta=1):
    parent, q = {s: None}, deque([s])
    while q and t not in parent:
        u = q.popleft()
        for e in net.adj[u]:
            v = net.to[e]
            if v not in parent and net.cap[e] >= delta:
                parent[v] = e
                q.append(v)
    if t not in parent:
        return None
    path, v = [], t
    while v != s:
        path.append(parent[v]); v = net.to[parent[v] ^ 1]
    return path[::-1]

Swap the deque for a stack and you have the depth-first rule. Replace it with a max-heap keyed on the best bottleneck seen so far (Dijkstra with min instead of plus) and you have the fattest-path rule. Wrap bfs_path in a loop over a threshold delta and you have capacity scaling:

def capacity_scaling(net, s, t):
    flow = rounds = 0
    delta = 1 << max(net.cap).bit_length()
    while delta >= 1:
        while (path := bfs_path(net, s, t, delta)) is not None:
            flow += augment(net, path)
            rounds += 1
        delta //= 2
    return flow, rounds

Worked example

Take the six-vertex network in the diagram. The depth-first rule, running on our adjacency order, found four augmenting paths:

After round 4 the residual graph reaches only s and b. The cut edges are s→a (10) and b→d (9), total 19, matching the flow. Note a→b: it crosses the cut from the T side back into S, so it contributes nothing to the cut's capacity, and the proof says it must carry zero flow, which it does. Breadth-first search also needed four rounds here, fattest-path needed three, and capacity scaling five. All found the same value, but different flows: fattest-path routed 6 units along s→a→d→c→t, where the depth-first rule routed 5 there and 1 along s→a→d→t. Max flow values are unique; the flows that reach them usually are not.

Any final flow also decomposes into at most m paths (plus cycles, if any). Following positive-flow edges from s and peeling off each path's bottleneck recovers a list such as {s→b→d→t: 9, s→a→c→t: 4, s→a→d→c→t: 5, s→a→d→t: 1}. That decomposition is what you hand to people who need routes, not edge totals.

How many augmentations: the rule decides

With integer capacities every bottleneck is at least 1, so the method stops after at most |f*| augmentations, where |f*| is the max-flow value. Each round costs O(m) for the search, giving O(m·|f*|). That bound is pseudo-polynomial, and it can actually be reached. Take the diamond s→a, s→b, a→t, b→t with capacity C and a middle edge a→b of capacity 1. A rule that alternates s→a→b→t and s→b→a→t (the second path cancels the middle edge's flow) gains exactly 1 per round. Simulating that choice gave 20 rounds for C = 10, 2,000 for C = 1,000 and 2,000,000 for C = 10⁶. Two shortest paths would have finished in two rounds.

Each rule has its own guarantee:

• Shortest path (Edmonds-Karp). Residual distances from s never shrink, and each edge can be the bottleneck only O(n) times, so the method needs O(nm) augmentations, or O(nm²) time, independent of capacities.
• Fattest path. If R is the flow still missing, the residual graph's max flow is R and it splits into at most m paths, so one path carries at least R/m. The fattest path therefore removes at least a 1/m fraction of what remains. After k rounds R ≤ |f*|(1 − 1/m)^k ≤ |f*|·e^−k/m, and with integers it stops within about m·ln|f*| + 1 rounds.
• Capacity scaling. At the end of a phase with threshold Δ, no residual path has bottleneck ≥ Δ, so the remaining flow is below m·Δ. Every augmentation in the next phase carries at least Δ/2, so a phase needs at most 2m rounds. With log U phases this gives O(m log U) augmentations and O(m² log U) time.

Measured on random graphs (seed 7, source vertex 0, sink n − 1), counting augmentations:

The worst-case bounds are far above these counts (m·ln|f*| is about 76,000 for the second row, yet fattest-path needed 10). The ranking is still useful: depth-first search grows with the capacities, and the other three barely notice them. Remember that the counts measure rounds, not time. A fattest-path round costs O(m log n) because of its heap, while a BFS round costs O(m).

When the method does not terminate

Termination depends on capacities being integers, or rationals, which you can scale to integers. With irrational capacities a bad path rule may never stop. In the commonly cited version of Uri Zwick's 1995 construction, three special edges get capacities 1, r and 1 with r = (√5 − 1)/2, which satisfies r² = 1 − r. Every other edge gets a large capacity M. A carefully repeated sequence of four augmenting paths leaves the special edges with residuals of the form r^k, r^k+1 and 0, so the cycle can repeat forever. The flow converges to 3 + 2r ≈ 4.24 even though the true maximum is 2M + 1. The method does not just run long; it converges to the wrong answer.

Floating-point capacities cause the practical version of the same problem. Bottlenecks like 1e-17 keep the loop alive, and rounding can leave a "saturated" edge with residual 1e-16, which corrupts the cut. Breadth-first selection terminates regardless of the numbers, since its bound counts edges and not capacities, but for exact answers you should still scale to integers (cents, megabits, milliseconds) before solving.

Integrality and what it buys

The method also yields a structural theorem: if every capacity is an integer, some maximum flow is integral, and Ford-Fulkerson finds one, because every bottleneck it ever pushes is an integer. Much of its modelling power comes from this fact. Bipartite matching (unit capacities from s to the left side, left to right, and right side to t) cannot end up "half-assigning" a worker. Edge-disjoint path counting cannot split a path in two. Unit-capacity networks also make the O(m·|f*|) bound practical, since |f*| ≤ n. For large matchings, though, Hopcroft-Karp is the specialised tool.

Operational guidance

Re-solving after a change. The residual graph is reusable state. If a capacity rises, keep the current flow, add the new room to the forward residual edge, and keep augmenting; typically only a few rounds are needed. If a capacity falls below the flow on that edge, the excess x has to come off: lower that edge's flow by x, cancel x units along a flow-carrying path from s to its tail and from its head to t, then resume augmenting, which may win the flow back through another route. Maximum flow in practice walks through this re-solve workflow.

Validate every result. Check capacity and conservation on each edge, then check that the cut read from the residual graph has capacity equal to the flow value. That check runs in O(m), needs no reference solver, and catches every bug that matters.

Engineering details. Use an iterative search (a recursive depth-first search hits CPython's default recursion limit of 1,000 frames on long paths). Use 64-bit capacities and an explicit value for "infinity" that cannot overflow when summed. Avoid per-round allocation of the parent array. Once graphs reach millions of edges, switch to Dinic or push-relabel. Ford-Fulkerson remains the method you can explain, verify and modify.

Failure modes

• Pseudo-polynomial blow-up: plain depth-first search on large capacities, as in the 18,993-round row above. Fix: use BFS or scaling.
• Missing reverse edges: augmenting without residual twins gives a greedy answer below the maximum. On the worked example, pushing 8 units along s→a→d→t first and never cancelling them ended at 16 instead of 19 in our test.
• Merging antiparallel edges incorrectly: u→v and v→u are both real edges. Pairing each with its own twin avoids subtle double counting.
• Wrong cut side: computing reachability in the original graph instead of the residual graph reports every edge as cut.
• Float residue: tiny positive residuals keep the loop running or misplace the cut. Scale to integers.

Trade-offs

What to do next

1. Implement Network with paired edges and bfs_path, then reproduce the worked example: value 19, cut {s→a, b→d}.
2. Add the stack and heap variants and count rounds on your own graphs. Expect DFS to grow with capacities.
3. Write the O(m) validator (capacity, conservation, cut equals flow) and run it after every solve.
4. Build the C-capacity diamond and simulate the alternating paths to see the |f*| bound with your own eyes.
5. Scale any real-valued capacities to integers before solving.
6. Model one real problem, such as an assignment or a circulation with demands, and read the min cut as the bottleneck report.