Randomized Classical Hafnian Estimators

• Randomized classical hafnian estimators are polynomial-time stochastic algorithms that compute the hafnian of even-dimensional symmetric nonnegative matrices, effectively counting perfect matchings.
• They construct a random skew-symmetric matrix from the input and use its determinant as an unbiased estimator, with error controlled by graph expansion properties.
• These estimators find key applications in combinatorial enumeration and quantum linear optics, with empirical results showing manageable variance for typical graph ensembles.

Randomized classical hafnian estimators are polynomial-time stochastic algorithms for approximating the hafnian of an even-dimensional symmetric nonnegative matrix. The hafnian $\haf(A)$, for $A\in\mathbb R^{n\times n}$ with $n=2m$ and $A_{ii}=0$, enumerates the perfect matchings in the graph with adjacency matrix $A$, and plays a critical role in combinatorial enumeration and quantum linear optics. Central among known classical estimators are the Barvinok and Godsil–Gutman estimators, both of which construct a random skew-symmetric matrix from $A$ and take its determinant as an unbiased estimator of $\haf(A)$ (Rudelson et al., 2014, Uvarov et al., 2023).

1. Mathematical Definition of the Hafnian and Estimators

Given $A\in\mathbb R^{2m\times 2m}$ with $A_{ij}=A_{ji}\ge 0$ and $A_{ii}=0$, the hafnian is defined as

$A\in\mathbb R^{n\times n}$0

or, equivalently, as the sum over all perfect matchings $A\in\mathbb R^{n\times n}$1 of $A\in\mathbb R^{n\times n}$2: $A\in\mathbb R^{n\times n}$3 Both estimators construct a random matrix $A\in\mathbb R^{n\times n}$4 with $A\in\mathbb R^{n\times n}$5 for $A\in\mathbb R^{n\times n}$6 and $A\in\mathbb R^{n\times n}$7. The determinant $A\in\mathbb R^{n\times n}$8 is then an unbiased estimator: $A\in\mathbb R^{n\times n}$9

• Barvinok estimator: $n=2m$0, i.i.d.
• Godsil–Gutman estimator: $n=2m$1 i.i.d. uniform on $n=2m$2

Unbiasedness stems from the independence and zero-mean nature of the $n=2m$3, ensuring that only perfect-matching monomials contribute in expectation (Rudelson et al., 2014, Uvarov et al., 2023).

2. Algorithmic Workflow and Computational Complexity

The standard algorithm for these estimators comprises the following steps:

1. Generate independent random variables $n=2m$4 for $n=2m$5 (Gaussian for Barvinok, sign for Godsil–Gutman).
2. Form the matrix $n=2m$6 ($n=2m$7, $n=2m$8).
3. Compute $n=2m$9 (runtime $A_{ii}=0$0).

A single sample $A_{ii}=0$1 delivers a random variable concentrated near $A_{ii}=0$2 in favorable regimes. Repeating the process and averaging yields success probability increasing exponentially in the number of samples. The overall complexity for any desired polynomially small error remains polynomial in $A_{ii}=0$3 when the variance is only polynomially large (Rudelson et al., 2014, Uvarov et al., 2023).

3. Variance, Error Bounds, and Expansion Criteria

Although $A_{ii}=0$4 always holds, the variance may be large: $A_{ii}=0$5 The exact combinatorial form for $A_{ii}=0$6 is

$A_{ii}=0$7

where $A_{ii}=0$8 runs over perfect 2-matchings, $A_{ii}=0$9 (Barvinok), $A$0 (Godsil–Gutman), and the cycles are all even (Uvarov et al., 2023).

Variance is manageable under combinatorial "strong expansion" of the large-variance graph $A$1: for a suitable threshold $A$2, $A$3 is the graph with edges where $A$4. $A$5 is strongly expanding up to level $A$6 with parameter $A$7 if for each $A$8,

$A$9

where $A$0 is the vertex boundary and $A$1 is the number of connected components. Under these conditions and with sufficient spectral gap, the error is subexponential in $A$2: $A$3 with high probability, improving to $A$4 error if a spectral gap exists (Rudelson et al., 2014).

