Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Log inSign up
Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Blog
Job Board
Q3AS

© 2026 Aqora Quantum S.A.S.

TermsPrivacyLegal Notice

arXiv:2610.27872

Moments of Random Permanents and Hafnians: A Replica Approach

Paolo Facchi, Angelo Mariano, Giorgio Parisi, Saverio Pascazio, Antonello Scardicchio

cond-mat.dis-nncond-mat.stat-mechquant-phSubmitted Oct 20, 2026• Announced Oct 21, 2026• 0 citations
arXivHTMLPDF

Abstract

Using contour-integral representations of the permanent and of the Hafnian, we turn the computation of their moments for random matrices into a replica calculation, analogous to that of the moments of the partition function of a mean-field spin glass. For the $2k$-th moment of $N\times N$ matrices with independent entries, the dominant saddle point is replica symmetric for $k>2$, and we obtain in closed form the exponential growth rate, the prefactor and the first $1/N$ correction, in agreement with exact enumerations. The growth rate depends only on whether the entries are real or complex; higher cumulants affect only the prefactor and the $1/N$ term. At $k=2$ the replica-symmetric solution becomes unstable, separating a power-law bulk of the distribution from a large-deviation tail, which we compute. The results apply to boson sampling with single photons and with squeezed light.

Comments
19 pages, 6 figures, 4 tables

Related papers

Permanent of random matrices from representation theory: moments, numerics, concentration, and comments on hardness of boson-sampling

Sepehr Nezami

Computing the distribution of permanents of random matrices has been an outstanding open problem for several decades. In quantum computing, "anti-concentration" of this distribution is an unproven input for the proof of hardness of the task of boson-sampling. We study permanents of random i.i.d. complex Gaussian matrices, and more broadly, submatrices of random unitary matrices. Using a hybrid representation-theoretic and combinatorial approach, we prove strong lower bounds for all moments of the permanent distribution. We provide substantial evidence that our bounds are close to being tight and constitute accurate estimates for the moments. Let $U(d)^{k\times k}$ be the distribution of $k\times k$ submatrices of $d\times d$ random unitary matrices, and $G^{k\times k}$ be the distribution of $k\times k$ complex Gaussian matrices. (1) Using the Schur-Weyl duality (or the Howe duality), we prove an expansion formula for the $2t$-th moment of $|Perm(M)|$ when $M$ is drawn from $U(d)^{k\times k}$ or $G^{k\times k}$. (2) We prove a surprising size-moment duality: the $2t$-th moment of the permanent of random $k\times k$ matrices is equal to the $2k$-th moment of the permanent of $t\times t$ matrices. (3) We design an algorithm to exactly compute high moments of the permanent of small matrices. (4) We prove lower bounds for arbitrary moments of permanents of matrices drawn from $G^{ k\times k}$ or $U(k)$, and conjecture that our lower bounds are close to saturation up to a small multiplicative error. (5) Assuming our conjectures, we use the large deviation theory to compute the tail of the distribution of log-permanent of Gaussian matrices for the first time. (6) We argue that it is unlikely that the permanent distribution can be uniquely determined from the integer moments and one may need to supplement the moment calculations with extra assumptions to prove the anti-concentration conjecture.

arXiv:2104.06423•Apr 13, 2021•0 citations

Permanents of matrix ensembles: computation, distribution, and geometry

Igor Rivin

