Hopcroft-Karp, in depth: layered BFS, disjoint shortest augmenting paths, the square-root phase bound, and an iterative implementation

Matchings and augmenting paths

A bipartite graph has vertices split into a left set L and a right set R, with every edge going between them. A matching is a set of edges sharing no endpoint; a vertex touched by the matching is matched, otherwise free. Typical readings: workers and tasks, pods and nodes that satisfy their constraints, students and project slots, reviewers and papers.

An alternating path uses unmatched and matched edges in turn. An augmenting path is an alternating path that starts at a free left vertex and ends at a free right vertex. It has one more unmatched edge than matched edges, so flipping every edge on it (matched becomes unmatched and vice versa) grows the matching by exactly one. Berge's theorem says a matching is maximum precisely when no augmenting path exists, which gives every augmenting-path algorithm its stopping rule.

Kuhn's algorithm finds one augmenting path per search, and each search can cost O(E), giving O(VE). Hopcroft-Karp's insight is to find many augmenting paths per pass, and specifically a maximal set of vertex-disjoint shortest ones. Restricting to shortest paths is what makes the analysis work: after each pass, the shortest augmenting path is strictly longer than before.

One phase: layer, then augment

The algorithm runs in phases. Each phase has two steps.

1. BFS layering. Start a breadth-first search simultaneously from every free left vertex at distance 0. From a left vertex, follow any edge to a right vertex. From a right vertex, if it is matched, follow its matched edge back to the left side, which places that left vertex one layer deeper. If it is free, an augmenting path has been found; record its length k and stop expanding beyond that layer. If no free right vertex is reachable, the matching is maximum and the algorithm ends.
2. DFS augmentation. From each free left vertex, search depth-first, moving only from layer i to layer i+1, and only accepting a free right vertex when the current left vertex sits on the final layer. When a path reaches a free right vertex, flip it. When a vertex turns out to be a dead end, mark it so no later search in this phase revisits it.

The layering guarantees every path found has length exactly k, the shortest possible. The dead-end marking and the fact that each edge is scanned through a per-vertex cursor make the whole DFS step O(E), so a phase costs O(E) in total. The paths found are vertex-disjoint because after an augmentation the flipped matched edges point backwards in the layering and cannot be used again in the same phase.

Implementation in Python

The implementation below uses 0-based left vertices 0..n_left-1, right vertices 0..n_right-1, and an adjacency list adj[u] of right neighbours. The DFS is iterative with an explicit stack, because a recursive version exceeds Python's default recursion limit of 1,000 on long augmenting paths.

from collections import deque

INF = float("inf")

def hopcroft_karp(n_left, n_right, adj, match_l=None):
    """Maximum bipartite matching. Returns (size, match_l, match_r, dist)."""
    match_l = list(match_l) if match_l else [-1] * n_left   # warm start allowed
    match_r = [-1] * n_right
    for u, v in enumerate(match_l):
        if v != -1:
            match_r[v] = u
    dist = [INF] * n_left

    def bfs():
        q = deque()
        for u in range(n_left):
            dist[u] = 0 if match_l[u] == -1 else INF
            if dist[u] == 0:
                q.append(u)
        limit = INF                       # layer where a free right vertex appears
        while q:
            u = q.popleft()
            if dist[u] > limit:
                continue                  # never build layers past the shortest length
            for v in adj[u]:
                w = match_r[v]
                if w == -1:
                    limit = min(limit, dist[u])
                elif dist[w] == INF:
                    dist[w] = dist[u] + 1
                    q.append(w)
        return limit

    def augment(root, limit, it):
        stack, via = [root], []           # via[i] = right vertex leading to stack[i+1]
        while stack:
            u = stack[-1]
            pushed = False
            while it[u] < len(adj[u]):
                v = adj[u][it[u]]
                it[u] += 1
                w = match_r[v]
                if w == -1:
                    if dist[u] == limit:  # free right vertex on the final layer
                        via.append(v)
                        for uu, vv in zip(stack, via):
                            match_l[uu], match_r[vv] = vv, uu
                        return True
                elif dist[w] == dist[u] + 1:
                    via.append(v)
                    stack.append(w)
                    pushed = True
                    break
            if not pushed:
                dist[u] = INF             # dead end for the rest of this phase
                stack.pop()
                if via:
                    via.pop()
        return False

    size = sum(1 for v in match_l if v != -1)
    while True:
        limit = bfs()
        if limit == INF:
            return size, match_l, match_r, dist
        it = [0] * n_left                 # per-vertex edge cursor, reset each phase
        for u in range(n_left):
            if match_l[u] == -1 and augment(u, limit, it):
                size += 1

Three details carry the complexity bound. The limit cutoff stops the BFS from building layers past the shortest augmenting length. The cursor it[u] means an edge examined once in a phase is never examined again in that phase. And setting dist[u] = INF on a dead end removes the vertex from the layered graph so other searches skip it. Drop any of the three and the code still returns a correct maximum matching, just more slowly, which is exactly why such bugs survive testing.

Worked example

Take left vertices a, b, c, d and right vertices 1, 2, 3, 4 with edges a-1, a-2, b-1, c-2, c-3, d-3, d-4. Process left vertices and adjacency lists in that order.

