Blossom Algorithm, in depth: weighted matching, odd-set duals, the primal-dual loop and certificates you can check

Why the bipartite LP fails on odd cycles

Write x_e = 1 when edge e is in the matching. A matching is a 0/1 vector where every vertex has at most one chosen edge, and we want to maximise the total weight. The obvious relaxation lets x take fractional values with the same degree constraints. On bipartite graphs that relaxation is exact: every vertex of the polytope is integral, which is why the Hungarian algorithm can work purely with vertex potentials.

On general graphs it is not exact, and a triangle shows why. Give the three edges of a triangle weight 1. Any real matching takes one edge, weight 1. The relaxation can set every x_e = 1/2: each vertex sees two half-edges, so the degree constraints hold, and the value is 1.5. The fractional point lives on the odd cycle, exactly where augmenting-path search also breaks.

Edmonds fixed both problems with one family of constraints. For every odd set B of vertices, a matching can use at most (|B| - 1)/2 edges with both ends inside B, because one vertex of B is always left over. Adding those odd-set constraints gives the matching polytope, and Edmonds proved that its vertices are exactly the matchings. On the triangle the constraint for B = {a, b, c} caps the inside edges at 1, and the fractional point is gone. There are exponentially many odd sets, but the algorithm only prices the ones it has shrunk as blossoms.

Duals and the optimality certificate

The dual of that LP has a variable y_v for every vertex and z_B for every odd set. Every edge must be covered: y_u + y_v plus the z_B of every set containing both ends must be at least w(u, v). The dual objective is the sum of all y plus the sum of z_B(|B| - 1)/2, and by weak duality it bounds every matching from above. Call the difference between an edge's cover and its weight its slack; an edge with zero slack is tight.

A matching M and a dual (y, z) prove each other optimal when four complementary slackness conditions hold. Every matched edge is tight. Every vertex with y_v > 0 is matched. Every set with z_B > 0 is full, holding exactly (|B| - 1)/2 matched edges. And every value is non-negative with every edge covered. When all four hold, the matching weight equals the dual objective, so neither can be improved. The triangle's certificate is y = 0 everywhere and z = 1 on {a, b, c}: each edge is covered exactly, the set is full with one matched edge, and the dual value is 1 x (3 - 1)/2 = 1.

The primal-dual loop

The algorithm keeps the dual feasible at all times and grows the matching only through tight edges. Initialise y_v to half the largest edge weight (so every edge is covered), z to zero and the matching to empty. Then repeat stages until no free vertex can be improved:

stage:
    label every free vertex OUTER (a tree root); everything else unlabelled
    loop:
        for each tight edge (u, v) with u OUTER:
            v unlabelled            -> v INNER, v's mate OUTER   (grow the tree)
            v OUTER, other tree     -> augment along root..u-v..root; end stage
            v OUTER, same tree      -> shrink the odd cycle into blossom B (OUTER)
        if nothing tight remains:
            delta = min(delta1, delta2, delta3, delta4)   # see table
            OUTER vertices: y -= delta     INNER vertices: y += delta
            OUTER blossoms: z += 2*delta   INNER blossoms: z -= 2*delta
            if delta == delta1: stop       # a free vertex reached y = 0
            if delta == delta4: expand the INNER blossom whose z hit 0
    expand every OUTER blossom whose z is 0, then start the next stage

The dual step is the heart of it. Tree edges join an outer to an inner vertex, so lowering one end and raising the other keeps them tight, while edges leaving outer vertices lose slack until a new one becomes tight. The z updates keep edges inside a blossom tight. There are O(V) stages, and careful bookkeeping gives O(V³) total, which is the bound most libraries implement.

The four dual adjustments

delta1 exists only when leaving a vertex unmatched is allowed; in the perfect-matching variant the y are free and stages end only by augmentation. The halving in delta3 and delta4 is why implementations use doubled integer weights: every dual stays a multiple of one half, and doubling keeps comparisons exact.

Worked example: greedy 13, optimum 17

Take the six-vertex graph in the figure: a triangle a-b-c with weights 6, 7 and 6, and a tail c-d (5), d-e (4), e-f (6). Greedy by weight picks b-c (7), then e-f (6), and then nothing else fits: a's only neighbours b and c are taken, and d's neighbours c and e are too. Total 13, with two vertices stranded.

The optimum is a-b + c-d + e-f = 17. It gives up the heaviest edge, which is exactly the kind of non-local trade a primal-dual method makes and a greedy one cannot. You do not have to trust that number. The prices y = (a 3, b 3, c 4, d 1, e 3, f 3) cover every edge: a-b 6, b-c 7, c-a 7 against weight 6, c-d 5, d-e 4, e-f 6. The three matched edges are tight, every vertex with a positive price is matched, and the prices sum to 17. Weak duality says no matching beats 17, and we have one with 17, so it is optimal. Notice that this certificate needs no blossom price; the triangle's odd-set constraint is not binding here because the tail gives c a better partner.

Calling and testing a library

Unless you are building a solver, call a tested implementation. In Python, NetworkX ships max_weight_matching, an O(V³) weighted blossom, and min_weight_matching. Read the semantics carefully: maxcardinality=True returns the heaviest matching among those of maximum size, and in NetworkX 3.7 min_weight_matching returns the lightest matching among maximum-cardinality matchings, which on a graph with a perfect matching is minimum-weight perfect matching. Older releases documented different behaviour, so pin the version and test it.

import itertools, random
import networkx as nx