We report on a computational and experimental study of permanents. On the computational side, we use the GPU to greaatly accelerate the computation of permanents over $\mathbb{C},$ $\mathbb{R},$ $\mathbb{F}_p$ and $\mathbb{Q}.$ First, for Haar-distributed unitary matrices~$U$, the permanent $\perm(U)$ follows a circularly-symmetric complex Gaussian distribution $\mathcal{CN}(0,σ^2)$ -- we confirm this via a number of tests for $n$ up to~23 with $50{,}000$ samples. The DFT matrix permanent is an extreme outlier for every prime $n\ge 7$. In contrast, for Haar-random \emph{orthogonal} matrices~$O$, the permanent $\perm(O)$ is approximately real Gaussian but with positive excess kurtosis that decays as~$O(1/n)$, indicating slower convergence. For matrices with Gaussian entries (GUE, GOE, Ginibre), the permanent follows an $α$-stable distribution with stability index $α\approx 1.0$--$1.4$, well below the Gaussian value $α=2$. We test Aaronson's conjecture that $|\perm(X)|^2$ is asymptotically lognormal for Gaussian~$X$: it is plausible for the complex Ginibre and GOE ensembles, but appears to fail for GUE and real Ginibre, where the $α$-stable tails prevent convergence. Anti-concentration, however, holds for all Gaussian ensembles and is more robust than for Haar unitaries. Secondly, we study the permanent along geodesics on the unitary group. For the geodesic from the identity to the $n$-cycle permutation matrix, we find a universal scaling function $f(t)=\frac{1}{n}\ln|\perm(γ(t))|$ that is independent of~$n$ in the large-$n$ limit, with a midpoint value \[ \perm(γ({\textstyle\frac12})) = (-1)^{(n-1)/2}\cdot 2e^{-n}\bigl(1+\tfrac{1}{3n}+O(n^{-2})\bigr) \] for odd~$n$ and zero for even~$n$. We also study the geodesic forom the identity to the DFT matrix.

quant-pharXiv:2602.10141•Feb 10, 2026

Approximating the Permanent of a Random Matrix with Polynomially Small Mean: Zeros and Universality

Frederic Koehler, Pui Kuen Leung

We study algorithms for approximating the permanent of a random matrix when the entries are slightly biased away from zero. This question is motivated by the goal of understanding the classical complexity of linear optics and \emph{boson sampling} (Aaronson and Arkhipov '11; Eldar and Mehraban '17). Barvinok's interpolation method enables efficient approximation of the permanent, provided one can establish a sufficiently large zero-free region for the polynomial $\mathrm{per}(zJ + W)$, where $J$ is the all-ones matrix and $W$ is a random matrix with independent mean-zero entries. We show that when the entries of $W$ are standard complex Gaussians, all zeros of the random polynomial $\mathrm{per}(zJ + W)$ lie within a disk of radius $\tilde{O}(n^{-1/3})$, which yields an approximation algorithm when the bias of the entries is $\tildeΩ(n^{-1/3})$. Previously, there were no efficient algorithms at biases smaller than $1/\mathrm{polylog}(n)$, and it was unknown whether there typically exist zeros $z$ with $|z| \ge 1$. As a complementary result, we show that the bulk of the zeros, namely $(1 - ε)n$ of them, have magnitude $Θ(n^{-1/2})$. This prevents our interpolation method from contradicting the conjectured average-case hardness of approximating the permanent. We also establish analogous zero-free regions for the hardcore model on general graphs with complex vertex fugacities. In addition, we prove universality results establishing zero-free regions for random matrices $W$ with i.i.d. subexponential entries.

cs.DSCross-listedarXiv:2604.01367•Apr 3, 2026•0 citations

Efficient computation of permanents, with applications to Boson sampling and random matrices

P. H. Lundow, Klas Markström

In order to find the outcome probabilities of quantum mechanical systems like the optical networks underlying Boson sampling, it is necessary to be able to compute the permanents of unitary matrices, a computationally hard task. Here we first discuss how to compute the permanent efficiently on a parallel computer, followed by algorithms which provide an exponential speed-up for sparse matrices and linear run times for matrices of limited bandwidth. The parallel algorithm has been implemented in a freely available software package, also available in an efficient serial version. As part of the timing runs for this package we set a new world record for the matrix order on which a permanent has been computed. Next we perform a simulation study of several conjectures regarding the distribution of the permanent for random matrices. Here we focus on permanent anti-concentration conjecture, which has been used to find the classical computational complexity of Boson sampling. We find a good agreement with the basic versions of these conjectures, and based on our data we propose refined versions of some of them. For small systems we also find noticeable deviations from a proposed strengthening of a bound for the number of photons in a Boson sampling system.

arXiv:1904.06229•Journal of Computational Physics•Jan 19, 2022•21 citations

Photon-number moments and cumulants of Gaussian states

Yanic Cardin, Nicolás Quesada

