Matroid Intersection, in depth: exchange graphs, shortest augmenting paths, the rank certificate and weighted versions

Two constraints at once

A matroid M on ground set E is a family of independent sets that contains the empty set, is closed under taking subsets, and satisfies the augmentation rule: if A and B are independent and A is smaller, some element of B can be added to A. The rank r(S) is the size of the largest independent subset of S. In code you never store the family. You call an independence oracle, a function that takes a set and answers yes or no.

Given two matroids M1 and M2 on the same ground set, a common independent set is independent in both. The task is to find one of maximum size, or of maximum weight. Three standard instances show the range:

Bipartite matching is the familiar special case. If you know augmenting paths for matching, the matroid algorithm is the same idea with the oracle deciding which swaps are legal. For three matroids the problem becomes NP-hard, because Hamiltonian path in a directed graph is the intersection of a graphic matroid and two partition matroids. Two is the boundary.

The exchange graph

Suppose I is common independent. We want to grow it by one, but adding a single element may break one of the matroids. The fix is a chain of swaps. Build a directed graph on the ground set:

• Sources X1: elements x not in I with I + x independent in M1.
• Sinks X2: elements x not in I with I + x independent in M2.
• An arc from y in I to x not in I when I - y + x is independent in M1.
• An arc from x not in I to y in I when I - y + x is independent in M2.

A path x0, y1, x1, y2, ..., xk from a source to a sink alternates between elements outside and inside I. Flipping membership along it adds k + 1 elements and removes k, so the set grows by one. Each M1 arc says one swap is safe in M1, and each M2 arc says one swap is safe in M2. The subtle point is that individually safe swaps are not always jointly safe. The guarantee holds only when the path is shortest, with no shortcuts between its vertices. A breadth-first search from the sources returns a shortest path, so use BFS and never an arbitrary DFS path.

If an element is in both X1 and X2, the path has one vertex and you simply add it. If no source can reach a sink, I is already maximum. That claim is the min-max theorem below, and it is what makes the algorithm trustworthy.

The algorithm in Python

The implementation below is generic: pass the ground set and two oracle functions. It returns the common independent set and the set R of elements reachable from the sources in the final exchange graph, which is the optimality certificate.

from collections import deque

def matroid_intersection(ground, indep1, indep2, start=()):
    """Largest common independent set. indep1/indep2 take a set, return bool."""
    I = set(start)
    while True:
        out = [x for x in ground if x not in I]
        sources = {x for x in out if indep1(I | {x})}
        sinks = {x for x in out if indep2(I | {x})}
        arcs = {v: [] for v in ground}
        for y in I:
            for x in out:
                J = (I - {y}) | {x}
                if indep1(J):
                    arcs[y].append(x)      # swap y for x keeps M1 happy
                if indep2(J):
                    arcs[x].append(y)      # swap y for x keeps M2 happy
        parent = {x: None for x in sources}
        queue, end = deque(sources), None
        while queue:                       # BFS gives a shortest path
            v = queue.popleft()
            if v in sinks:
                end = v
                break
            for w in arcs[v]:
                if w not in parent:
                    parent[w] = v
                    queue.append(w)
        if end is None:
            return I, set(parent)          # R: everything reachable from sources
        while end is not None:             # flip membership along the path
            I ^= {end}
            end = parent[end]

Rebuilding the graph every round is deliberate: after an augmentation any arc can change, so caching arcs across rounds is a classic source of wrong answers.

Worked example: a rainbow spanning tree

Take a graph on vertices A, B, C, D with five coloured edges: a = A-B red, b = B-C blue, c = A-C green, d = C-D red, e = B-D blue. We want a spanning tree, three edges, with no repeated colour. Greedy in the order a, b, c, d, e takes a and b. Then c closes the triangle A-B-C, d repeats red, and e repeats blue. Greedy stops at two edges, and D is left uncovered.

Now build the exchange graph for I = {a, b}. The sources are d and e: either joins the forest without a cycle, but each repeats a colour. The only sink is c, the unused green, which would close a cycle. Every swap of a or b for an outside edge keeps a forest, so all six M1 arcs exist. The colour arcs exist only where colours match up: c can replace a or b, d can replace a (both red), and e can replace b (both blue).

BFS from {d, e} reaches a from d and b from e, then reaches the sink c. Take the path d, a, c: add d, remove a, add c. The result {b, c, d} is B-C, C-D and A-C, a star at C that spans all four vertices in blue, red and green. The other shortest path, e, b, c, gives {a, c, e}, which is also valid. Which one you get depends only on BFS order. Running the code with start={'a', 'b'} reproduces this, and a further round finds no source, so it stops at size three.