Phase 1. Every left vertex is free, so all four sit at distance 0. Every right vertex is free too, so the BFS finds augmenting paths of length 1 and sets limit = 0. The DFS from a takes edge a-1. The DFS from b tries vertex 1, now matched to a, but a is at distance 0, not 1, so b is a dead end. c takes 2, d takes 3. The matching is {a-1, c-2, d-3}, size 3, after one phase.

Phase 2. Only b is free. The BFS puts b at distance 0, reaches 1, follows the matched edge to a (distance 1), reaches 2 and then c (distance 2), reaches 3 and then d (distance 3), and finally reaches 4, which is free, so limit = 3. The DFS from b walks the same layers and flips the path b-1-a-2-c-3-d-4: four unmatched edges become matched and three matched edges become unmatched.

The result {b-1, a-2, c-3, d-4} is perfect, and phase 3's BFS finds no free left vertex, so the algorithm stops. Kuhn's algorithm on the same input would also succeed, but it would have spent one full search per left vertex; on large graphs, the batching of many short paths into one O(E) phase is the entire advantage.

Why the phase count is O(√V)

Why at most about 2√V phases? Two facts do the work.

• Lengths grow. After a phase that augments along a maximal set of vertex-disjoint shortest paths of length k, every remaining augmenting path has length greater than k. Since augmenting paths have odd length, the shortest length rises by at least 2 per phase.
• Few long paths remain. Let M be the current matching and M* a maximum one. The symmetric difference of M and M* contains at least |M*| - |M| vertex-disjoint augmenting paths for M. If each has length greater than √V, they together use more than (|M*| - |M|) √V vertices, and there are only V vertices, so |M*| - |M| < √V.

After √V phases the shortest augmenting path exceeds √V, so fewer than √V augmentations remain, and each later phase performs at least one. Total phases: at most about 2√V, each O(E), giving O(E √V). On dense graphs where E approaches V², that is O(V^2.5). In practice the phase count is usually far below the bound; sparse and random graphs commonly finish in a small number of phases, which is why measured run time often looks close to linear.

The certificate: a minimum vertex cover

A maximum matching is more useful with a certificate. König's theorem says that in a bipartite graph the size of a maximum matching equals the size of a minimum vertex cover, a set of vertices touching every edge. The last BFS of Hopcroft-Karp, the one that finds no augmenting path, explores everything reachable from free left vertices by alternating paths. Call that reachable set Z. Then the cover is the left vertices not in Z plus the right vertices that are in Z.

def min_vertex_cover(n_left, adj, match_l, dist):
    z_left = {u for u in range(n_left) if dist[u] != INF}      # reached by final BFS
    z_right = {v for u in z_left for v in adj[u]}
    cover_left = [u for u in range(n_left) if u not in z_left]
    return cover_left, sorted(z_right)                          # total size == matching size

Change the example by deleting edge d-4. Phase 2's BFS from b now reaches 1, a, 2, c, 3 and d and stops without a free right vertex. Z is {b, a, c, d} on the left and {1, 2, 3} on the right, so the cover is {1, 2, 3}: three vertices, matching the size-3 maximum matching. That is a proof, checkable by anyone, that no assignment can place all four. Operationally it names the bottleneck: tasks 1 to 3 are the only ones any of these workers can take. The König vertex cover article covers the theorem itself in depth.

Engineering it

Some practical points matter more than the asymptotics.

• Warm starts. The algorithm works from any valid initial matching. A greedy pass that matches each left vertex to its first free neighbour typically finds most of the final matching in O(E), and when the graph changes slightly, starting from the previous matching (minus edges that disappeared) usually needs only a phase or two.
• Data layout. In a compiled language, store adjacency in compressed sparse row form (one offsets array and one targets array) and use integer arrays for dist, match_l and match_r. Pure Python handles graphs with a few hundred thousand edges comfortably; beyond that, use a library.
• Libraries. NetworkX offers networkx.bipartite.hopcroft_karp_matching(G, top_nodes) and networkx.bipartite.to_vertex_cover; SciPy's scipy.sparse.csgraph.maximum_bipartite_matching takes a sparse biadjacency matrix and is implemented with Hopcroft-Karp. Pass top nodes explicitly to NetworkX, because a disconnected graph has no unique bipartition.
• Equivalent formulations. Hopcroft-Karp is Dinic's max-flow algorithm specialized to a unit-capacity network: source to every left vertex, every right vertex to sink, all capacities 1. If you already have Dinic, it gives the same bound. Plain Edmonds-Karp on that network does not.
• Reductions. Minimum path cover in a DAG equals the vertex count minus a maximum matching in the split graph, so the same code schedules chains of dependent jobs.

Failure modes

What to do next

1. Implement the code above and test it against a brute-force matcher on random graphs with up to eight vertices per side.
2. Add a phase counter and confirm it stays small on your real data.
3. Verify every result with the König cover: the cover must touch every edge and have the same size as the matching.
4. Benchmark against scipy.sparse.csgraph.maximum_bipartite_matching and switch if your graphs are large.
5. Add a greedy initial matching and warm starts if the graph changes incrementally.
6. If costs or preferences enter the problem, move on to weighted assignment or stable matching.