Counting Perfect Matchings and Hamiltonian Cycles Faster

In this paper, we further improve the algorithm of Björklund, Kaski, and Williams [11], removing the $\log \log n$ term in the exponent. We also show how to extend the complexity bound for permanent to a natural extension called hafnian.

For a $2n\times 2n$ symmetric matrix $A$ over a commutative ring $R$, the $R$-Hafnian of $A$ is defined as

where $P_{2n}$ is the family of partitions of $[2n]$ into $n$ pairs. The permanent of an $n\times n$ matrix can be reduced to the hafnian of a $2n\times 2n$ matrix via the following basic relation:

Previously the best known algorithm for hafnian ran in time $\tilde{O}(2^n)$ (for a $2n\times 2n$ matrix), first developed by Björklund [5] and alternatively by Cygan and Pilipczuk [15].

The Chinese remainder theorem and a simple estimate of prime products yield the following corollary for integer-valued matrices.

Multivariate Multipoint Evaluation. Our algorithm is inspired by progress in multivariate multipoint evaluation. Kedlaya and Umans [20] introduced a tabulation-based approach (combined with the Chinese remainder theorem) that later became a key ingredient in fast polynomial composition and factorization. More recently, a sequence of works [4, 2, 3] developed the use of Hasse derivatives and Hermite interpolation to extract more information per evaluation point. We adapt these ideas to our sparse tabulation framework.

Permanents. There exist faster algorithms for computing the permanent in other settings. For sparse matrices, Cygan and Pilipczuk [15] gave a $2^{n-\mathrm{\Omega }(n/d)}$ time algorithm, where $d$ is the average degree of non-zero entries per row. Björklund and Williams [12] gave a $2^{n-\mathrm{\Omega }(n/d^{3/4})}$ time algorithm for $d$-regular bipartite graphs, and a $2^{n-\mathrm{\Omega }(n/r)}$ time algorithm that runs over a finite ring with $r$ elements. Björklund, Husfeldt, and Lyckberg [9] gave a $2^{n-\mathrm{\Omega }(n/(p\log p))}$ time algorithm for computing the permanent modulo a prime power $p^{\lfloor \lambda n/p\rfloor }$, for any constant .

Hamiltonian cycles. There exist faster algorithms for counting Hamiltonian cycles in other settings as well. Björklund, Kaski, and Koutis [10] gave an $O({(2-\delta )}^n)$-time algorithm for counting Hamiltonian cycles modulo moderate prime powers. In the general setting, it is somewhat surprising that our counting algorithm also yields the fastest known algorithm for deciding Hamiltonicity. This differs from the case of permanents: the support of the permanent corresponds to perfect matchings in bipartite graphs, whose existence can be decided in polynomial time. Faster decision algorithms are known in special cases, including Björklund’s $O({1.66}^n)$-time algorithm for undirected graphs [6] and the $O({1.888}^n)$-time algorithm of Cygan, Kratsch, and Nederlof for directed bipartite graphs [14].

For simplicity, we first sketch the case of computing the permanent.

Our improvement comes from combining three ideas:

1. 1.

   Reduce to smaller instances. We reduce the computation on an $n\times n$ matrix to many instances on $k\times k$ matrices. Taking $k$ around $\sqrt{n}$ is what creates room for an improvement in the exponent, provided we can answer the reduced instances fast. This kind of reduction (“self-reduction”) for permanents and Hamiltonian cycles was introduced by Björklund [7]. More concretely, for a parameter $k$, the reduction produces about $2^{n-k}\mathrm{poly}(n)$ instances on $k\times k$ matrices.

2. 2.

   Tabulate only on a sparse set of points. A direct lookup table for all smaller matrices would be far too large. Instead, we employ the fact that permanent is a low-degree polynomial – we tabulate only on a carefully chosen sparse subset of points with the following key property: for every query point, there exists a low-degree univariate curve passing through it whose other points all lie inside the tabulated subset. Then we can recover the value at the query point by interpolating along that curve. (Over finite fields, such subsets are called Kakeya sets. Björklund, Kaski, and Williams [11] were the first to leverage this idea for multivariate polynomial multipoint evaluation in our setting, and we build on their approach.)

