Hi everyone!

Some time ago Monogon wrote an article about Edmonds blossom algorithm to find the maximum matching in an arbitrary graph. Since the problem has a very nice algebraic approach, I wanted to share it as well. I'll start with something very short and elaborate later on.

tl;dr. The Tutte matrix of the graph $$$G=(V, E)$$$ is

\begin{equation} T(x_{12}, \dots, x_{(n-1)n}) = \begin{pmatrix} 0 & x_{12} e_{12} & x_{13} e_{13} & \dots & x_{1n} e_{1n} \newline -x_{12} e_{12} & 0 & x_{23} e_{23} & \dots & x_{2n} e_{2n} \newline -x_{13} e_{13} & -x_{23} e_{23} & 0 & \dots & x_{3n} e_{3n} \newline \vdots & \vdots & \vdots & \ddots & \vdots \newline -x_{1n} e_{1n} & -x_{2n} e_{2n} & -x_{3n} e_{3n} & \dots & 0 \end{pmatrix} \end{equation}

Here $$$e_{ij}=1$$$ if $$$(i,j)\in E$$$ and $$$e_{ij}=0$$$ otherwise, $$$x_{ij}$$$ are formal variables. Key facts:

1. Graph has a perfect matching if and only if $$$\det T \neq 0$$$ when considered as polynomial of $$$x_{ij}$$$.
2. Rank of $$$T$$$ is the number of vertices in the maximum matching.
3. Maximal linearly independent subset of rows corresponds to the subset of vertices on which it is a perfect matching.
4. If graph has a perfect matching, $$$(T^{-1})_{ij} \neq 0$$$ iff there exists a perfect matching which includes the edge $$$(i,j)$$$.
5. After such $$$(i,j)$$$ is found, to fix it in the matching one can eliminate $$$i$$$-th and $$$j$$$-th rows and columns of $$$T^{-1}$$$ and find next edge.

Randomization comes when we substitute $$$x_{ij}$$$ with random values. It can be proven that conditions above still hold with high probability. This provides us with $$$O(n^3)$$$ algorithm to find maximum matching in general graph. For details, dive below.

Let's start with bipartite case. For bipartite graph Edmonds matrix $$$E$$$ is defined as $$$E_{ij}=1$$$ if there is an edge between $$$i$$$-th vertex of the first part and $$$j$$$-th vertex of the second part and $$$E_{ij}=0$$$ otherwise. Let $$$S_n$$$ be a set of all permutations of $$$n$$$ elements. Then

\begin{equation} \begin{matrix} \text{per }E=\sum\limits_{\sigma \in S_n} \prod\limits_{i=1}^n E_{i \sigma_i}, \newline \det E=\sum\limits_{\sigma \in S_n} \text{sgn }\sigma \prod\limits_{i=1}^n E_{i\sigma_i} \end{matrix} \end{equation}