4. Asymptotic Analysis and Sample Complexity

The efficiency of the estimators hinges on the growth of the relative standard deviation $A$5. For the complete graph $A$6, asymptotics yield

$A$7

so $A$8 for Barvinok ($A$9) and $\haf(A)$0 for Godsil–Gutman ($\haf(A)$1) (Uvarov et al., 2023). For random graphs $\haf(A)$2 and typical inputs, $\haf(A)$3 grows sub-$\haf(A)$4.

The number of samples $\haf(A)$5 to achieve relative error $\haf(A)$6 with constant confidence is

$\haf(A)$7

If $\haf(A)$8, then $\haf(A)$9 and a fully polynomial randomized approximation scheme (FPRAS) is achieved. However, if $A\in\mathbb R^{2m\times 2m}$0 is exponentially large (as in worst-case graph constructions), exponential sampling is required and classical estimation is infeasible (Uvarov et al., 2023).

5. Empirical Performance and Numerical Experiments

Empirical tests confirm the theoretical picture. For “typical” random graphs (Erdős–Rényi, moderate to large $A\in\mathbb R^{2m\times 2m}$1), both estimators yield polynomially growing variance and sample complexity. For the complete graph, observed $A\in\mathbb R^{2m\times 2m}$2 matches $A\in\mathbb R^{2m\times 2m}$3 scaling. Numerical experiments with Godsil–Gutman’s estimator for Gaussian Boson Sampling reveal that increasing the number of samples $A\in\mathbb R^{2m\times 2m}$4 rapidly improves accuracy, consistent with theoretical $A\in\mathbb R^{2m\times 2m}$5 error decay. Low-order correlation functions in simulated Gaussian Boson Sampling can be estimated efficiently with polynomial resources for adjacency matrices with nonnegative entries (Uvarov et al., 2023).

Special graph structures, such as unions of short cycles or bridged components, demonstrate exponentially growing variance for both estimators, as predicted by the combinatorial expansion of the second moment. These cases limit the practical applicability of the estimators to certain graph classes unless preprocessing, variance reduction, or graph decomposition is employed (Uvarov et al., 2023).

6. Extensions, Limitations, and Open Directions

Both estimators extend naturally to weighted nonnegative matrices, provided suitable scaling (e.g., doubly-stochastic normalization) is applied to maintain the expansion property and avoid large outliers. The analysis holds for weighted matching counts under analogous expansion and no-large-entry hypotheses (Rudelson et al., 2014). In the absence of strong expansion, such as in sparse graphs of small minimum degree, estimators can exhibit exponential fluctuations, rendering them ineffective.

Unlike MCMC-based permanent estimators, the randomized hafnian estimators require no sophisticated Markov-chain mixing: each sample is sharply concentrated around the mean in favorable regimes (Rudelson et al., 2014). Identifying specific “worst-case” graph structures or introducing variance-reducing randomizations—such as alternative ensembles for $A\in\mathbb R^{2m\times 2m}$6—remains an open research avenue (Uvarov et al., 2023). Further refinement may involve graph decomposition, spectral preprocessing, or biased sampling to minimize higher moments of $A\in\mathbb R^{2m\times 2m}$7.

7. Applications in Combinatorics and Quantum Optics

Randomized classical hafnian estimators provide efficient, polynomial-time algorithms for the approximate enumeration of perfect matchings in graphs, up to a subexponential factor $A\in\mathbb R^{2m\times 2m}$8, or even $A\in\mathbb R^{2m\times 2m}$9 in the presence of a spectral gap (Rudelson et al., 2014). For nonnegative kernels in Gaussian Boson Sampling, they enable accurate simulation of correlation functions for moderate system sizes, suggesting no exponential quantum-classical separation for physically relevant observables in this regime (Uvarov et al., 2023). The connection between combinatorial expansion, variance control, and quantum sampling complexity positions these estimators as a central tool in contemporary research at the interface of combinatorics, probability, and quantum computation.