Max Flow in Depth: Ford-Fulkerson, Residual Graphs, Edmonds-Karp and the Minimum Cut, with Working Code

Flow networks and what a flow is

A flow network is a directed graph with a source s, a sink t, and a capacity c(u, v) of zero or more on each edge. A flow assigns a number f(u, v) to each edge and must satisfy two rules. The capacity constraint: 0 ≤ f(u, v) ≤ c(u, v), so no edge carries more than it can. Conservation: for every vertex other than s and t, flow in equals flow out, so nothing is created or stored in the middle. The value of the flow is the net amount leaving s, which by conservation equals the net amount arriving at t.

Think of a flow as s-to-t paths, each carrying some amount, stacked without overloading any edge. The problem is choosing the paths well, because a greedy choice can block better ones.

Why greedy fails and reverse edges fix it

Take a diamond: s connects to a and b, a connects to b, and both a and b connect to t, every edge with capacity 1. The maximum flow is 2: s-a-t and s-b-t. But a naive algorithm that repeatedly finds any path with spare capacity might first choose s-a-b-t. Now s-a, a-b and b-t are full; the only spare edges are s-b and a-t, and there is no forward path from s to t. Greedy stops at 1.

The fix is to allow the algorithm to undo earlier decisions. For every unit of flow on u→v, add a residual edge v→u with capacity equal to that flow. Pushing flow along v→u means cancelling flow on u→v. In the diamond, after s-a-b-t the residual graph contains s→b, b→a (the reverse of a→b) and a→t, so the path s-b-a-t exists. Pushing 1 along it cancels the a→b flow; the net result is exactly s-a-t plus s-b-t, and the value is 2.

Formally, the residual capacity is r(u, v) = c(u, v) - f(u, v) + f(v, u). The residual graph contains every edge with r > 0, and an augmenting path is any s-to-t path in it. Its bottleneck is the smallest residual capacity along the path, and pushing that much increases the flow value by the bottleneck without violating any constraint.

The Ford-Fulkerson method

Ford-Fulkerson is a method more than an algorithm, because it leaves open how to find the path:

flow = 0
while there is an augmenting path P from s to t in the residual graph:
    b = min residual capacity on P          # the bottleneck
    for each edge (u, v) on P:
        r(u, v) -= b                         # use forward capacity
        r(v, u) += b                         # make it undoable
    flow += b
return flow

With integer capacities each iteration raises the flow by at least 1, so the method terminates after at most |f*| iterations, where f* is the maximum flow, each costing O(E) for a search. That bound depends on the capacity values, not the graph size: a four-vertex graph with capacities of one million and one unlucky middle edge of capacity 1 can take two million iterations if the search keeps choosing paths through the middle edge in alternating directions. With irrational capacities a bad path sequence may never terminate. So choose a path rule whose bound does not depend on capacities.

Integrality is a valuable side effect. If all capacities are integers, every bottleneck is an integer, so the maximum flow found is integral on every edge. That is what lets flow solve combinatorial problems such as matching, where an edge must carry either 0 or 1.

Edmonds-Karp: shortest augmenting paths

Edmonds-Karp makes one choice: always augment along a shortest path, measured in edge count, found by breadth-first search. The key lemma is that the BFS distance from s to every vertex in the residual graph never decreases between iterations. Each augmentation saturates at least one edge on the path (the bottleneck edge), and an edge (u, v) that was critical cannot become critical again until the distance to u has grown by at least 2. Distances are bounded by V, so each edge is critical O(V) times, giving O(VE) augmentations of O(E) each: O(VE²) total, independent of the capacities.

The implementation below stores edges in flat arrays with each edge's reverse at the index xor 1, which is compact, cache-friendly and makes the reverse lookup a single bit flip. It runs as written on Python 3.

from collections import deque

class FlowNetwork:
    def __init__(self, n):
        self.n = n
        self.adj = [[] for _ in range(n)]   # adj[u] = list of edge ids
        self.to, self.cap = [], []           # edge e and its reverse e ^ 1

    def add_edge(self, u, v, c):
        self.adj[u].append(len(self.to)); self.to.append(v); self.cap.append(c)
        self.adj[v].append(len(self.to)); self.to.append(u); self.cap.append(0)

    def max_flow(self, s, t):
        flow = 0
        while True:
            parent_edge = [-1] * self.n
            parent_edge[s] = -2
            q = deque([s])
            while q and parent_edge[t] == -1:           # BFS in the residual graph
                u = q.popleft()
                for e in self.adj[u]:
                    v = self.to[e]
                    if self.cap[e] > 0 and parent_edge[v] == -1:
                        parent_edge[v] = e
                        q.append(v)
            if parent_edge[t] == -1:
                return flow                              # no augmenting path left
            path, v = [], t
            while v != s:
                e = parent_edge[v]
                path.append(e)
                v = self.to[e ^ 1]
            push = min(self.cap[e] for e in path)       # bottleneck
            for e in path:
                self.cap[e] -= push
                self.cap[e ^ 1] += push
            flow += push

    def min_cut(self, s):
        seen, q = {s}, deque([s])
        while q:
            u = q.popleft()
            for e in self.adj[u]:
                if self.cap[e] > 0 and self.to[e] not in seen:
                    seen.add(self.to[e]); q.append(self.to[e])
        cut = [(u, self.to[e]) for u in seen for e in self.adj[u]
               if e % 2 == 0 and self.to[e] not in seen]   # even ids are original edges
        return seen, cut

