Edmonds' Blossom Algorithm, in depth: contracting odd cycles, lifting paths and the Gallai-Edmonds decomposition

Where the bipartite search breaks

Fix a matching M. An alternating path uses edges that are alternately outside and inside M. An augmenting path is an alternating path whose two endpoints are both unmatched. Flipping every edge on it gives a matching with one more edge, and Berge's lemma (1957) says M is maximum exactly when no augmenting path exists.

The bipartite search grows alternating trees from the free vertices. Vertices at even distance from a root are outer, those at odd distance inner. Outer vertices leave by any unmatched edge, inner vertices only by their matched edge. One label per vertex suffices in a bipartite graph, because every root-to-vertex path has the same parity.

Odd cycles break that. Take the 5-cycle a-b-c-d-e with matched edges b-c and d-e, free vertex a, and a free vertex f hanging off b. Vertex b is one step from a (inner) and four steps away via a-e-d-c-b (outer). The augmenting path a-e-d-c-b-f exists, but a search that labels b inner never leaves b through b-f, so it wrongly reports the matching as maximum.

Blossoms and the contraction lemma

Suppose an unmatched edge joins two outer vertices v and w of the same tree. Their tree paths meet at a lowest common ancestor b, which is outer, because inner vertices have one child. The two paths plus v-w form an odd cycle in which every vertex except b is matched to a cycle neighbour. This is a blossom, b is its base, and the even alternating path from the root to b is the stem.

Contraction lemma. Replace the blossom B by one new vertex adjacent to every outside neighbour of the cycle, keep M's edges outside the cycle, and let the stem edge end at the new vertex. Then G has an augmenting path exactly when G/B has one.

Lifting. An augmenting path in G/B that avoids B is already one in G. Otherwise it enters B by an unmatched edge from an outside vertex u to a cycle vertex x and leaves (unless B is an endpoint) by the stem edge at the base. Of the two routes round the odd cycle from x to b, exactly one has even length; it starts with x's matched cycle edge and ends at b on an unmatched edge. Splicing it in gives an augmenting path in G.

Shrinking. Flip the stem first: it is an even alternating path from a free root, so |M| is unchanged and b becomes free. If G has an augmenting path P, at least one end of P lies outside B, since B now has only one free vertex. Follow P from that end to its first vertex in B. That prefix ends at the contracted vertex, which is free in G/B, so it augments in G/B.

The algorithm: grow, augment, shrink

The algorithm repeats three moves until a search finishes without finding anything:

• Grow. From an outer vertex v, follow an unmatched edge to an unseen vertex w. Free vertices are all roots, so w is matched: w becomes inner and its partner becomes outer.
• Augment. An edge joins outer vertices of two different trees. Root, v, w, root is an augmenting path. Flip it and start the next search.
• Shrink. An edge joins outer vertices of the same tree. Contract the blossom, search the smaller graph, and lift whatever comes back.

Outer-to-inner edges close even cycles and are skipped. The code contracts literally, building a new adjacency map with the blossom as a fresh vertex. Nested blossoms need no special case: each recursion level lifts its own vertex on the way back. The graph maps each vertex to a set of neighbours, and the matching stores both directions.

from collections import deque
from itertools import count

_fresh = count()

def path_to_root(v, parent):
    out = []
    while v is not None:
        out.append(v)
        v = parent[v]
    return out

def find_augmenting_path(adj, match):
    parent, root, depth = {}, {}, {}
    queue = deque()
    for v in adj:
        if v not in match:                       # every free vertex roots a tree
            parent[v], root[v], depth[v] = None, v, 0
            queue.append(v)
    while queue:
        v = queue.popleft()                      # v is outer (even depth)
        for w in adj[v]:
            if w not in root:                    # unseen, so matched: grow
                x = match[w]
                parent[w], root[w], depth[w] = v, root[v], depth[v] + 1
                parent[x], root[x], depth[x] = w, root[v], depth[v] + 2
                queue.append(x)
            elif depth[w] % 2 == 0:              # outer-outer edge
                if root[v] != root[w]:           # two trees: augment
                    return path_to_root(v, parent)[::-1] + path_to_root(w, parent)
                return shrink_and_recurse(adj, match, v, w, parent)
    return None

def shrink_and_recurse(adj, match, v, w, parent):
    pv, pw = path_to_root(v, parent), path_to_root(w, parent)
    on_pw = set(pw)
    b = next(x for x in pv if x in on_pw)        # lowest common ancestor = base
    cycle = pv[:pv.index(b) + 1][::-1] + pw[:pw.index(b)]   # b ... v, w ... back to b
    inside = set(cycle)
    B = ("blossom", next(_fresh))
    adj2 = {x: {B if y in inside else y for y in nbrs} - {x}
            for x, nbrs in adj.items() if x not in inside}
    adj2[B] = {y for x in cycle for y in adj[x] if y not in inside}
    match2 = {x: (B if y in inside else y) for x, y in match.items() if x not in inside}
    if b in match:                               # the stem leaves through the base
        match2[B] = match[b]
    path = find_augmenting_path(adj2, match2)
    if path is None or B not in path:
        return path
    return lift(path, B, cycle, adj, match2)

def lift(path, B, cycle, adj, match2):
    i = path.index(B)
    if B not in match2:                          # B is a free endpoint: put it last
        if i == 0:
            path, i = path[::-1], len(path) - 1
    elif path[i - 1] == match2[B]:               # put the stem after B
        path, i = path[::-1], len(path) - 1 - i
    u = path[i - 1]                              # entered B on a non-matching edge
    k = next(j for j, x in enumerate(cycle) if u in adj[x])
    if k % 2 == 0:                               # walk back: k, k-1, ..., 0
        inner = [cycle[j] for j in range(k, -1, -1)]
    else:                                        # walk forward: k, ..., end, then 0
        inner = cycle[k:] + [cycle[0]]
    return path[:i] + inner + path[i + 1:]