3. 3.

   Make each tabulated point more informative. Interpolating from plain point evaluations is limited by how many points we see on a curve. We use the recent idea from Bhargava, Ghosh, Guo, Kumar, and Umans [2, 3] to enrich each tabulated point with additional local information (captured via suitable high-order derivatives), and then use Hermite interpolation to recover higher-degree information along the curve. This is what allows the sparse tabulation approach to only require smaller Kakeya sets, thus working at the larger subproblem sizes we need. (We refer to this task as high-order derivative evaluation.)

Our main technical contribution is a dynamic programming algorithm that makes the required derivative access efficient for the permanent, and we develop analogous data structures for hafnian and Hamiltonian cycles. More concretely, for any constant $ϵ>0$ and a parameter $k$, our data structure tabulates over a Kakeya set of size $2^{O_{ϵ}(k^2)}$¹¹1Here $O_{ϵ}(\cdot )$ means that the hidden constant factor depends on $ϵ$. and takes $O(2^{ϵk})$ time to evaluate one $k\times k$ permanent, so the total time complexity is $O(2^{n-\mathrm{\Omega }(k)}+2^{O(k^2)})$. Balancing the savings from self-reduction with the costs of tabulation yields the final running time $2^{n-\mathrm{\Omega }(\sqrt{n})}$.

The same high-level framework extends to hafnian and Hamiltonian cycles. For hafnian, we derive a suitable self-reduction by modifying components of Björklund’s algorithm [5]. For Hamiltonian cycles, we design an efficient derivative-evaluation data structure based on a determinant characterization [13].

With the tabulation of information on Kakeya sets in $k\times k$ dimensional space, our algorithm essentially computes the hafnian and Hamiltonian cycles in $2^{n-\mathrm{\Omega }(k)}{n}^{O(1)}$ time. It seems that a better construction of a smaller Kakeya set of size $2^{o(k^2)}$ would lead to a faster algorithm. However, the resolution of the finite field Kakeya conjecture [16, 17] rules out such possibilities, showing that the size of a Kakeya set is at least $\mathrm{\Omega }({\delta }^{k^2})$ when the degree of the curve is not greater than $q/\delta$, which corresponds to the regime of our application. Thus, the current construction is essentially optimal for our purposes.

It seems that we have reached a limit with the current approach of self-reduction and Kakeya sets. It remains open whether the techniques on perturbing Ryser’s formula, which work well for sparse permanents [12] and modulo $p^k$ permanents [9], can be adapted to dense permanents to achieve a faster algorithm.

We use $\tilde{O}(f(n))$ to denote $O(f(n)\mathrm{poly}\log (f(n)))$.

Bold symbols like $𝒙$ denote vectors $𝒙=(x_1,\mathrm{\dots },x_n)$.

For any positive integer $n$, we use $[n]$ to denote the set $\{1,2,\mathrm{\dots },n\}$.

We use Iverson’s bracket notation. Let $P$ be a logical proposition, we let $⟦P⟧$ be $1$ if $P$ is true and $0$ otherwise.

With $A\bigsqcup B$, we denote the disjoint union of two sets $A$ and $B$.

For an $n\times m$ matrix $A$, for subsets $S\subseteq [n]$ and $T\subseteq [m]$, we use $A_{S,T}$ to denote the submatrix of $A$ with rows indexed by $S$ and columns indexed by $T$.

Let $\left(\genfrac{}{}{0pt}{}{n}{\downarrow m}\right)$ denote the partial sum of binomial coefficients, i.e.,

We need the estimate of the partial sum of binomials, see [18, Lemma 3.13] for a proof.

We need the following lemma for Hermite interpolation, see [26, Section 5.6] for a proof.

In particular, our algorithm uses the case where those distinct points are the entire finite field ${𝔽}_q$, and $e_i=r$ for all $i$.

