Ford–Fulkerson, residual graphs and reverse edges

Lesson · 13 min read · Graphs · Maximum flow and minimum cut

Push flow along any path with spare capacity, and let later paths undo earlier choices through reverse edges.

Greedy is not enough

The natural idea is to find any path from s to t with spare capacity on every edge, push as much as the tightest edge allows, and repeat. It can get stuck. Below, every capacity is 1. If the first path chosen is s → a → b → t, no second path exists among the remaining spare capacities, yet the maximum flow is 2 (s → a → t and s → b → t).

After pushing along s → a → b → t, the edges a→b and b→t are full and greedy stops at 1.

The residual graph

The fix is to allow later paths to undo earlier flow. For every edge u → v carrying flow f, the residual graph contains:

a forward edge u → v with the capacity still unused, c − f;
a reverse edge v → u with capacity f: pushing along it cancels some of the flow on u → v.

An augmenting path is a path from s to t in the residual graph. In the trap above, s → b → a → t is one: it uses the reverse edge b → a, which cancels the unit on a → b and reroutes it. The result is the two disjoint paths, flow 2.

The poor first path, then an augmenting path that travels the reverse edge b → a.

Step 1 / 4

total0
Every edge is labelled flow/capacity, all flows 0. Look for a path from s to t with spare capacity on every step.
Step 2 / 4

total1
Augmenting path s → a → b → t, bottleneck 1: push 1 more unit. Total flow 1.
Step 3 / 4

total2
Augmenting path s → b → a → t, bottleneck 1: push 1 more unit. The step b→a runs against existing flow: it travels a reverse edge, cancelling flow on a→b and sending it elsewhere. Total flow 2.
Step 4 / 4

total2
No augmenting path is left: the maximum flow is 2. The nodes still reachable from s through spare capacity (shaded) form one side of a minimum cut; the cut edges (red) are full and their capacities add to 2.

carrying flowflow cancelled through a reverse edge

The algorithm

Ford–Fulkerson

maxFlow(graph, s, t)

Start with zero flow on every edge.
While the residual graph has a path from s to t: let b be the smallest residual capacity on it (the bottleneck).
Push b along the path: add b to each forward edge used, subtract b from the flow on each edge whose reverse was used.
Return the total pushed.

Which path to take matters for speed. With arbitrary paths and large capacities the number of rounds can be as large as the flow value itself. Edmonds–Karp always takes a shortest augmenting path, found by breadth-first search, which bounds the number of rounds by O(n · m) whatever the capacities.

Edmonds–Karp, step by step

On the network from the previous lesson. Each frame pushes along the shortest augmenting path; labels show flow/capacity.

Edmonds–Karp: breadth-first search for each augmenting path, then the minimum cut it leaves behind.

Step 1 / 6

total0
Every edge is labelled flow/capacity, all flows 0. Look for a path from s to t with spare capacity on every step.
Step 2 / 6

total4
Augmenting path s → a → b → t, bottleneck 4: push 4 more units. Total flow 4.
Step 3 / 6

total10
Augmenting path s → a → d → t, bottleneck 6: push 6 more units. Total flow 10.
Step 4 / 6

total14
Augmenting path s → c → d → t, bottleneck 4: push 4 more units. Total flow 14.
Step 5 / 6

total19
Augmenting path s → c → d → b → t, bottleneck 5: push 5 more units. Total flow 19.
Step 6 / 6

total19
No augmenting path is left: the maximum flow is 19. The nodes still reachable from s through spare capacity (shaded) form one side of a minimum cut; the cut edges (red) are full and their capacities add to 19.

carrying flowreachable from s at the endreverse edge used

The last frame is the theorem made concrete. The nodes still reachable from s through spare capacity form one side of a cut; every edge leaving that side is full, so the cut's capacity equals the flow. No flow can beat a cut, so this flow is maximum.

Implementation

Edmonds–Karp on an adjacency matrix of residual capacities, which handles reverse edges and parallel edges with no extra bookkeeping:

from collections import deque

def max_flow(n, edges, source, sink):
    res = [[0] * n for _ in range(n)]       # residual capacities
    for u, v, c in edges:
        res[u][v] += c                       # parallel edges add up
    total = 0
    while True:
        prev = [-1] * n
        prev[source] = source
        q = deque([source])
        while q and prev[sink] == -1:        # shortest augmenting path
            u = q.popleft()
            for v in range(n):
                if prev[v] == -1 and res[u][v] > 0:
                    prev[v] = u
                    q.append(v)
        if prev[sink] == -1:
            return total                     # no augmenting path left
        b, v = float('inf'), sink
        while v != source:                   # bottleneck
            b = min(b, res[prev[v]][v])
            v = prev[v]
        v = sink
        while v != source:                   # push, and open the reverse edge
            res[prev[v]][v] -= b
            res[v][prev[v]] += b
            v = prev[v]
        total += b

#include <algorithm>
#include <climits>
#include <queue>
#include <vector>

long long solve(int n, const std::vector<std::vector<int>>& edges, int source, int sink) {
    std::vector<std::vector<long long>> res(n, std::vector<long long>(n, 0));   // residual capacities
    for (const auto& e : edges) res[e[0]][e[1]] += e[2];                         // parallel edges add up
    long long total = 0;
    for (;;) {
        std::vector<int> prev(n, -1);
        prev[source] = source;
        std::queue<int> q;
        q.push(source);
        while (!q.empty() && prev[sink] == -1) {                                 // shortest augmenting path
            int u = q.front();
            q.pop();
            for (int v = 0; v < n; ++v)
                if (prev[v] == -1 && res[u][v] > 0) { prev[v] = u; q.push(v); }
        }
        if (prev[sink] == -1) return total;
        long long b = LLONG_MAX;
        for (int v = sink; v != source; v = prev[v]) b = std::min(b, res[prev[v]][v]);
        for (int v = sink; v != source; v = prev[v]) { res[prev[v]][v] -= b; res[v][prev[v]] += b; }
        total += b;
    }
}