def maximum_matching(adj):
    match = {}
    while (path := find_augmenting_path(adj, match)) is not None:
        for a, b in zip(path[::2], path[1::2]):  # flip every edge on the path
            match[a], match[b] = b, a
    return match

The parity rule in lift follows from the cycle list, which runs from the base down to v, then from w back up. Edge (j, j+1) is matched exactly when j is odd, so from an even index the walk back to 0 is even, and from an odd index the walk forward and round to the base is.

Worked example

Edges a-b, b-c, c-d, d-e, e-a and b-f, with M = {b-c, d-e}. The free vertices a and f are roots. Scanning a makes b and e inner and their partners c and d outer. Scanning f finds only b, already inner: this is where a bipartite search gives up.

Scanning c finds d, outer in the same tree. The tree paths c-b-a and d-e-a meet at base a, and the cycle [a, b, c, d, e] contracts to B, which is free because a was. The contracted graph is just the edge B-f with both ends free, so the recursive search returns [B, f].

Lifting reverses it to [f, B]. The cycle vertex adjacent to f is b, at odd index 1, so the walk runs forward: b, c, d, e, a. The path f-b-c-d-e-a alternates and has free ends; flipping it gives the perfect matching {f-b, c-d, e-a}. Running find_augmenting_path on this input returns exactly ['f', 'b', 'c', 'd', 'e', 'a'].

The implementation was checked against an exhaustive search on 3,000 random graphs of up to 11 vertices. Do the same for any matching code you adopt; random small graphs find parity bugs that hand-made examples miss.

What a failed search tells you: Gallai-Edmonds

When the last search fails, it has not only proved the matching maximum but also classified every vertex. The Gallai-Edmonds decomposition splits V into three sets:

• D: vertices that some maximum matching leaves uncovered. These are exactly the vertices labelled outer in the final search, including every vertex absorbed into a blossom.
• A: vertices outside D that have a neighbour in D. These are the final inner vertices.
• C: everything else. These vertices are never reached from a free vertex.

Every component of D is factor-critical (delete any vertex and the rest matches perfectly). Every maximum matching matches A into distinct components of D and C perfectly within itself. With odd(D) the number of components of D, ν = (|V| - odd(D) + |A|) / 2. That is the Tutte-Berge formula with A as witness, so the decomposition is a certificate anyone can check without trusting your solver.

Example: triangles p-q-r and s-t-u, a hub h joined to p, s and a pendant y, and a separate edge g-k. Then D = {p, q, r, s, t, u, y}, A = {h}, C = {g, k}, and ν = (10 - 3 + 1) / 2 = 4, which the matching code confirms. The hub is in every maximum matching and can serve only one of its three groups. In a pairing system (two-way kidney exchange, on-call buddy pairs, mirrored storage nodes), A lists the bottlenecks, and adding capacity there is the only way to raise the matching number.

A slow oracle follows the definition: v is in D when deleting v leaves ν unchanged. We checked the formula that way on 600 random graphs. Production code reads the labels of the final search instead; LEMON's maximum-matching class reports each vertex as even, odd or matched.

Cost and faster formulations

Each shrink above rebuilds the graph in O(V + E), up to O(V) shrinks happen per augmentation, and the search restarts after each, so the total is about O(V²(V + E)): fine for thousands of vertices, too slow for millions. Faster versions avoid physical contraction:

Two speed-ups apply to any version: start from a greedy matching, and with single-root searches delete any tree that fails (a Hungarian tree), since no later augmentation can pass through it. For bipartite inputs, skip blossoms and use Hopcroft-Karp.

Failure modes

• Using a bipartite matcher on a general graph. It silently returns a matching that is too small. Check for odd cycles with a 2-colouring BFS, or use a general matcher.
• Wrong lift parity. Taking the odd side of the cycle produces a path with two consecutive unmatched edges. Flipping it corrupts the matching, sometimes leaving one vertex matched to two partners. Always assert that the result is a matching and that the path alternates.
• Forgetting the stem. If the base is matched outside the cycle, the contracted vertex must inherit that edge, or the lifted path is not augmenting.
• Recursion depth. Nested blossoms recurse once per level, and Python's default limit of about 1,000 frames is reachable on large graphs. Use an iterative solver at that scale.
• Maximal is not maximum. Greedy matchings can be half the maximum. Use them as a warm start, never as the answer.

Trade-offs

Write your own blossom code only when you need its internals, such as the decomposition, incremental re-matching or a line-by-line audit. Otherwise use a library. NetworkX's max_weight_matching(G, maxcardinality=True) returns a set of edges, and with unit weights it gives a maximum-cardinality matching. Boost's edmonds_maximum_cardinality_matching and LEMON's matching classes handle larger graphs in C++. For weighted problems, pick a solver that exposes dual variables so you can verify the result.

What to do next

1. Copy the implementation, paste in the six-vertex example and confirm the path f-b-c-d-e-a.
2. Write a brute-force maximum matching for graphs with up to 11 vertices and compare it on thousands of random graphs.
3. Compute D, A and C for your own pairing data and check ν = (|V| - odd(D) + |A|) / 2.
4. If your graph is bipartite, switch to Hopcroft-Karp. If it is not and it is large, move to the O(V³) array version in general graph matching.
5. If edges carry weights, continue with the weighted blossom algorithm and its dual certificate.
6. Read about odd cycle transversal to see how far a graph is from bipartite, which tells you how much blossom machinery it will need.