Optimizing and benchmarking the computation of the permanent of general matrices

TL;DR

The paper addresses the challenge of computing the permanent for general (including rectangular) matrices by introducing matrix-permanent, an open-source library that implements the fastest known exact algorithms—combinatoric (naïve), Ryser, and Glynn—and automatically switches among them based on input shape and characteristics. It extends Glynn’s formula to rectangular matrices, documents a robust automatic tuning workflow (including SVM-based boundaries) to dispatch the optimal algorithm, and provides a header-only C++ library plus a Python extension for easy integration. Comprehensive benchmarks on diverse hardware (Apple M1 Pro and Intel Xeon) and matrix classes (random binary and real, identity, ones, and Cauchy) demonstrate that the naïve combinatorial method becomes infeasible quickly, while Ryser excels for rectangular cases and Glynn dominates near-square cases, with precision analyses highlighting stability limits and overflow risks. The work delivers a practical, adaptable tool for researchers across graph theory, quantum many-body physics, and probabilistic point processes, and it is integrated with PyCI for applications in quantum simulations. Overall, matrix-permanent fills a critical gap by providing an efficient, tunable, open-source solution for exact permanent computation across matrix shapes, with broad implications for computational mathematics and physics.

Abstract

Evaluating the permanent of a matrix is a fundamental computation that emerges in many domains, including traditional fields like computational complexity theory, graph theory, many-body quantum theory and emerging disciplines like machine learning and quantum computing. While conceptually simple, evaluating the permanent is extremely challenging: no polynomial-time algorithm is available (unless $\textsc{P} = \textsc{NP}$). To the best of our knowledge there is no publicly available software that automatically uses the most efficient algorithm for computing the permanent. In this work we designed, developed, and investigated the performance of our software package which evaluates the permanent of an arbitrary rectangular matrix, supporting three algorithms generally regarded as the fastest while giving the exact solution (the straightforward combinatoric algorithm, the Ryser algorithm, and the Glynn algorithm) and, optionally, automatically switching to the optimal algorithm based on the type and dimensionality of the input matrix. To do this, we developed an extension of the Glynn algorithm to rectangular matrices. Our free and open-source software package is distributed via Github, at https://github.com/theochem/matrix-permanent.

Figures (6)

• Figure 1: The fastest algorithm by matrix order, $n$, and the degree of rectangularity, $\tfrac{m}{n}$. The values reported indicate the performance (as a factor of execution time) that would be lost were the second-fastest algorithm used instead of the fastest algorithm.
• Figure 2: The space of fastest algorithms for computing the permanent of a matrix. The separating hyperplanes between the algorithms are depicted by the dotted lines. The hyperplanes intersect at matrix order $n = 13$ and rectangularity $m/n = 0.29$
• Figure 3: A comparison of the algorithm's execution time for evaluating the permanent of real-valued matrices with elements randomly selected from the interval [-1, 1]. The varying matrix aspect ratios are shown by the coloured lines for each algorithm. For the naiv̈e (combinatoric) algorithm, matrices with more than 14 columns were not considered because the time ($> 15$ minutes) exceeded that for all other algorithms ($< 1$ second) by three orders of magnitude.
• Figure 4: Accuracy of the Glynn algorithm for $n$-by-$n$ identity matrix with double (a) and integer (b) type. The error is assessed using Eq. (\ref{['eq:digits_lost']}). No other algorithms are displayed as the accuracy remained stable within the testing period.
• Figure 5: Accuracy of the Ryser and Glynn algorithm for $n$-by-$n$ ones matrix (every element is one) with double (a) and integer (b) type. The error is assessed using Eq. (\ref{['eq:digits_lost']}). Overflow occurs at the same moment for all matrix aspect ratios for the Glynn algorithm.