Why it is optimal: the rank certificate

Edmonds' matroid intersection theorem states that the largest common independent set has size equal to the minimum, over all subsets A of E, of r1(A) + r2(E - A). One direction is easy. Any common independent I splits into the part inside A, which is independent in M1, and the part outside, which is independent in M2, so |I| is at most r1(A) + r2(E - A) for every A.

The algorithm supplies the other direction. When no path exists, let R be the set reachable from the sources. Then |I| = r1(E - R) + r2(R). This gives you a checkable proof of optimality: compute two ranks and compare. For bipartite matching it reduces to Konig's theorem that maximum matching equals minimum vertex cover. That is the same min-max shape as LP duality, and in fact the matroid intersection polytope is described by exactly these rank inequalities.

In production, compute the certificate on a sample of runs, or on every run when the oracle is cheap, and alert when it fails. A mismatch means a buggy oracle, not a bad algorithm, and it is the cheapest bug detector you will get.

Weighted matroid intersection

With weights w, you want the maximum-weight common independent set. Give every exchange-graph vertex a length: -w(x) for x outside I and +w(y) for y inside I. Adding x gains w(x) and removing y loses w(y), so a path of minimum total length is the most profitable augmentation. Among minimum-length paths, take one with the fewest arcs. Lengths can be negative, so use Bellman-Ford rather than BFS. The algorithm maintains a strong invariant: after k augmentations, I is a maximum-weight common independent set among those of size k.

best, I = set(), set()
for k in range(1, rank_bound + 1):
    D = exchange_graph(I)                 # same arcs as before
    P = bellman_ford_path(D, sources, sinks,
                          length=lambda v: w[v] if v in I else -w[v],
                          tie_break="fewest_arcs")
    if P is None:
        break                             # no common independent set of size k
    I = I ^ set(P)
    if weight(I) > weight(best):
        best = set(I)                     # keep the best over all sizes

For a minimum-cost common basis, such as the cheapest rainbow spanning tree, set w to minus the cost and keep the iterate of full size rather than the best over all sizes. Frank's weight-splitting algorithm from 1981 is an alternative that maintains dual weights instead of searching with negative lengths.

Cost and oracle engineering

Cost is measured in oracle calls, because the oracle usually dominates. A round tests |I| times (n - |I|) swaps against each matroid, which is at most about 2rn calls, where r is the answer size. There are at most r rounds, so the simple algorithm makes O(r squared n) oracle calls plus BFS work. Cunningham showed in 1986 that augmenting in phases of equal shortest-path length, as Hopcroft-Karp does for matching, lowers this to O(n r to the 1.5) calls. That is why Hopcroft-Karp is the right mental model.

The incremental forms turn each round into roughly linear work per element of I, and they are where most of the speed in a real implementation comes from.

Operational guidance

• Prefer a specialised solver when one exists. If both matroids are partition matroids, you have bipartite matching or flow. Use Hopcroft-Karp or a max-flow library, which will be much faster than a generic oracle loop.
• Test oracles alone first. Check subset-closure and augmentation on random small sets by brute force. A non-matroid oracle silently breaks the optimality guarantee.
• Warm-start. Seed start with a greedy solution. Greedy is usually close, so only a few augmentations remain.
• Use exact arithmetic for linear matroids. Floating-point rank tests flip with tolerance. Use rationals, or a finite field if the application allows it.

Failure modes

• Taking any path, not a shortest one. A path with a shortcut can produce a set that is dependent in M1 or M2. Always use BFS, or the fewest-arcs tie-break in the weighted case.
• Arc directions swapped. M1 arcs go from inside I to outside, and M2 arcs go from outside to inside. Reversing one family gives plausible output that is not maximum. The rank certificate catches it.
• A non-matroid constraint. Budgets, knapsack limits and "no two adjacent" rules are not matroids. The code still runs and returns a common independent set, but it is not optimal.
• Three constraints. Adding a third matroid makes the problem NP-hard. Either fold one constraint into another matroid, use an approximation, or move to integer programming.
• Stale caches. Reusing arcs or oracle state from an earlier I after a flip.

Trade-offs

What to do next

1. Run the code on the five-edge example with start={'a', 'b'} and check that it flips a three-vertex path.
2. Write a brute-force rank function and assert |I| = r1(E - R) + r2(R) on 100 random small instances.
3. Model one of your own two-constraint selection problems as two oracles, and test each oracle for the matroid axioms before trusting the result.
4. If both oracles are partition matroids, switch to a matching solver and compare answers.
5. Replace the naive oracles with the incremental forms in the table and measure oracle calls per round.
6. Read the matroid greedy theorem in practice to see when a single constraint lets greedy finish the job.