Recent demonstrations of quantum computational advantage have been driven largely by sampling problems. A prominent model, boson sampling, involves sampling from the output distribution of a linear optical network. However, its classical hardness hinges on two plausible yet less-studied conjectures: the average-case hardness of approximating Gaussian permanents, and the permanent anti-concentration conjecture (PACC). The PACC is a purely mathematical assertion regarding the distributional properties of random Gaussian matrices. While the typical magnitude of the permanent has been established for discrete random matrices, the complex Gaussian case, which governs transition amplitudes in linear optical networks, has remained open. Here, we establish a weak anti-concentration bound by upper-bounding the probability that a random Gaussian permanent is superexponentially smaller than its standard deviation. Tightening this bound to an inverse-polynomial fraction would prove the original PACC. As a corollary, we establish the typical magnitude of Gaussian permanents, on par with Tao and Vu's seminal result for Bernoulli matrices. Combined with the Aaronson-Arkhipov framework, our result implies that classically simulating boson sampling to within a superexponentially small total variation distance would collapse the polynomial hierarchy, assuming the remaining conjectures hold.
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 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, R, Fp and
Q. First, for Haar-distributed unitary matrices~
U, the permanent $\perm(U)$ follows a circularly-symmetric complex Gaussian distribution
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≥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
α≈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.
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×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.
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), 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
per(zJ+W) lie within a disk of radius
O~(n−1/3), which yields an approximation algorithm when the bias of the entries is
Ω~(n−1/3). Previously, there were no efficient algorithms at biases smaller than
1/polylog(n), and it was unknown whether there typically exist zeros
z with
∣z∣≥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.
Exact calculation and even multiplicative error estimation of matrix permanent are challenging for both classical and quantum computers. Regarding the permanents of random Gaussian matrices, the additive error estimation is closely linked to boson sampling and achieving multiplicative error estimation requires exponentially many samplings. Our formula for matrix permanents and its corresponding quantum expression have enabled better estimation of the average additive error for random Gaussian matrices compared to Gurvits' classical sampling algorithm. The well-known Ryser formula has been converted into a quantum permanent estimator. When dealing with real random Gaussian square matrices of size
N, the quantum estimator can approximate the matrix permanent with an additive error smaller than $\ensuremath{\epsilon}{(\sqrt{\mathrm{e}N})}^{N}$, where $\ensuremath{\epsilon}$ is the estimation precision. In contrast, Gurvits' classical sampling algorithm has an estimation error of $\ensuremath{\epsilon}{(2\sqrt{N})}^{N}$, which is exponentially larger (
1.2N) than the quantum method. As expected, the quantum additive error bound fails to reach the multiplicative error bound of ${(2\ensuremath{\pi}N)}^{1/4}\ensuremath{\epsilon}{(\sqrt{N/\mathrm{e}})}^{N}$. Additionally, the quantum permanent estimator can be up to quadratically faster than the classical estimator when using quantum phase estimation-based amplitude estimation.