Push-Relabel Algorithm, in depth: building a highest-label solver with gap and global relabelling, then proving its answer

What the solver must return

Before writing code, decide what the solver must return, because that decides whether you need the second phase.

Phase 1 ends when no vertex below height n is active. At that point the cut is final, but some excess may still be stranded inside the graph. The solver below simply keeps going until every vertex is inactive. Excess that cannot reach t climbs above n and drains back to s. You then get a valid flow on every edge without a separate phase-2 routine. That is simpler and costs a little extra time on graphs with a lot of stranded excess.

Data layout

Store edges in flat arrays: to[e] and cap[e], with each forward edge at an even index and its residual partner at the next odd index. The partner of e is then e ^ 1, with no lookup. Each vertex keeps a list of its edge ids, or a CSR offset range for large graphs. Active vertices sit in buckets indexed by height. A top pointer tracks the highest non-empty bucket, and a count of vertices per height makes the gap test O(1).

A complete highest-label solver

from collections import deque

class MaxFlow:
    def __init__(self, n):
        self.n, self.head, self.to, self.cap = n, [[] for _ in range(n)], [], []

    def add_edge(self, u, v, c):
        for a, b, w in ((u, v, c), (v, u, 0)):        # edge e and partner e ^ 1
            self.head[a].append(len(self.to)); self.to.append(b); self.cap.append(w)

    def _global_relabel(self, s, t):
        n = self.n; h = [2 * n] * n
        for root, base in ((t, 0), (s, n)):           # exact distance to t, else n + distance to s
            if h[root] != 2 * n: continue
            h[root] = base; q = deque([root])
            while q:
                v = q.popleft()
                for e in self.head[v]:
                    u = self.to[e]
                    if h[u] == 2 * n and self.cap[e ^ 1] > 0:   # residual arc u -> v
                        h[u] = h[v] + 1; q.append(u)
        h[s] = n
        return h

    def max_flow(self, s, t):
        n, to, cap, head = self.n, self.to, self.cap, self.head
        ex = [0] * n
        for e in head[s]:                              # saturate source arcs
            ex[to[e]] += cap[e]; ex[s] -= cap[e]; cap[e ^ 1] += cap[e]; cap[e] = 0
        work, budget = 0, 6 * n + len(to)
        def rebuild():
            h = self._global_relabel(s, t)
            count = [0] * (2 * n + 1); bucket = [[] for _ in range(2 * n + 1)]
            for v in range(n):
                count[h[v]] += 1
                if v not in (s, t) and ex[v] > 0 and h[v] < 2 * n: bucket[h[v]].append(v)
            return h, count, bucket, [0] * n, 2 * n
        h, count, bucket, cur, top = rebuild()
        while top >= 0:
            if not bucket[top]: top -= 1; continue
            u = bucket[top].pop()
            if h[u] != top or ex[u] == 0: continue     # stale entry
            while ex[u] > 0:
                if cur[u] == len(head[u]):             # relabel
                    old = h[u]
                    h[u] = min([h[to[e]] + 1 for e in head[u] if cap[e] > 0] + [2 * n])
                    count[old] -= 1; count[h[u]] += 1; cur[u] = 0
                    work += 12 + len(head[u])
                    if count[old] == 0 and 0 < old < n:   # gap: O(n) scan here
                        for v in range(n):
                            if old < h[v] < n:
                                count[h[v]] -= 1; h[v] = n + 1; count[n + 1] += 1; cur[v] = 0
                                if v != u and ex[v] > 0: bucket[n + 1].append(v)
                        top = max(top, n + 1)
                    if h[u] >= 2 * n: break
                    continue
                e = head[u][cur[u]]; v = to[e]
                if cap[e] > 0 and h[u] == h[v] + 1:   # push
                    d = min(ex[u], cap[e])
                    cap[e] -= d; cap[e ^ 1] += d; ex[u] -= d; ex[v] += d
                    if v not in (s, t) and ex[v] == d: bucket[h[v]].append(v)
                else:
                    cur[u] += 1
            if h[u] < 2 * n:                           # pushes may have landed above old top
                top = max(top, h[u])
                if ex[u] > 0: bucket[h[u]].append(u)
            if work > budget:
                work = 0; h, count, bucket, cur, top = rebuild()
        return ex[t]

Three details carry the correctness. The current-arc pointer resets only on relabel. A vertex is queued when its excess goes from zero to positive, and s and t are never queued. Stale bucket entries are skipped by checking the height on pop. The gap lift here scans all n vertices, which is O(n) per gap. Production solvers keep a doubly linked list of all vertices per height, so the lift touches only the vertices it moves.

Worked example: watching the gap fire

Use four vertices: s, a, b, t with arcs s-a capacity 4, a-b capacity 4 and b-t capacity 1. The maximum flow is obviously 1. The interesting part is what happens to the 3 units that cannot get through.

