How to get different maximum matchings

I have a large bipartite graph and I can find a maximum matching quickly using Hopcroft-Karp. But what I would really like is a few hundred different maximum matchings for the same graph. How can I get those?

Here is an example small bipartite graph with a maximum matching shown.

Similar graphs can be made with:

import igraph as ig
from scipy.sparse import random, find
from scipy import stats
from numpy.random import default_rng
import numpy as np
from igraph import Graph, plot
np.random.seed(7)
rng = default_rng()
rvs = stats.poisson(2).rvs
S = random(20, 20, density=0.35, random_state=rng, data_rvs=rvs)
triples = [*zip(*find(S))]
edges = [(triple[0], triple[1]+20) for triple in triples]
print(edges)
types = [0]*20+[1]*20
g = Graph.Bipartite(types, edges)
matching = g.maximum_bipartite_matching()
layout = g.layout_bipartite()
visual_style={}
visual_style["vertex_size"] = 10
visual_style["bbox"] = (600,300)
plot(g, bbox=(600, 200), layout=g.layout_bipartite(), vertex_size=20, vertex_label=range(g.vcount()), 
    vertex_color="lightblue", edge_width=[5 if e.target == matching.match_of(e.source) else 1.0 for e in g.es], edge_color=["red" if e.target == matching.match_of(e.source) else "black" for e in g.es]
    )

There’s an enumeration algorithm due to Fukuda and Matsui (“Finding all the perfect matchings in bipartite graphs”), (pdf) which was improved for non-sparse graphs by Uno (“Algorithms for enumerating all perfect, maximum and maximal matchings in bipartite graphs”) at the cost of more implementation complexity.

Given the graph G, we find a matching M (e.g., with Hopcroft–Karp) to pass along with G to the root of a recursive enumeration procedure. On input (G, M), if M is empty, then the procedure yields M. Otherwise, the procedure chooses an arbitrary e ∈ M. A maximum matching in G either contains e or not. To enumerate the matchings that contain e, delete e’s endpoints from G to obtain G′, delete e from M to obtain M′, make a recursive call for (G′, M′), and add e to all of the matchings returned. To enumerate the matchings that don’t contain e, delete e from G to obtain G′′ and look for an augmenting path with respect to (G′′, M′). If we find a new maximum matching M′′ thereby, recur on (G′′, M′′).

With Python you can implement this procedure using generators and then grab as many matchings as you like.

def augment_bipartite_matching(g, m, u_cover=None, v_cover=None):
    level = set(g)
    level.difference_update(m.values())
    u_parent = {u: None for u in level}
    v_parent = {}
    while level:
        next_level = set()
        for u in level:
            for v in g[u]:
                if v in v_parent:
                    continue
                v_parent[v] = u
                if v not in m:
                    while v is not None:
                        u = v_parent[v]
                        m[v] = u
                        v = u_parent[u]
                    return True
                if m[v] not in u_parent:
                    u_parent[m[v]] = v
                    next_level.add(m[v])
        level = next_level
    if u_cover is not None:
        u_cover.update(g)
        u_cover.difference_update(u_parent)
    if v_cover is not None:
        v_cover.update(v_parent)
    return False


def max_bipartite_matching_and_min_vertex_cover(g):
    m = {}
    u_cover = set()
    v_cover = set()
    while augment_bipartite_matching(g, m, u_cover, v_cover):
        pass
    return m, u_cover, v_cover


def max_bipartite_matchings(g, m):
    if not m:
        yield {}
        return
    m_prime = m.copy()
    v, u = m_prime.popitem()
    g_prime = {w: g[w] - {v} for w in g if w != u}
    for m in max_bipartite_matchings(g_prime, m_prime):
        assert v not in m
        m[v] = u
        yield m
    g_prime_prime = {w: g[w] - {v} if w == u else g[w] for w in g}
    if augment_bipartite_matching(g_prime_prime, m_prime):
        yield from max_bipartite_matchings(g_prime_prime, m_prime)


# Test code

import itertools
import random


def erdos_renyi_random_bipartite_graph(n_u, n_v, p):
    return {u: {v for v in range(n_v) if random.random() < p} for u in range(n_u)}


def is_bipartite_matching(g, m):
    for v, u in m.items():
        if u not in g or v not in g[u]:
            return False
    return len(set(m.values())) == len(m)


def is_bipartite_vertex_cover(g, u_cover, v_cover):
    for u in g:
        if u in u_cover:
            continue
        for v in g[u]:
            if v not in v_cover:
                return False
    return True


def is_max_bipartite_matching(g, m, u_cover, v_cover):
    return (
        is_bipartite_matching(g, m)
        and is_bipartite_vertex_cover(g, u_cover, v_cover)
        and len(m) == len(u_cover) + len(v_cover)
    )


def brute_force_count_bipartite_matchings(g, k):
    g_edges = [(v, u) for u in g for v in g[u]]
    count = 0
    for m_edges in itertools.combinations(g_edges, k):
        m = dict(m_edges)
        if len(m) == k and is_bipartite_matching(g, m):
            count += 1
    return count


def test():
    g = erdos_renyi_random_bipartite_graph(7, 7, 0.35)
    m, u_cover, v_cover = max_bipartite_matching_and_min_vertex_cover(g)
    assert is_max_bipartite_matching(g, m, u_cover, v_cover)

    count = 0
    for m_prime in max_bipartite_matchings(g, m):
        assert is_bipartite_matching(g, m_prime)
        assert len(m_prime) == len(m)
        count += 1
    assert brute_force_count_bipartite_matchings(g, len(m)) == count


for i in range(100):
    test()