Our algorithm uses the Chinese remainder theorem to reduce the problem to small finite fields.

See [26, Section 10.3] for a proof.

We also need an estimate on the product of primes.

See [20, Lemma 2.4] for a proof.

We borrow the self-reduction lemma of Hamiltonian cycles from [7, Lemma 4].

We borrow the definition and construction of Kakeya sets mentioned in [11].

For the convenience of the reader, we provide the construction below, which is originally from [22].

This terminology comes from the intuition in characteristic zero. In that case, computing $P(𝒇(t))mod{t}^r$ is equivalent to computing all the derivatives of $P(𝒇(t))$ up to order $r$.

We rephrase the idea of [2, 3] to reveal information from derivatives.

First, we list the technical ingredients we borrow from Björklund.

We first recap a basic concept introduced in [5, Section 3.1]. For a commutative ring $R$, the set-partition algebra $R[U_m]$ is defined as follows. Intuitively, $U_m$ may be viewed as a partial semigroup whose elements are the subsets of $[m]$, with the operation given by disjoint union. Thus, the product of two subsets $U,V\subseteq [m]$ is defined precisely when $U$ and $V$ are disjoint. The algebra $R[U_m]$ is then the $R$-algebra associated with this partial semigroup, analogous to the group algebra $R[G]$ associated with a group $G$.

More explicitly, every element $r\in R[U_m]$ can be written uniquely as a formal sum

where $r_X\in R$. Addition in $R[U_m]$ is defined componentwise, while multiplication is given by

where the inner sum ranges over all ordered decompositions of $X$ as a disjoint union $Y\bigsqcup Z$.

We identify $R=R[U_0]$ and the inclusion $R[U_m]\subseteq R[U_{m+1}]$ via the inclusion $[m]\subset [m+1]$. Computationally, an element $r\in R[U_m]$ can be stored as the $2^m$ coefficients $r_X\in R$ where $X$ goes through all subsets of $[m]$. The multiplication in $R[U_m]$ is known as the subset convolution, which can be done in $\tilde{O}(2^m)$ $R$-operations [8].

Then, Björklund [5, Section 3.2] introduced a sequence of matrices $B^{(i)}\in R{[U_i]}^{(2n-2i)\times (2n-2i)}$ (for $0\le i\le n$) starting with $B^{(0)}=A$, together with a sequence called squeeze factors ${\beta }^{(i)}\in R[U_i]$ (for $1\le i\le n$). We will not need the precise definition of these matrices and factors, but we invoke the following two facts.

The above lemma has an efficient algorithmic counterpart.²²2Lemma 14 is implicit in Björklund’s original paper. In the original text Björklund merely stated how to sequentially compute $B^{(i)}$ and ${\beta }^{(i)}$ for all $1\le i\le n$ and gave a time complexity analysis, however, in the last paragraph of [5, Section 3.4], the analysis of a single squeeze step is given.

We take several steps to obtain a self-reduction for hafnian. For $1\le \mathrm{\ell }\le n$, we denote $b^{(\mathrm{\ell })}$ to be the prefix product $b^{(\mathrm{\ell })}={\beta }^{(1)}\mathrm{\cdots }{\beta }^{(\mathrm{\ell })}$, by repeatedly applying Lemma 13, we obtain

Then the next step involves extracting the ideas of ranked Möbius transform and inversion from the subset convolution algorithm [8, Section 2].

We consider the ranked Möbius transform, for which for every $X\subseteq [\mathrm{\ell }]$, we introduce the polynomial ${\widehat{b}}_X^{(\mathrm{\ell })}\in R[T]$ defined as

Similarly, we define ${\widehat{B}}_X^{(\mathrm{\ell })}\in R{[T]}^{(2n-2\mathrm{\ell })\times (2n-2\mathrm{\ell })}$ as

For a polynomial $f\in R[T]$, let $[T^d]f$ denote the coefficient of $T^d$ in $f$. We first prove a general statement and then apply it to equation (1).