We develop closed-form expressions for the moments and cumulants of Gaussian states when measured in the photon-number basis. We express the photon-number moments of a Gaussian state in terms of the loop Hafnian, a function that when applied to a ( 0 , 1 ) -matrix representing the adjacency of a graph, counts the number of its perfect matchings. Similarly, we express the photon-number cumulants in terms of the Montrealer, a newly introduced matrix function that when applied to a ( 0 , 1 ) -matrix counts the number of Hamiltonian cycles of that graph. Based on these graph-theoretic connections, we show that the calculation of photon-number moments and cumulants are #P-hard. Moreover, we provide an exponential time algorithm to calculate Montrealers (and thus cumulants), matching well-known results for Hafnians. We then demonstrate that when a uniformly lossy interferometer is fed in every input with identical single-mode Gaussian states with zero displacement, all the odd-order cumulants but the first one are zero. Finally, we employ the expressions we derive to study the distribution of cumulants up to the fourth order for different input states in a Gaussian boson sampling setup where K identical states are fed into an ℓ -mode interferometer. We analyze the dependence of the cumulants as a function of the type of input state, squeezed, lossy squeezed, squashed, or thermal, and as a function of the number of non-vacuum inputs. We find that thermal states perform much worse than other classical states, such as squashed states, at mimicking the photon-number cumulants of lossy or lossless squeezed states.

arXiv:2212.06067•Quantum•Nov 13, 2024•8 citations

Performance of randomized estimators of the Hafnian of a non-negative matrix

Alexey Uvarov, Dmitry Vinichenko

Gaussian boson samplers aim to demonstrate quantum advantage by performing a sampling task believed to be classically hard. The probabilities of individual outcomes in the sampling experiment are determined by the Hafnian of an appropriately constructed symmetric matrix. For non-negative matrices, there is a family of randomized estimators of the Hafnian based on generating a particular random matrix and calculating its determinant. While these estimators are unbiased (the mean of the determinant is equal to the Hafnian of interest), their variance may be so high as to prevent an efficient estimation. Here we investigate the performance of two such estimators, which we call the Barvinok and Godsil-Gutman estimators. We find that, in general, both estimators perform well for adjacency matrices of random graphs, demonstrating a slow growth of variance with the size of the problem. Nonetheless, there are simple examples where both estimators show high variance, requiring an exponential number of samples. In addition, we calculate the asymptotic behavior of the variance for the complete graph. Finally, we simulate the Gaussian boson sampling using the Godsil-Gutman estimator and show that this technique can successfully reproduce low-order correlation functions.

