Matrix Hafnians: Theory, Complexity, and Applications

• Matrix hafnians are combinatorial functions defined for even-dimensional symmetric matrices that sum over all perfect matchings using the matrix entries as weights.
• They are applied in fields like quantum optics and combinatorics, serving to compute probabilities in Gaussian Boson Sampling and enumerate graph matchings.
• Computing matrix hafnians is #P-hard, and advanced algorithms use techniques such as inclusion–exclusion and generating functions to address their computational complexity.

A matrix hafnian is a combinatorial matrix function that, for an even-dimensional symmetric matrix, encodes the sum over all perfect matchings of the complete graph with edge weights given by the matrix entries. Hafnians generalize the permanent (with applications in bipartite graphs) to nonbipartite, undirected settings and are fundamental in a variety of fields, including combinatorics, quantum optics (notably Gaussian Boson Sampling), and statistical physics. Efficient computation and sharp bounds for hafnians remain central subjects in computational complexity and algorithm design, with connections to #P-completeness and quantum advantage proposals.

1. Formal Definition and Key Properties

Given a symmetric matrix $A\in\mathbb{C}^{2n\times 2n}$, the hafnian is defined by

$$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$

where the sum runs over the $(2n)!/(2^n n!)$ unordered partitions of $\{1,\dots,2n\}$ into $n$ disjoint unordered pairs. The function is zero unless $A$ is of even order, and is fully symmetric in the matrix entries.

The hafnian generalizes the notion of the permanent, with the precise algebraic relationship: $$\mathrm{per}(B) = \mathrm{haf}\left(\begin{array}{cc} 0 & B \ B^T & 0 \end{array}\right)$$ for any $n\times n$ matrix $B$ (Barvinok, 2014, Barvinok, 2016).

For general $N\times N$ (not necessarily even) matrices, the loop hafnian or $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$0 is defined as the sum over all "single-pair matchings" (partitions into pairs and singletons), allowing for loops (diagonal terms) (Bulmer et al., 2021, Tarasov, 21 Jul 2025). For zero-diagonal $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$1, $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$2 reduces to $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$3.

2. Connections to Graph Theory and Combinatorics

For an unweighted, undirected graph $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$4 on $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$5 vertices with adjacency matrix $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$6, the hafnian $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$7 counts the number of perfect matchings in $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$8, i.e., sets of $$\mathrm{haf}(A) = \sum_{\text{pair partitions } \{i_1,j_1\},\dots,\{i_n,j_n\}=\{1,\dots,2n\}} \prod_{k=1}^n a_{i_kj_k}$$9 edges covering each vertex exactly once and with no shared vertices (Rudelson et al., 2014, Efimov, 2021). For graphs with edge weights, the hafnian becomes the sum over all perfect matchings of the product of edge weights in the matching.

In combinatorics, hafnians naturally enumerate certain chord diagrams, dimer coverings, and appear in the enumeration of matchings in regular and random graphs. Exact formulas and combinatorial interpretations exist for special classes of matrices, such as certain Toeplitz types and two-parameter models, yielding closed expressions in terms of matchings of subtemplate graphs (Efimov, 2021, Efimov, 2019).

3. Computational Complexity and Algorithms

Computing the hafnian of a general symmetric matrix is #P-hard, mirroring the complexity of the permanent (Barvinok, 2016, Cardin et al., 2022). For an $(2n)!/(2^n n!)$0 matrix, classical algorithms require exponential time. The fastest known exact algorithms, e.g., the Björklund–Gupt–Quesada method, achieve $(2n)!/(2^n n!)$1 time by exploiting inclusion–exclusion and power trace techniques, and are strongly parallelizable (Björklund et al., 2018, Bulmer et al., 2021). The algorithm's cubic prefactor per $(2n)!/(2^n n!)$2 subset is asymptotically optimal, unless one assumes major breakthroughs in permanent computation.

For quasi-polynomial-time (in $(2n)!/(2^n n!)$3) approximate computation, Barvinok established that if $(2n)!/(2^n n!)$4 meets entrywise analytic conditions (e.g., $(2n)!/(2^n n!)$5 for all $(2n)!/(2^n n!)$6), then $(2n)!/(2^n n!)$7 can be deterministically approximated within any prescribed relative error $(2n)!/(2^n n!)$8 in $(2n)!/(2^n n!)$9 time by Taylor expansion of a suitable analytic generating function and complex-zeros analysis (Barvinok, 2014, Barvinok, 2016). This framework can also be extended to the computation of the logarithm of the hafnian for matrices with entries bounded in $\{1,\dots,2n\}$0 for arbitrary fixed $\{1,\dots,2n\}$1 (Barvinok, 2016).

For random nonnegative matrices with strong expansion (highly connected underlying graphs), a randomized estimator via determinants of random skew-symmetric Gaussian matrices provides a polynomial-time scheme with subexponential error guarantees, under expansion assumptions analogous to those required for probabilistic permanent estimators (Rudelson et al., 2014).