Instrumented with counters, the solver records 4 pushes, 2 relabels and 1 gap on this run, and returns 1. Without the gap heuristic, a and b would bounce the 3 units between them, climbing one level per relabel until they rose above s. On a graph with n vertices that climb is up to n levels per vertex, and it is the most common reason a textbook implementation is slow. Note also that after step 4 no active vertex is below height n, so the cut {s, a, b} versus {t} is final. That is where phase 1 would have stopped.

Tuning knobs

The 6n + m budget is a reasonable starting point chosen for this page, not a published constant. Tune it on graphs drawn from your own workload. Count pushes, relabels, gaps and global relabels. A high relabel-to-push ratio means heights are drifting, and global relabels should run more often.

Extracting the cut and verifying the answer

A max-flow answer is cheap to verify, and you should verify it rather than trust it. Keep a copy of the original capacities before solving. Then check capacity bounds, conservation at every vertex except s and t, and finally that the residual cut has exactly the flow value. If all three pass, the max-flow min-cut theorem says the answer is optimal, as explained in the max-flow min-cut theorem article.

def verify(g, orig, s, t, value):
    bal = [0] * g.n
    for e in range(0, len(g.to), 2):                  # forward edges only
        f = orig[e] - g.cap[e]
        assert 0 <= f <= orig[e], "capacity violated"
        bal[g.to[e ^ 1]] -= f; bal[g.to[e]] += f
    assert all(bal[v] == 0 for v in range(g.n) if v not in (s, t)), "conservation"
    assert bal[t] == value
    seen = {s}; q = [s]                               # S side of the min cut
    while q:
        u = q.pop()
        for e in g.head[u]:
            if g.cap[e] > 0 and g.to[e] not in seen:
                seen.add(g.to[e]); q.append(g.to[e])
    assert t not in seen
    cut = sum(orig[e] for e in range(0, len(g.to), 2)
              if g.to[e ^ 1] in seen and g.to[e] not in seen)
    assert cut == value, "cut differs from flow"
    return seen                                       # the minimum cut's S side

The returned set is the minimum cut, which is all that segmentation and project-selection users need. For differential testing, generate thousands of small random graphs, including zero capacities, parallel arcs, and graphs where t is unreachable. Compare the value against a slow, obviously correct Edmonds-Karp, and run verify on every result. The code on this page passed that test on 5,000 random graphs with 2 to 14 vertices. Keep the harness in CI, because heuristic code regresses quietly.

Operational guidance

• Use a library unless you need the internals. The Boost Graph Library's push_relabel_max_flow and Google OR-Tools' max-flow solver are push-relabel implementations. Write your own when you must modify the algorithm, embed it, or run it on unusual hardware.
• Profile on real graphs. Vision grids, bipartite assignment and network topologies behave very differently. Record pushes, relabels and gaps alongside wall time.
• Bound the run. Put a wall-clock or work limit on the solver in services, and return a typed error rather than hanging a request.
• Parallelise by vertex. Push and relabel touch only one vertex's arcs and its neighbours' heights. That locality is why push-relabel is the usual base for multi-threaded and GPU max-flow. The cost is synchronisation on shared excess and heights, so measure before assuming a speed-up.
• Compare against Dinic. On many sparse graphs Dinic's algorithm is competitive and simpler. Keep both behind one interface and pick per workload.

Failure modes

• Gap lift without re-queueing. Lifted vertices that hold excess must go into the new bucket. If they don't, their old entries are discarded as stale and the excess is lost. The verifier reports this as a conservation failure.
• Gap applied at height n or above. The test must be 0 < old < n. Above n, heights measure distance back to s, and lifting there breaks the labelling.
• Global relabel that ignores s-side vertices. A vertex that cannot reach t still needs a valid height. Use n plus its distance to s, or excess is stranded.
• Stale height counts after a global relabel. Rebuild counts, buckets and current arcs together, as rebuild() does.
• Overflow. Saturating every source arc can pile the total source capacity onto one vertex. Use 64-bit excess.
• Not raising top after a relabel. A relabelled vertex can push to neighbours above the old top, and those vertices are then never popped. An earlier draft of the code on this page had exactly this bug. The differential test caught it on random graphs even though the hand-traced examples came out right.

Trade-offs

What to do next

1. Add push, relabel and gap counters to max_flow, run the s-a-b-t graph, and confirm 4 pushes, 2 relabels and 1 gap.
2. Add verify and an Edmonds-Karp oracle, then run 5,000 random graphs.
3. Replace the O(n) gap scan with per-height linked lists and measure on a large grid graph.
4. Sweep the global relabel budget from n to 20n + m and plot time against budget on your own workload.
5. Add a phase-1 stop flag for cut-only callers and check the cut matches the full run.
6. Read min-cost flow next. Cost scaling reuses push and relabel with reduced costs in place of heights.