def brute_max_weight(G):
    edges = list(G.edges(data="weight"))
    best = 0
    for r in range(len(edges) + 1):
        for combo in itertools.combinations(edges, r):
            ends = [x for u, v, _ in combo for x in (u, v)]
            if len(ends) == len(set(ends)):
                best = max(best, sum(w for _, _, w in combo))
    return best

rng = random.Random(7)
for trial in range(1000):
    G = nx.gnp_random_graph(rng.randint(2, 8), 0.5, seed=rng.randint(0, 10**9))
    for u, v in G.edges:
        G[u][v]["weight"] = rng.randint(1, 20)
    M = nx.max_weight_matching(G)
    assert nx.is_matching(G, M)
    assert sum(G[u][v]["weight"] for u, v in M) == brute_max_weight(G), trial
print("1000 random graphs agree")

That property test passed against NetworkX 3.7 when this article was written, and on the worked example all three calls return {a-b, c-d, e-f} with weight 17. Keep a test like this in your repository: it is cheap, and it catches the realistic bugs, which are in graph construction and weight conventions, not in the library.

Two cases need care. For maximum weight at maximum cardinality, pass maxcardinality=True; on the path a-b (1), b-c (10), c-d (1) the plain call returns b-c alone (10) while the cardinality-first call returns a-b + c-d (2). For minimum-weight perfect matching on a graph that might not have one, check the returned size equals V/2 before using the answer; otherwise you silently get a smaller matching.

Checking any answer with its duals

If you write your own solver, or port one to C++ or a GPU, make it emit its final duals and verify them. The check is linear in the input and turns a heavy algorithm into something you can trust in production:

def check_certificate(edges, matching, y, z):
    """edges {(u, v): w}; matching set of (u, v); y {v: price}; z {frozenset: price}."""
    matched = {v for e in matching for v in e}
    assert all(val >= 0 for val in y.values()) and all(val >= 0 for val in z.values())
    for (u, v), w in edges.items():
        cover = y[u] + y[v] + sum(zb for B, zb in z.items() if u in B and v in B)
        assert cover >= w, f"edge {u}-{v} under-covered"
        if (u, v) in matching or (v, u) in matching:
            assert cover == w, f"matched edge {u}-{v} not tight"
    for v, val in y.items():
        assert val == 0 or v in matched, f"{v} has positive dual but is free"
    for B, zb in z.items():
        inside = sum(1 for u, v in matching if u in B and v in B)
        assert zb == 0 or inside == (len(B) - 1) // 2, f"blossom {set(B)} not full"
    primal = sum(edges.get((u, v), edges.get((v, u))) for u, v in matching)
    dual = sum(y.values()) + sum(zb * ((len(B) - 1) // 2) for B, zb in z.items())
    assert primal == dual
    return primal

On the triangle with y = 0 and z = 1 on {a, b, c} it returns 1; on the worked example with the prices above it returns 17; given greedy's matching it fails with "a has positive dual but is free".

Faster solvers and approximations

O(V³) is fine to tens of thousands of vertices on dense graphs. Beyond that, the choice of implementation matters more than the theory.

• Gabow (1990) reaches O(V(E + V log V)) with careful priority queues, which matters on sparse graphs.
• Blossom V (Kolmogorov, 2009) is the standard C++ code for minimum-weight perfect matching, with greedy initialisation and per-tree dual updates; it is often far faster than textbook code. Check its licence before shipping it.
• Sparse blossom (Higgott and Gidney, 2023), used by PyMatching 2, grows regions from defects instead of building a dense complete graph, which is what lets surface-code decoding keep up in real time.
• Approximation: greedy is a 1/2-approximation, and Duan and Pettie gave a (1 - eps) approximation in near-linear time. If a 2 percent loss is acceptable, these are much simpler to run at scale. See approximation algorithms for the general trade.

Failure modes

• Floating-point weights. Tightness tests compare slacks with zero. With floats, an edge can be 1e-12 short of tight and the search stalls or returns a wrong matching. Scale to integers before solving.
• Wrong problem variant. Maximum-weight, maximum-weight at maximum cardinality and minimum-weight perfect are different problems; calling the first when you need the third returns a valid, smaller matching with no error.
• Missing perfect matching. Odd vertex counts, or a Tutte-Berge obstruction, mean no perfect matching exists. Check the size of the result, or add dummy vertices with explicit costs if leaving a vertex unmatched has a real price.
• Dense reductions. Christofides and postman reductions build a complete graph on the odd-degree vertices: over a billion edges at 50,000 of them. Use a sparse solver.
• Expansion bugs. Hand-written solvers most often fail when an inner blossom is expanded and its vertices are relabelled. Fuzz against brute force on graphs of 8 to 12 vertices with many triangles.

Trade-offs

The weighted blossom algorithm is exact and certifiable; its costs are implementation complexity and cubic worst-case time. A MIP solver with lazily added odd-set cuts is easier to extend with side constraints but slower. Greedy is fast and guarantees only one half. In production, use a tested library for the exact problem, a check in CI, and an approximation only where the measured loss is acceptable. For a full application, the Chinese postman problem shows minimum-weight perfect matching doing the heavy lifting.

What to do next

1. Write down which variant you need: maximum weight, maximum weight at maximum cardinality, or minimum-weight perfect.
2. Convert weights to integers and confirm the graph is not bipartite before choosing a general solver.
3. Run the brute-force property test above against the exact library version you deploy.
4. On the worked example, confirm your solver returns 17, not greedy's 13.
5. If you write or port a solver, emit duals and run the certificate checker on every output in tests.
6. For more than about 10,000 vertices, benchmark a sparse implementation such as Blossom V or sparse blossom before scaling hardware.