Table: Complexity of Key Hafnian Algorithms

Algorithm/Estimator	Time Complexity	Error/Regime
Björklund–Gupt–Quesada exact method	$\{1,\dots,2n\}$2	Exact, general matrices
Barvinok quasi-poly approx (analytic)	$\{1,\dots,2n\}$3	Relative error $\{1,\dots,2n\}$4; entrywise analytic conditions
Randomized (Barvinok estimator)	Poly($\{1,\dots,2n\}$5)	Subexp error, strong expansion

For loop hafnians, similarly structured algorithms exist, and additional optimized schemes exploit collision structure and repeated row/column patterns to reduce prefactors in Gaussian Boson Sampling simulations (Bulmer et al., 2021, Quesada et al., 2019).

4. Analytic, Combinatorial, and Generating Function Techniques

Hafnians admit expansion and reduction identities, including a row-expansion (analogous to Laplace expansion for determinants) and a sum-of-matrices identity: $\{1,\dots,2n\}$6 where $\{1,\dots,2n\}$7 is the principal submatrix indexed by $\{1,\dots,2n\}$8, and $\{1,\dots,2n\}$9 is its complement (Efimov, 2021, Efimov, 2019).

Generating function methods, notably via Gaussian integration, link (loop) hafnians directly to exponential moments of quadratic forms (Tarasov, 21 Jul 2025). For a symmetric matrix $n$0, the generating function

$n$1

generates all loop hafnians of block expansions of $n$2; the coefficient of $n$3 yields $n$4 for repeated row/column blocks as appearing in Gaussian state output probability amplitudes (Tarasov, 21 Jul 2025).

For certain structured families, such as special Toeplitz matrices, closed-form analytic expressions for hafnians are available, with complexity reduced to $n$5 or $n$6 for matrix dimension $n$7 (Efimov, 2019, Efimov, 2021).

5. Inequalities, Extremal Problems, and Asymptotics

Matrix hafnians satisfy nontrivial upper and lower bounds, with significant interest in analogues of van der Waerden’s and Tverberg’s theorems for permanents.

• Upper bounds: A generalized Laplace expansion combined with subset-convolution yields multiplicative upper bounds for hafnians, including, e.g.,

$n$8

which can outperform Fischer’s permanent-based bound, especially for matrices with sign cancellations or complex structure (Roos, 2019).

• Lower bounds: Extension of hyperbolic polynomial techniques enables explicit lower bounds on the minimal hafnian over convex hulls of symmetric permutation matrices with zero diagonal. For the uniform matrix $n$9,

$A$0

which achieves asymptotic value $A$1 as $A$2 (Friedland, 2011).

While local minimality for this bound is established for all $A$3-matchings, global minimality fails for large $A$4 due to explicit expander graph constructions, making the precise identification of global extremizers open for many regimes.

6. Applications in Quantum Optics and Gaussian Boson Sampling

Matrix hafnians and, crucially, loop hafnians are central in the calculation of output amplitudes and probabilities for quantum-optical experiments based on Gaussian states and photon-number measurements. In Gaussian Boson Sampling (GBS), the probability amplitude for a given detection pattern is proportional to the loop hafnian of a matrix derived from the state’s covariance matrix and displacement (for general detection patterns with collisions or loss) (Bulmer et al., 2021, Cardin et al., 2022, Quesada et al., 2019). For the ideal zero-displacement case and no collisions, the hafnian suffices.

Recent developments include:

• Explicit demonstration of nonnegativity and monotonicity properties for hafnians of physical GBS-derivable matrices with block structure $A$5, where $A$6 is Hermitian and positive semidefinite (Bradler et al., 2018).
• Development of fast classical algorithms and eigenvalue-sieve acceleration schemes enabling GBS simulation with much larger photon/mode numbers than previously possible, closing the quantum-classical computational divide for several experimental regimes (Bulmer et al., 2021).
• Tailored loop hafnian algorithms, exploiting repeated structure, enable orders-of-magnitude acceleration for heralded state preparation, underpinning simulation studies of non-Gaussian state generation under loss and practical device imperfections (Quesada et al., 2019).

7. Extensions and Open Problems

Significant extensions include multidimensional hafnians (hyperhafnians), direct connections to moment-generating functions in Gaussian analysis, and structural links to other combinatorial matrix functions such as the Montrealer (encoding cumulants and Hamiltonian cycles) (Cardin et al., 2022).

Open problems are abundant: optimal global lower bounds for hafnians on convex polytopes of fixed degree, further reduction of exponential prefactors in practical algorithms, and closed combinatorial characterizations for new classes of structured matrices. The boundary between classical and quantum computation for hafnian-based problems remains a dynamic research front (Bulmer et al., 2021).

---

In summary, matrix hafnians provide the fundamental, computationally hard, and algorithmically rich bridge between graph matchings and modern quantum-combinatorial phenomena. Their study integrates analytic, algebraic, computational, and physical methods and continues to inform the frontiers of both mathematical theory and quantum technology.