import java.util.*;

class Solution {
    public long maxFlow(int n, int[][] edges, int source, int sink) {
        long[][] res = new long[n][n];                     // residual capacities
        for (int[] e : edges) res[e[0]][e[1]] += e[2];     // parallel edges add up
        long total = 0;
        while (true) {
            int[] prev = new int[n];
            Arrays.fill(prev, -1);
            prev[source] = source;
            ArrayDeque<Integer> q = new ArrayDeque<>();
            q.add(source);
            while (!q.isEmpty() && prev[sink] == -1) {     // shortest augmenting path
                int u = q.poll();
                for (int v = 0; v < n; v++)
                    if (prev[v] == -1 && res[u][v] > 0) { prev[v] = u; q.add(v); }
            }
            if (prev[sink] == -1) return total;
            long b = Long.MAX_VALUE;
            for (int v = sink; v != source; v = prev[v]) b = Math.min(b, res[prev[v]][v]);
            for (int v = sink; v != source; v = prev[v]) { res[prev[v]][v] -= b; res[v][prev[v]] += b; }
            total += b;
        }
    }
}

function maxFlow(n, edges, source, sink) {
  const res = Array.from({ length: n }, () => new Array(n).fill(0));   // residual capacities
  for (const [u, v, c] of edges) res[u][v] += c;                       // parallel edges add up
  let total = 0;
  for (;;) {
    const prev = new Array(n).fill(-1);
    prev[source] = source;
    const q = [source];
    for (let head = 0; head < q.length && prev[sink] === -1; head++) {   // shortest augmenting path
      const u = q[head];
      for (let v = 0; v < n; v++) if (prev[v] === -1 && res[u][v] > 0) { prev[v] = u; q.push(v); }
    }
    if (prev[sink] === -1) return total;
    let b = Infinity;
    for (let v = sink; v !== source; v = prev[v]) b = Math.min(b, res[prev[v]][v]);
    for (let v = sink; v !== source; v = prev[v]) { res[prev[v]][v] -= b; res[v][prev[v]] += b; }   // reverse edge
    total += b;
  }
}

use std::collections::VecDeque;

fn solve(n: i64, edges: &[[i64; 3]], source: i64, sink: i64) -> i64 {
    let (n, s, t) = (n as usize, source as usize, sink as usize);
    let mut res = vec![vec![0i64; n]; n];                 // residual capacities
    for e in edges {
        res[e[0] as usize][e[1] as usize] += e[2];        // parallel edges add up
    }
    let mut total = 0;
    loop {
        let mut prev = vec![usize::MAX; n];
        prev[s] = s;
        let mut q = VecDeque::from([s]);
        while let Some(u) = q.pop_front() {               // shortest augmenting path
            if prev[t] != usize::MAX {
                break;
            }
            for v in 0..n {
                if prev[v] == usize::MAX && res[u][v] > 0 {
                    prev[v] = u;
                    q.push_back(v);
                }
            }
        }
        if prev[t] == usize::MAX {
            return total;
        }
        let mut b = i64::MAX;
        let mut v = t;
        while v != s {
            b = b.min(res[prev[v]][v]);
            v = prev[v];
        }
        v = t;
        while v != s {
            res[prev[v]][v] -= b;
            res[v][prev[v]] += b;                         // reverse edge
            v = prev[v];
        }
        total += b;
    }
}

import qualified Data.Map.Strict as M
import Data.Sequence (Seq (..), (|>))
import qualified Data.Sequence as Q

solve :: Int -> [[Int]] -> Int -> Int -> Int
solve n edges s t = go (M.fromListWith (+) [((u, v), c) | [u, v, c] <- edges]) 0
  where
    cap res u v = M.findWithDefault 0 (u, v) res
    -- breadth-first search for a shortest augmenting path; returns predecessors
    bfs res = loop (Q.singleton s) (M.singleton s s)
      where
        loop Empty prev = prev
        loop (u :<| q) prev
          | M.member t prev = prev
          | otherwise =
              let next = [v | v <- [0 .. n - 1], not (M.member v prev), cap res u v > 0]
              in loop (foldl (|>) q next) (foldr (\v -> M.insert v u) prev next)
    go res total = case M.lookup t prev of
        Nothing -> total
        Just _ -> go res' (total + b)
      where
        prev = bfs res
        path = takeWhile (\(u, v) -> u /= v) [(prev M.! v, v) | v <- iterate (prev M.!) t]
        b = minimum [cap res u v | (u, v) <- path]
        res' = foldr (\(u, v) r -> M.insertWith (+) (v, u) b (M.insert (u, v) (cap r u v - b) r)) res path

The same code is the reference solution to Maximum flow through a network, where it is compiled and run against the tests in all six languages.

TimeO(n · m²)

SpaceO(n²)

With a matrix, each breadth-first search is O(n²), so O(n³ · m) in this form; fine for the few hundred nodes interview problems use. Dinic's algorithm, O(n² · m), is the next step up.