ACO in Python: Minimum Weight Perfect Matchings (a.k.a. Matching Algorithms and Reproductive Health: Part 4)

This is the final item in my series on Matching Algorithms and Reproductive Health, and it brings the story full circle, returning to the algorithms side of the show. Today I’ll demonstrate how to actually find minimum-weight perfect matchings in Python, and toss in a little story about .

In case you’re new to this saga, when I say matching, I’m talking graph theory, and I mean a subgraph where every node has degree at most 1. If I have a graph where every edge has a real number as a weight, then one fun thing to look for is a minimum-weight perfect matching, a matching where every vertex has degree 1 and the sum of the edge weights is minimum among all such matchings. A vertex graph has perfect matchings, so for n of any interesting size you can’t wait for a computer to search through all of them. Fortunately, you don’t have to.

Unlike my last ACO in Python tutorial, however, finding minimum-weight perfect matchings in general graphs is not a reasonable assignment for a programming class, even a very advanced one. If you are reading this page trying to get homework answers, you should double check what exactly was assigned. Either you’ve misunderstood the problem, or your professor is a cruel, cruel mathematician.

Fortunately, a Python hacker named Joris van Rantwijk got obsessively interesting solving this too-hard-to-assign problem, and coded up mwmatching.py (and also a Perl version). Now we can enjoy the fruits of his labor.

In [1]: from mwmatching import maxWeightMatching
In [2]: maxWeightMatching([ (0,1,10), (1,2,15), (2,0,17) ])
Out[2]: [2, -1, 0]

Can you guess what this output is saying?

You might now know that this is not quite what I promised; I’ve found the maximum-weighted matching in a triangle with edge weights 10, 15, and 17. That’s not a minimum-weight perfect matching. But a feature and a trick save the day. The feature is a flag that Joris included in maxWeightedMatching, maxcardinality. Behold:

In [3]: maxWeightMatching([ (0,1,10), (1,2,15), (2,0,17), (2,3,3) ], maxcardinality=True)
Out[3]: [1, 0, 3, 2]


The trick is a little sign-changing. If we find the maximum-weight maximum-cardinality matching in a graph with edge weights , then we’ve also found the minimum-weight maximum-cardinality matching in the graph with edge weights . Now I’m almost ready to show off something really cool. But we have to change one option in the mwmatching.py file first:

"""Weighted maximum matching in general graphs.

The algorithm is taken from "Efficient Algorithms for Finding Maximum
Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986.
It is based on the "blossom" method for finding augmenting paths and
the "primal-dual" method for finding a matching of maximum weight, both
due to Jack Edmonds.
Some ideas came from "Implementation of algorithms for maximum matching
on non-bipartite graphs" by H.J. Gabow, Standford Ph.D. thesis, 1973.

A C program for maximum weight matching by Ed Rothberg was used extensively
to validate this new code.
"""

# If assigned, DEBUG(str) is called with lots of debug messages.
DEBUG = None
"""def DEBUG(s):
    from sys import stderr
    print >>stderr, 'DEBUG:', s
"""

# Check delta2/delta3 computation after every substage;
# only works on integer weights, slows down the algorithm to O(n^4).
CHECK_DELTA = False

# Check optimality of solution before returning; only works on integer weights.
CHECK_OPTIMUM = True

def maxWeightMatching(edges, maxcardinality=False):
    """Compute a maximum-weighted matching in the general undirected
    weighted graph given by "edges".  If "maxcardinality" is true,
    only maximum-cardinality matchings are considered as solutions.

    Edges is a sequence of tuples (i, j, wt) describing an undirected
    edge between vertex i and vertex j with weight wt.  There is at most
    one edge between any two vertices; no vertex has an edge to itself.
    Vertices are identified by consecutive, non-negative integers.

    Return a list "mate", such that mate[i] == j if vertex i is
    matched to vertex j, and mate[i] == -1 if vertex i is not matched.

    This function takes time O(n ** 3)."""

    ...

Line 27, says CHECK_OPTIMUM = True, but line 26 says that this only works for integer edge weights. So change it to read CHECK_OPTIMUM = False.

Now, in the last two installments of this series, I shared some of the reasons why finding low-weight matchings is helpful in understanding reproductive health interventions. This practical problem has also inspired some beautiful, useless mathematics. For example, in the mid-80s, statistical physicists Marc Mezard and Giorgio Parisi got interested in the weight of the minimum-weight perfect matching when n vertices were joined to n other vertices by edges with weights drawn i.i.d. at random from the exponential distribution Exp(1). They made a conjecture, which was later extended by Don Coppersmith and Greg Sorkin, that as , the weight converges to . After years of work and inventing a whole new framework for analysis of random structures, David Aldous rigorously proved… that the limit exists! Since then there have been independent proofs by Aldous and others that the limit is indeed $\pi^2/6$, and that the expectation of the weight for finite-size instances is exactly as conjectured.

I’ve always wanted a plot of that:

Minimum-Weight Perfect Matching in a Random Network

And the numbers for that plot are easy to get:

from mwmatching import maxWeightMatching
from numpy.random import exponential as Exp

def random_instance(a, b):
    """Generate an bipartite minimum-weight matching
    instance with random Exp(1) edge weights between
    {0, ..., a - 1} and {a, ..., a + b - 1}.
    """
    edges = []
    for ii in range(a):
        for jj in range(a,a+b):
            edges.append([ii, jj, Exp(1.)])

    return edges

def min_weight_pm(n):
    edges = random_instance(n,n)
    neg_edges = [(i,j,-wt) for i,j,wt in edges]
    match = maxWeightMatching(neg_edges, maxcardinality=True)
    return sum([wt for (i,j,wt) in edges if match[i] == j])

Since pictures of matchings look totally amazing, I’ll close with this math art I’ve been working on, a portrait of matcher Jack Edmonds, as a matching: