Ford-Fulkerson Method

A Bit of History

As covered in Part 15, L. R. Ford Jr. and D. R. Fulkerson published this method in their 1956 RAND Corporation report "Maximal Flow Through a Network," directly motivated by a classified 1955 analysis of Soviet railway capacity. Their method wasn't just an algorithm — it came paired with the max-flow min-cut theorem's proof, making it one of the rare cases in this series where an algorithm and its correctness proof were published together as two faces of the same insight.

Working Principle

The Ford-Fulkerson method (deliberately not called a fully-specified "algorithm" — it leaves one choice open, discussed below) repeatedly finds any augmenting path from \(s\) to \(t\) in the current residual graph (Part 15's construction), pushes flow equal to that path's bottleneck capacity, and updates the residual graph accordingly. It terminates exactly when no augmenting path remains — at which point, by the max-flow min-cut argument from Part 15, the current flow is provably maximum.

def ford_fulkerson_pseudocode(capacity, source, sink, vertices):
    """capacity: dict[(u,v)] -> capacity. Builds residual capacities as it goes."""
    residual = dict(capacity)
    for u, v in list(capacity):
        residual.setdefault((v, u), 0)   # reverse residual edges start at 0

    def find_augmenting_path():
        # any path-finding method works here -- DFS, BFS, etc.
        parent = {source: None}
        stack = [source]
        while stack:
            u = stack.pop()
            if u == sink:
                break
            for v in vertices:
                if residual.get((u, v), 0) > 0 and v not in parent:
                    parent[v] = u
                    stack.append(v)
        return parent if sink in parent else None

    max_flow = 0
    while (parent := find_augmenting_path()):
        # find bottleneck capacity along the path
        path, v = [], sink
        while v is not None:
            path.append(v)
            v = parent[v]
        path.reverse()
        bottleneck = min(residual[(path[i], path[i+1])] for i in range(len(path)-1))

        for i in range(len(path) - 1):
            u, v = path[i], path[i+1]
            residual[(u, v)] -= bottleneck
            residual[(v, u)] += bottleneck   # enable "undoing" this flow later

        max_flow += bottleneck

    return max_flow

Worked Example

A small network: \(S \to A\) (cap 3), \(S \to B\) (cap 2), \(A \to T\) (cap 2), \(B \to T\) (cap 3), \(A \to B\) (cap 1).

flowchart LR
    S -->|3| A
    S -->|2| B
    A -->|2| T
    B -->|3| T
    A -->|1| B
            

First augmenting path \(S \to A \to T\), bottleneck \(\min(3,2)=2\): push 2 units, saturating \(A \to T\). Second augmenting path \(S \to B \to T\), bottleneck \(\min(2,3)=2\): push 2 more units. Third, using the leftover \(A \to B\) capacity: \(S \to A \to B \to T\), bottleneck \(\min(1, 1, 1) = 1\) (only 1 unit of \(S \to A\) capacity remains, and \(B \to T\) has 1 unit of capacity left after the second path). Total max flow: \(2 + 2 + 1 = 5\).

Complexity Analysis

Each augmenting path increases the flow by at least 1 (assuming integer capacities), and the maximum possible flow is bounded by the sum of capacities out of the source, \(F\):

$$\text{Time: } O(E \cdot F) \qquad \text{Space: } O(V + E)$$

This bound depends on the numeric value \(F\), not just the graph's size — a genuinely unusual complexity class called pseudo-polynomial, since \(F\) can be exponentially large relative to the number of bits needed to write the capacities down.

A Pathological Case

Implementation

from collections import defaultdict

def ford_fulkerson(capacity, source, sink):
    """capacity: dict[(u,v)] -> capacity (int). Returns max flow value."""
    residual = defaultdict(int)
    graph = defaultdict(set)
    for (u, v), cap in capacity.items():
        residual[(u, v)] += cap
        graph[u].add(v)
        graph[v].add(u)   # allow traversing reverse residual edges too

    def dfs_augmenting_path():
        parent = {source: None}
        stack = [source]
        while stack:
            u = stack.pop()
            if u == sink:
                return parent
            for v in graph[u]:
                if residual[(u, v)] > 0 and v not in parent:
                    parent[v] = u
                    stack.append(v)
        return None

    max_flow = 0
    while True:
        parent = dfs_augmenting_path()
        if parent is None:
            break
        path = []
        v = sink
        while v is not None:
            path.append(v)
            v = parent[v]
        path.reverse()

        bottleneck = min(residual[(path[i], path[i+1])] for i in range(len(path) - 1))
        for i in range(len(path) - 1):
            u, v = path[i], path[i + 1]
            residual[(u, v)] -= bottleneck
            residual[(v, u)] += bottleneck

        max_flow += bottleneck

    return max_flow

capacity = {("S","A"): 3, ("S","B"): 2, ("A","T"): 2, ("B","T"): 3, ("A","B"): 1}
print(ford_fulkerson(capacity, "S", "T"))   # 5

#include <unordered_map>
#include <unordered_set>
#include <vector>
#include <iostream>
using namespace std;

int fordFulkerson(unordered_map<string, unordered_map<string,int>>& residual,
                   unordered_map<string, unordered_set<string>>& graph,
                   string source, string sink) {
    int maxFlow = 0;

    while (true) {
        unordered_map<string, string> parent{{source, ""}};
        vector<string> stack{source};
        bool found = false;

        while (!stack.empty() && !found) {
            string u = stack.back(); stack.pop_back();
            if (u == sink) { found = true; break; }
            for (auto& v : graph[u]) {
                if (residual[u][v] > 0 && !parent.count(v)) {
                    parent[v] = u;
                    stack.push_back(v);
                }
            }
        }
        if (!parent.count(sink)) break;

        vector<string> path;
        for (string v = sink; !v.empty(); v = parent[v]) path.push_back(v);
        int bottleneck = INT_MAX;
        for (size_t i = path.size() - 1; i > 0; i--)
            bottleneck = min(bottleneck, residual[path[i]][path[i-1]]);

        for (size_t i = path.size() - 1; i > 0; i--) {
            residual[path[i]][path[i-1]] -= bottleneck;
            residual[path[i-1]][path[i]] += bottleneck;
        }
        maxFlow += bottleneck;
    }
    return maxFlow;
}

int main() {
    unordered_map<string, unordered_map<string,int>> residual = {
        {"S", {{"A",3},{"B",2}}}, {"A", {{"T",2},{"B",1}}}, {"B", {{"T",3}}}
    };
    unordered_map<string, unordered_set<string>> graph = {
        {"S", {"A","B"}}, {"A", {"S","T","B"}}, {"B", {"S","A","T"}}, {"T", {"A","B"}}
    };
    cout << "Max flow: " << fordFulkerson(residual, graph, "S", "T") << endl;  // 5
    return 0;
}

import java.util.*;

class FordFulkerson {
    static int maxFlow(Map<String, Map<String, Integer>> residual,
                        Map<String, Set<String>> graph, String source, String sink) {
        int flow = 0;

        while (true) {
            Map<String, String> parent = new HashMap<>();
            parent.put(source, null);
            Deque<String> stack = new ArrayDeque<>();
            stack.push(source);
            boolean found = false;

            while (!stack.isEmpty() && !found) {
                String u = stack.pop();
                if (u.equals(sink)) { found = true; break; }
                for (String v : graph.getOrDefault(u, Set.of())) {
                    if (residual.get(u).getOrDefault(v, 0) > 0 && !parent.containsKey(v)) {
                        parent.put(v, u);
                        stack.push(v);
                    }
                }
            }
            if (!parent.containsKey(sink)) break;

            List<String> path = new ArrayList<>();
            for (String v = sink; v != null; v = parent.get(v)) path.add(v);
            int bottleneck = Integer.MAX_VALUE;
            for (int i = path.size() - 1; i > 0; i--)
                bottleneck = Math.min(bottleneck, residual.get(path.get(i)).getOrDefault(path.get(i-1), 0));

            for (int i = path.size() - 1; i > 0; i--) {
                String u = path.get(i), v = path.get(i - 1);
                residual.get(u).merge(v, -bottleneck, Integer::sum);
                residual.computeIfAbsent(v, k -> new HashMap<>()).merge(u, bottleneck, Integer::sum);
            }
            flow += bottleneck;
        }
        return flow;
    }

    public static void main(String[] args) {
        Map<String, Map<String,Integer>> residual = new HashMap<>();
        residual.put("S", new HashMap<>(Map.of("A",3,"B",2)));
        residual.put("A", new HashMap<>(Map.of("T",2,"B",1)));
        residual.put("B", new HashMap<>(Map.of("T",3)));
        residual.put("T", new HashMap<>());

        Map<String, Set<String>> graph = Map.of(
            "S", Set.of("A","B"), "A", Set.of("S","T","B"),
            "B", Set.of("S","A","T"), "T", Set.of("A","B")
        );
        System.out.println("Max flow: " + maxFlow(residual, graph, "S", "T"));  // 5
    }
}

Real-World Applications

Exercises

1. Run Ford-Fulkerson by hand on the worked example using a different order of augmenting paths, and confirm you still reach the same maximum flow value of 5.
2. Explain, using the residual-graph argument from Part 15, why the algorithm's termination (no augmenting path left) guarantees optimality rather than just "no more obvious improvement."
3. Construct a small network with irrational or very large capacities and describe (without fully simulating) why an unlucky path choice could make convergence extremely slow.
4. Challenge: Modify the implementation to use BFS instead of DFS for finding augmenting paths, and compare the number of iterations needed on a larger random network — this is exactly the Edmonds-Karp refinement, previewed next.

Limitations