Note that cap holds residual capacity, not the original capacity. Because every added edge gets its own pair, the flow on an original edge e is always the residual capacity of its partner, cap[e ^ 1], even when the graph also has an explicit opposite edge v→u.

Worked example: six vertices, four augmentations

Use vertices s, a, b, c, d, t and nine edges: s→a 10, s→c 10, a→b 4, a→c 2, a→d 8, c→d 9, d→b 6, b→t 10, d→t 10. Running the code above, with edges added in that order, performs four augmentations:

After iteration 4 BFS from s can reach only c: s→a is saturated, and c's only outgoing edge c→d is saturated too. The algorithm returns 19. Paths grow from three edges to four, which is the non-decreasing-distance lemma in action.

Reading the minimum cut off the residual graph

An s-t cut splits the vertices into a set S containing s and a set T containing t; its capacity is the total capacity of original edges from S to T. Every unit of flow must cross from S to T, so any flow is at most any cut. The max-flow min-cut theorem says the best of each meet: the maximum flow value equals the minimum cut capacity. There is more on the theorem in the max-flow min-cut theorem article; the constructive half is the one you use in code.

When Ford-Fulkerson stops, let S be every vertex reachable from s in the residual graph. t is not in S, or there would still be an augmenting path. Every original edge from S to T must be saturated, otherwise its head would be reachable, and every edge from T back to S must carry zero flow, otherwise its reverse residual edge would make its tail reachable. So the net flow across the cut equals the cut's capacity, and both are optimal. In the example min_cut returns S = {s, c} and the edges s→a and c→d, with 10 + 9 = 19.

Often the cut is the answer you actually want: the set of links whose failure separates two sites, or, in image segmentation with pixels joined to foreground and background terminals, the lowest-energy labelling.

Modelling: turning problems into networks

Most of the skill in using max flow is in building the graph. The recurring patterns are few:

• Bipartite matching. Source to each left vertex with capacity 1, each allowed pair with capacity 1, each right vertex to sink with capacity 1. The max flow is the maximum matching size, and integrality guarantees each left vertex uses at most one edge. For large sparse instances the specialised Hopcroft-Karp algorithm runs in O(E√V).
• Multiple sources or sinks. Add a super-source with edges to each source, capacity equal to its supply, and a super-sink likewise.
• Vertex capacities. Split each vertex v into v_in and v_out joined by an edge whose capacity is the vertex limit; incoming edges go to v_in and outgoing edges leave v_out. This is also how you count vertex-disjoint paths: give every split edge capacity 1.
• Edge-disjoint paths. Give every edge capacity 1; the max flow is the number of edge-disjoint s-t paths, which by Menger's theorem equals the minimum number of edges whose removal disconnects t from s.
• Costs. Minimising cost per unit of flow is a different problem, covered in min-cost max flow.

Failure modes and implementation pitfalls

• Forgetting reverse edges, or giving them the original capacity instead of 0. Without them the algorithm is greedy and can return too little, whenever an early path blocks better ones; with the wrong initial capacity it returns too much.
• Shallow copies. In Java, int[][] residual = capacity.clone() copies the outer array only; both arrays share rows, so building the residual graph destroys the input capacities, and a second call on the same object returns 0. Copy each row, or build the residual structure separately.
• Adjacency matrices on sparse graphs. A V-by-V matrix makes each BFS O(V²) regardless of E, and costs V² memory. At 100,000 vertices that is ten billion cells. Use edge lists as in the code above.
• Parallel and antiparallel edges. Matrix implementations merge two u→v edges into one entry, which is fine if you add their capacities but wrong if you overwrite, and they conflate an original v→u edge with the reverse of u→v. Edge lists keep them separate.
• Overflow. Use a sentinel for infinity that exceeds the total source capacity but stays well inside the integer type; summing several 'infinite' capacities in 32 bits wraps negative.
• Floating-point capacities. Rounding leaves residuals like 1e-17 and the loop finds useless paths. Scale to integers or use an epsilon.

Performance and choosing an algorithm

Edmonds-Karp is the right first implementation: short, easy to verify, and adequate for graphs up to tens of thousands of edges. Its weakness is that each BFS finds one path. Dinic's algorithm builds a BFS level graph once per phase and then finds a blocking flow of many shortest paths, giving O(V²E) in general and O(E√V) on unit-capacity bipartite graphs. Push-relabel algorithms work locally on vertices, reach O(V³), and are often fastest on large dense instances.

Measure before switching: a simple Edmonds-Karp on edge arrays often beats a sloppy Dinic. And always verify results in tests by checking conservation on every vertex, capacity limits on every edge, and that the flow value equals the capacity of the cut you recover; together the three checks prove optimality.

What to do next

1. Type in the FlowNetwork class and reproduce the worked example: flow 19, cut {s, c}.
2. Build a graph where the blocking path is the unique shortest one: all capacities 1, edges s→a, a→b, b→t, a→c, c→d, d→t, s→e, e→f, f→b. Confirm the code returns 2, and 1 once you delete the reverse-edge update.
3. Write a checker that validates capacity, conservation and flow value equals cut capacity, and run it on random graphs.
4. Model a small bipartite matching problem and a vertex-capacity problem with the patterns above.
5. Implement Dinic's level graph and blocking flow, then benchmark both on a random graph with 10,000 vertices and 100,000 edges.
6. Review the BFS and DFS fundamentals if the level-graph argument feels shaky.