Physical Review A•Apr 17, 2024•0 citations
Using contour-integral representations of the permanent and of the Hafnian, we turn the computation of their moments for random matrices into a replica calculation, analogous to that of the moments of the partition function of a mean-field spin glass. For the 2k2k2k-th moment of N×NN\times NN×N matrices with independent entries, the dominant saddle point is replica symmetric for k>2k>2k>2, and we obtain in closed form the exponential growth rate, the prefactor and the first 1/N1/N1/N correction, in agreement with exact enumerations. The growth rate depends only on whether the entries are real or complex; higher cumulants affect only the prefactor and the 1/N1/N1/N term. At k=2k=2k=2 the replica-symmetric solution becomes unstable, separating a power-law bulk of the distribution from a large-deviation tail, which we compute. The results apply to boson sampling with single photons and with squeezed light.
Computing the distribution of permanents of random matrices has been an outstanding open problem for several decades. In quantum computing, "anti-concentration" of this distribution is an unproven input for the proof of hardness of the task of boson-sampling. We study permanents of random i.i.d. complex Gaussian matrices, and more broadly, submatrices of random unitary matrices. Using a hybrid representation-theoretic and combinatorial approach, we prove strong lower bounds for all moments of the permanent distribution. We provide substantial evidence that our bounds are close to being tight and constitute accurate estimates for the moments. Let U(d)k×kU(d)^{k\times k}U(d)k×k be the distribution of k×kk\times kk×k submatrices of d×dd\times dd×d random unitary matrices, and Gk×kG^{k\times k}Gk×k be the distribution of k×kk\times kk×k complex Gaussian matrices. (1) Using the Schur-Weyl duality (or the Howe duality), we prove an expansion formula for the 2t2t2t-th moment of ∣Perm(M)∣|Perm(M)|∣Perm(M)∣ when MMM is drawn from U(d)k×kU(d)^{k\times k}U(d)k×k or Gk×kG^{k\times k}Gk×k. (2) We prove a surprising size-moment duality: the 2t2t2t-th moment of the permanent of random k×kk\times kk×k matrices is equal to the 2k2k2k-th moment of the permanent of t×tt\times tt×t matrices. (3) We design an algorithm to exactly compute high moments of the permanent of small matrices. (4) We prove lower bounds for arbitrary moments of permanents of matrices drawn from Gk×kG^{ k\times k}Gk×k or U(k)U(k)U(k), and conjecture that our lower bounds are close to saturation up to a small multiplicative error. (5) Assuming our conjectures, we use the large deviation theory to compute the tail of the distribution of log-permanent of Gaussian matrices for the first time. (6) We argue that it is unlikely that the permanent distribution can be uniquely determined from the integer moments and one may need to supplement the moment calculations with extra assumptions to prove the anti-concentration conjecture.
We report on a computational and experimental study of permanents. On the computational side, we use the GPU to greaatly accelerate the computation of permanents over C,\mathbb{C},C, R,\mathbb{R},R, Fp\mathbb{F}_pFp​ and Q.\mathbb{Q}.Q. First, for Haar-distributed unitary matrices~UUU, the permanent $\perm(U)$ follows a circularly-symmetric complex Gaussian distribution CN(0,σ2)\mathcal{CN}(0,σ^2)CN(0,σ2) -- we confirm this via a number of tests for nnn up to~23 with 50,00050{,}00050,000 samples. The DFT matrix permanent is an extreme outlier for every prime n≥7n\ge 7n≥7. In contrast, for Haar-random \emph{orthogonal} matrices~OOO, the permanent $\perm(O)$ is approximately real Gaussian but with positive excess kurtosis that decays as~O(1/n)O(1/n)O(1/n), indicating slower convergence. For matrices with Gaussian entries (GUE, GOE, Ginibre), the permanent follows an ααα-stable distribution with stability index α≈1.0α\approx 1.0α≈1.0--1.41.41.4, well below the Gaussian value α=2α=2α=2. We test Aaronson's conjecture that $|\perm(X)|^2$ is asymptotically lognormal for Gaussian~XXX: it is plausible for the complex Ginibre and GOE ensembles, but appears to fail for GUE and real Ginibre, where the ααα-stable tails prevent convergence. Anti-concentration, however, holds for all Gaussian ensembles and is more robust than for Haar unitaries. Secondly, we study the permanent along geodesics on the unitary group. For the geodesic from the identity to the nnn-cycle permutation matrix, we find a universal scaling function $f(t)=\frac{1}{n}\ln|\perm(γ(t))|$ that is independent of~nnn in the large-nnn limit, with a midpoint value \[ \perm(γ({\textstyle\frac12})) = (-1)^{(n-1)/2}\cdot 2e^{-n}\bigl(1+\tfrac{1}{3n}+O(n^{-2})\bigr) \] for odd~nnn and zero for even~nnn. We also study the geodesic forom the identity to the DFT matrix.
We study algorithms for approximating the permanent of a random matrix when the entries are slightly biased away from zero. This question is motivated by the goal of understanding the classical complexity of linear optics and \emph{boson sampling} (Aaronson and Arkhipov '11; Eldar and Mehraban '17). Barvinok's interpolation method enables efficient approximation of the permanent, provided one can establish a sufficiently large zero-free region for the polynomial per(zJ+W)\mathrm{per}(zJ + W)per(zJ+W), where JJJ is the all-ones matrix and WWW is a random matrix with independent mean-zero entries. We show that when the entries of WWW are standard complex Gaussians, all zeros of the random polynomial per(zJ+W)\mathrm{per}(zJ + W)per(zJ+W) lie within a disk of radius O~(n−1/3)\tilde{O}(n^{-1/3})O~(n−1/3), which yields an approximation algorithm when the bias of the entries is Ω~(n−1/3)\tildeΩ(n^{-1/3})Ω~(n−1/3). Previously, there were no efficient algorithms at biases smaller than 1/polylog(n)1/\mathrm{polylog}(n)1/polylog(n), and it was unknown whether there typically exist zeros zzz with ∣z∣≥1|z| \ge 1∣z∣≥1. As a complementary result, we show that the bulk of the zeros, namely (1−ε)n(1 - ε)n(1−ε)n of them, have magnitude Θ(n−1/2)Θ(n^{-1/2})Θ(n−1/2). This prevents our interpolation method from contradicting the conjectured average-case hardness of approximating the permanent. We also establish analogous zero-free regions for the hardcore model on general graphs with complex vertex fugacities. In addition, we prove universality results establishing zero-free regions for random matrices WWW with i.i.d. subexponential entries.
In order to find the outcome probabilities of quantum mechanical systems like the optical networks underlying Boson sampling, it is necessary to be able to compute the permanents of unitary matrices, a computationally hard task. Here we first discuss how to compute the permanent efficiently on a parallel computer, followed by algorithms which provide an exponential speed-up for sparse matrices and linear run times for matrices of limited bandwidth. The parallel algorithm has been implemented in a freely available software package, also available in an efficient serial version. As part of the timing runs for this package we set a new world record for the matrix order on which a permanent has been computed. Next we perform a simulation study of several conjectures regarding the distribution of the permanent for random matrices. Here we focus on permanent anti-concentration conjecture, which has been used to find the classical computational complexity of Boson sampling. We find a good agreement with the basic versions of these conjectures, and based on our data we propose refined versions of some of them. For small systems we also find noticeable deviations from a proposed strengthening of a bound for the number of photons in a Boson sampling system.
We develop closed-form expressions for the moments and cumulants of Gaussian states when measured in the photon-number basis. We express the photon-number moments of a Gaussian state in terms of the loop Hafnian, a function that when applied to a ( 0 , 1 ) -matrix representing the adjacency of a graph, counts the number of its perfect matchings. Similarly, we express the photon-number cumulants in terms of the Montrealer, a newly introduced matrix function that when applied to a ( 0 , 1 ) -matrix counts the number of Hamiltonian cycles of that graph. Based on these graph-theoretic connections, we show that the calculation of photon-number moments and cumulants are #P-hard. Moreover, we provide an exponential time algorithm to calculate Montrealers (and thus cumulants), matching well-known results for Hafnians. We then demonstrate that when a uniformly lossy interferometer is fed in every input with identical single-mode Gaussian states with zero displacement, all the odd-order cumulants but the first one are zero. Finally, we employ the expressions we derive to study the distribution of cumulants up to the fourth order for different input states in a Gaussian boson sampling setup where K identical states are fed into an ℓ -mode interferometer. We analyze the dependence of the cumulants as a function of the type of input state, squeezed, lossy squeezed, squashed, or thermal, and as a function of the number of non-vacuum inputs. We find that thermal states perform much worse than other classical states, such as squashed states, at mimicking the photon-number cumulants of lossy or lossless squeezed states.
Gaussian boson samplers aim to demonstrate quantum advantage by performing a sampling task believed to be classically hard. The probabilities of individual outcomes in the sampling experiment are determined by the Hafnian of an appropriately constructed symmetric matrix. For non-negative matrices, there is a family of randomized estimators of the Hafnian based on generating a particular random matrix and calculating its determinant. While these estimators are unbiased (the mean of the determinant is equal to the Hafnian of interest), their variance may be so high as to prevent an efficient estimation. Here we investigate the performance of two such estimators, which we call the Barvinok and Godsil-Gutman estimators. We find that, in general, both estimators perform well for adjacency matrices of random graphs, demonstrating a slow growth of variance with the size of the problem. Nonetheless, there are simple examples where both estimators show high variance, requiring an exponential number of samples. In addition, we calculate the asymptotic behavior of the variance for the complete graph. Finally, we simulate the Gaussian boson sampling using the Godsil-Gutman estimator and show that this technique can successfully reproduce low-order correlation functions.