Where $$$\text{per} E$$$ and $$$\det E$$$ are the permanent and determinant of $$$E$$$, and $$$\text{sgn }\sigma$$$ is a sign of permutation $$$\sigma$$$. By definition, $$$\text{per }E$$$ is the number of perfect matchings in $$$G$$$, so to check if there is a perfect matching is equivalent to check if this permanent is non-zero. Calculating permanent is $$$\text{#}P$$$-complete task, but calculating determinant is possible in $$$O(n^3)$$$ with gaussian elimination.

To check if a permanent of some non-negative matrix is zero is far easier task than actually computing it, because it may be probabilistically reduced to the determinant computing. Instead of $$$E$$$ consider matrix $$$E^{(p)}$$$ such that $$$E^{(p)}_{ij} = E_{ij} x_{ij}$$$. For such matrix, \begin{equation} \text{per }E = 0 \iff \det E^{(p)}=0. \end{equation} Indeed, the permanent of non-negative matrix is non-zero iff at least one summand is non-zero and every summand in $$$\det E^{(p)}$$$ has a unique monomial associated with it, thus it is identical zero iff every summand is zero.

If we have a multivariate polynomial $$$P(x_1, \dots, x_k)$$$ of degree $$$n$$$ and we check if it's identical zero by substituting each $$$x_i$$$ with a random value chosen from a set $$$S$$$, the probability of obtaining false positive is, by Schwartz–Zippel lemma, at most $$$\frac{n}{|S|}$$$.

Let's move on to the general case. You may note that Tutte matrix of bipartite graph with $$$n$$$-size parts looks like

\begin{equation} T = \begin{pmatrix} 0 & E \newline -E^T & 0 \end{pmatrix} \end{equation}

Thus, $$$\det T = \det^2 E$$$ for such graph. Now, how should we interpret such determinant for generic skew-symmetric matrixes? To answer this question we should perceive permutation as a set of cycles. In this case, $$$\det T$$$ is some alternating sum over all directed cycle covers of $$$G$$$. Specifically, since $$$T_{ii}=0$$$, only derangements are considered. Note that $$$T_{ij}=-T_{ji}$$$, thus for any cycle $$$(i_1, i_2, \dots, i_k)$$$ it holds that \begin{equation} T_{i_1 i_2} T_{i_2 i_3} \dots T_{i_k i_1}=(-1)^k T_{i_1 i_k} \dots T_{i_3 i_2} T_{i_2 i_1} \end{equation}

That means that if $$$\sigma$$$ contains a cycle of odd length, its summand will be cancelled out by the same permutation with this specific cycle reversed. Other permutations may be joined into classes of equivalence by comparing undirected sets of cycles corresponding to them. Equivalence class of a permutation that is comprised of even cycles, $$$k$$$ of them having length greater than $$$2$$$, will contain $$$2^k$$$ permutations (for every cycle with length greater than $$$2$$$ its reversal will yield a new permutation) and products of $$$T_{1\sigma_1} \dots T_{n\sigma_n}$$$ will be the same for all such permutations due to aforementioned property.

In this way, $$$\det T\neq 0$$$ if and only if there is a cycle cover of $$$G$$$ comprised of even-length cycles. From any such cycle cover we may construct a perfect matching by taking alternating edges in each cycle, and in the same time perfect matching itself is a cycle cover of $$$2$$$-length cycles, thus it is equivalent to graph having a perfect matching.

A small digression about this determinant. One may note that for $$$(i_1, i_2, \dots, i_{2k})$$$ it holds that

\begin{equation} \begin{matrix} T_{i_1 i_2} T_{i_2 i_3} \dots T_{i_{2k} i_1} = (T_{i_1 i_2} T_{i_3 i_4} \dots T_{i_{2k-1}i_{2k}}) (T_{i_2 i_3} T_{i_4 i_5} \dots T_{i_{2k}i_{1}}) \end{matrix} \end{equation}

Thus, any permutation comprised of even-length cycles can be represented as a combination of two perfect matchings. It highlights that $$$\det T = \text{Pf }^2 T$$$ for even $$$n$$$, where $$$\text{Pf }T$$$ is the so-called Pfaffian of skew-symmetric matrix:

\begin{equation} \text{Pf }T = \frac{1}{2^{n/2}\left(\frac{n}{2}\right)!}\sum\limits_{\sigma \in S_n} \text{sgn }\sigma \cdot T_{\sigma_1 \sigma_2} T_{\sigma_3 \sigma_4} \dots T_{\sigma_{n-1} \sigma_n} \end{equation}

Every perfect matching is counted $$$2^{n/2} \left(\frac{n}{2}\right)!$$$ times in the sum, as you may obtain the same matching by changing the order of multipliers or changing $$$T_{ij}$$$ with $$$T_{ji}$$$. This is canceled by dividing the sum with this amount. After two arbitrary summands of Pfaffian are multiplied, they may be perceived as a product over edges in some cycle cover. And every specific cycle cover may be obtained as a product of two perfect matchings in $$$2^k$$$ way, given that the cover has $$$k$$$ cycles of length greater than $$$2$$$.

The determinant of Tutte matrix now allows us to check if a graph has a perfect matching. But what if there is no perfect matching and we want to recover the size of maximum matching? As was already mentioned, this size is equal to the rank $$$\text{rank }T$$$ of Tutte matrix.

Let $$$\text{rank }T=m$$$ and $$${a_1, \dots, a_m}$$$ is the set of linearly independent rows of $$$T$$$. By skew-symmetry, columns of $$$T$$$ with these indices are linearly independent as well, from which we may conclude that submatrix comprised of these rows and columns is non-singular, thus has non-zero determinant.

Indeed, let's assume that it is not true. Then we may nullify the first row of this submatrix with Gaussian elimination. On the other hand, since $$$m$$$ is the rank of matrix, and columns are linearly independent in $$$T$$$, every other column can be obtained as their linear combination, which means that after elimination of submatrix's first row, the whole $$$a_1$$$ row would be nullified, which contradicts the assumption of rows being linearly independent.

On the other hand, if $$$m$$$ is the size of maximum matching, corresponding submatrix will be non-singular, making rank of $$$T$$$ at least $$$m$$$. That being said, one can find a set of vertices on which maximum matching is achieved by finding a maximal set of linearly independent rows in Tutte matrix.

Finally, let's tackle the task of restoring perfect matching. It is known from linear algebra that

\begin{equation} T^{-1} = \frac{\text{adj }T}{\det T} \end{equation}

Where $$$\text{adj }T$$$ is adjugate matrix, that is $$$(\text{adj }T)_{ji}$$$ is $$$(-1)^{i+j}$$$ times the determinant of $$$T$$$ with $$$i$$$-th row and $$$j$$$-th column removed. In terms of sum over all cycle covers, deletion of $$$i$$$-th rown and $$$j$$$-th column is up to constant multiplier equivalent to turning all elements in $$$i$$$-th row and $$$j$$$-th column except for $$$T_{ij}$$$ into zero, effectively preventing $$$i$$$-th vertex to be followed by anything but $$$j$$$ in the cycle cover. Such sum would be non-zero if and only if there exists a cyclic cover of even-length cycles which passes through $$$i \to j$$$ edge.

This immediately gives us $$$O(n^4)$$$ algorithms to find the perfect matching as we may find any edge $$$(i, j)$$$ that is guaranteed to be included in at least one perfect matching and then removing $$$i$$$ and $$$j$$$ vertices from the graph (correspondingly, $$$i$$$-th and $$$j$$$-th rows and columns from $$$T$$$) and repeating the procedure of picking next such edge.

To recover maximum matching, one may firstly recover set of vertices on which it is a perfect matching by finding maximal linearly independent set of matrix rows and then find a perfect matching in a graph reduced to only these vertices.

Marcin Mucha and Piotr Sankowski from Warsaw university suggested $$$O(n^3)$$$ improvement in 2004. They noticed that instead of calculating inverse matrix from scratch every time, one may simply do Gaussian elimination on corresponding rows and columns.

For simplicity, let's say that we want to remove $$$1$$$-st row and $$$1$$$-st column of matrix $$$A$$$. Let

\begin{equation} \begin{matrix} A=\begin{pmatrix} a_{11} & v^T \newline u & B \end{pmatrix} & A^{-1}=\begin{pmatrix} \hat a_{11} & \hat v^T \newline \hat u & \hat B \end{pmatrix} \end{matrix} \end{equation}

What we actually want to do here is to make transition from $$$A^{-1}$$$ to $$$B^{-1}$$$. For this we note

\begin{equation} AA^{-1} = \begin{pmatrix} a_{11} \hat a_{11} + v^T \hat u & a_{11} \hat v^T + v^T \hat B \newline u \hat a_{11} + B \hat u & u \hat v^T + B \hat B \end{pmatrix} = \begin{pmatrix} I_1 & 0 \newline 0 & I_{n-1} \end{pmatrix} = I_n \end{equation}

From these matrix equations we get

\begin{equation} \begin{matrix} B \hat B = I_{n-1} — u \hat v^T,\newline (u \hat a_{11} + B \hat u)\hat v^T = 0 \end{matrix} \end{equation}

Second equation implies $$$B \hat u \hat v^T = - \hat a_{11} u \hat v^T$$$, dividing it by $$$-\hat a_{11}$$$ and adding to the first one, we obtain

\begin{equation} B\left(\hat B -\frac{\hat u \hat v^T}{\hat a_{11}}\right) = I_{n-1} \end{equation}

Thus $$$B^{-1} = \hat B - \hat u \hat v^T / \hat a_{11}$$$. In more generic terms, this expression is called a Schur complement of $$$\hat a_{11}$$$ and it can be obtained by Gaussian elimination of the first row and first column of $$$A^{-1}$$$. For erasure of $$$i$$$-th row and $$$j$$$-th column, one may eliminate these row and column with Gauss algorithm instead of computing inverse matrix from scratch. Ultimately, it results in the following algorithm:

1. Compute $$$B=T^{-1}$$$.
2. Find any $$$(i, j)$$$ such that $$$T_{ij} \neq 0$$$ and $$$B_{ij} \neq 0$$$.
3. Add $$$(i, j)$$$ to the current matching and eliminate $$$i$$$-th and $$$j$$$-th columns and rows of $$$B$$$ with Gaussian elimination.
4. Repeat steps 2-3 until the matching is found.

Due to linear algebra properties, the process will converge for any non-singular matrix. Thus the failure probability is still at most $$$\frac{n}{|S|}$$$ due to Schwartz–Zippel lemma and it may only occur if matrix rank turned out to be less than actual maximum matching size.