On the classical complexity of quantum interference of indistinguishable bosons

Boson Sampling idea of Aaronson & Arkhipov [1] links sampling from the output distribution of the many-body quantum interference of $N$ indistinguishable single bosons on a unitary linear $M$-port chosen at random from the Haar measure and a mathematical problem of estimating matrix permanents of matrices with elements being independent identically distributed complex Gaussian random variables. The relation is due the fact that in the so-called no-collision regime, with at least $M\gg N^2$, when an output port receives at most a singe boson [2], matrix elements of a unitary network are approximated by independent identically distributed complex Gaussian random variables. Under plausible conjectures a classical simulation of the above sampling task to a given error $ϵ$ with the computations polynomial in $N$ and $1/ϵ$ is impossible [1]. The amplitudes of $N$-boson interferences are given as matrix permanents of $N$-dimensional submatrices of a unitary network matrix [3, 4], whose computation is believed to be classically hard [5, 6]. The fastest known algorithms [7, 8] compute a matrix permanent of an arbitrary matrix in $O\left(N{2}^N\right)$ computations. In contrast, distinguishable bosons (classical particles) result in an output probability distribution expressed as matrix permanents of positive matrices, which can be estimated polynomially in $N,1/ϵ$ [9].

Boson Sampling, is among several proposals [10, 11, 12, 13] considered for the quantum supremacy demonstration, i.e., demonstration of a provable computational advantage of a quantum system over digital computers [14], the first step on the way of using the computational advantage promised by quantum mechanics [15, 16]. The importance of such a demonstration cannot be underestimated in view of the opposite hypothesis [17].

Boson Sampling can be easier than the universal quantum computation [18] implemented with linear optics, since it requires neither interaction between photons nor error correction schemes, the proof of principle experiments [19, 20, 21, 22] followed the initial idea, moreover, improvements and advances are constantly reported [23, 24, 25, 26, 27, 28, 29]. Moreover, alternative platforms include ion traps [30], superconducting qubits [31], neutral atoms in optical lattices [32] and dynamic Casimir effect [33]. On the other hand, inevitable experimental imperfections can allow for efficient classical simulations algorithms [34, 35, 36, 37, 38], with current optical platform offering only a finite-size window for the quantum supremacy demonstration with Boson Sampling due to scaling up of boson losses with the size of a planar optical network [37].

Before any demonstration of quantum supremacy is attempted, it is important to know how exactly the classical computations required to sample from the output distribution of Boson Sampling scale up with the size of the system, i.e., $N,M$. A significant reduction in the number of classical computations was demonstrated by using a Markov Chain Monte Carlo classical simulation algorithm [39] and subsequent analytical estimate [40] with the threshold $N$ for the quantum advantage elevated from the original $N\approx 30$ [1] to $N\approx 50$ bosons. This result is essentially the threshold of classical computability of a single matrix permanent by Glynn’s method [8], which allows simultaneous computation of matrix permanents appearing in the chain rule of conditional probabilities (lemma 2 of Ref. [40]) with the largest size of a matrix equal to the total number of bosons. However Glynn’s formula ignores reduction of the complexity when boson bunching at output has non-vanishing probability, thus though the approach of Ref. [40] is general, the estimate on the classical complexity can be further reduced beyond the no-collision regime.

A finite boson density regime $\rho =\frac{N}{M}=O\left(1\right)$ can be experimentally advantageous, since using a much smaller network, with only $M=O\left(N\right)$ ports instead of at least $M\gg N^2$ [1], reduces the effect of boson losses on the computational hardness [36, 37, 38]. On the other hand, bunching of bosons at output of a unitary network, increasing with the boson density, affects computational hardness of the corresponding output probabilities [41, 42] due to repeated columns (corresponding to output ports of a network) in the matrices giving quantum interference amplitudes and the fact that the computational hardness of a general matrix permanent depends on the matrix rank [43].

To estimate how exactly boson density affects the computational hardness of many-body bosonic interference on a linear unitary multiport is the goal of the present work. First, a modification of Glynn’s method [8] is proposed, featuring the same speed-up for the matrix permanents of matrices with repeated columns/rows as found previously in Ryser’s method [42]. The classical hardness of a probability at the output of Boson Sampling operating at a non-vanishing boson density regime $\rho =O\left(1\right)$ depends on the distribution of bosons over the output ports (the output configuration). Therefore, to quantify the classical complexity of the output probability distribution, the lower and upper bounds on the number of computations in the modified Glynn’s method are used. The bounds depend on the total number of output ports occupied by bosons, the crucial fact is that in a Haar-random network the distribution of the total number of ports occupied by bosons has a bell-shaped form, with the tails bounded by those of a binomial distribution. This fact allows to state lower and upper bounds with probability arbitrarily close to $1$ in the Haar measure (theorem 1 below). Moreover, since a modified Glynn method is employed for computation of matrix permanents, the algorithm of Ref. [40] (applicable uniformly over all regimes of boson density $0\le \rho \le 1$) can be used for estimating the number of computations for classical sampling from the output distribution of Boson Sampling for any boson density $0\le \rho \le 1$. Thus, our lower and upper bounds also apply to the classical complexity of sampling from the output distribution of Boson Sampling for arbitrary $M\ge N$.

We consider quantum interference of $N$ perfectly indistinguishable single bosons on a unitary network $U$ with $M$ input and output ports, below fixing the input ports to be $k=1,\dots ,N$. We use the notations: $\rho =N/M$ for the density of bosons,

for the probability [1] of detecting bosons in a multi-set of output ports $l_1\le \dots \le l_N$ corresponding to output configuration $m=(m_1,\dots ,m_M)$, $m_1+\dots +m_M=N$, $U[1,\dots ,N|l_1,\dots ,l_N]$ for the submatrix of $U$ on the rows $1,\dots ,N$ and columns $l_1,\dots ,l_N$, $per(\dots )$ for the matrix permanent [44].

For a general multi-set $l_1,\dots ,l_N$ with coinciding elements, the corresponding submatrix $U[1,\dots ,N|l_1,\dots ,l_N]$ is rank-deficient, which reduces the computational complexity of its matrix permanent [43]. An estimate of the number of classical computations for evaluation of such a matrix permanent with account of the speed up due to reduced matrix rank was found before [42] by analysing Ryser’s algorithm [7]. Below, however, we will need the number of computations according to Glynn’s algorithm [8], since this algorithm is used for analysing the sampling complexity in Ref. [40]. The algorithm itself has to be modified to obtain the same speed up as in Ref. [42].

Let us assume that only $n$ output ports are occupied by bosons and set (without loosing the generality) $m_{n+1}=\dots m_M=0$. Introducing for each $l=1,\dots ,n$ an auxiliary complex variable $x_l$ taking $m_l+1$ values

we can rewrite the matrix permanent in the output probability $p\left(m\right)$ Eq. (1) as follows

Indeed, the r.h.s. of Eq. (2) is a sum over all permutations in the product of $N$ elements of $U[1,\dots N|l_1^{\prime },\dots ,l_N^{\prime }]$ with the multi-set of columns $l_1^{\prime },\dots ,l_N^{\prime }$ corresponding to a configuration $m^{\prime }=(m_1^{\prime },\dots ,m_n^{\prime },0,\dots ,0)$, where $m_l^{\prime }=m_l+{\Delta }_l(m_l+1)$ for a whole number ${\Delta }_l\ge 0$, since there is at least one factor $x_l$ for all $l=1,\dots ,n$ on the r.h.s of Eq. (2) and each auxiliary variable $x_l$ satisfies

where $rem(s,q)$ is the remainder of division $s/q$. On the other hand, since the inner sum on the r.h.s. of Eq. (2) contributes factors consisting of products of exactly $N$ elements, i.e., $x_{l_1^{\prime }}\dots x_{l_N^{\prime }}{U}_{1l_1^{\prime }}\dots U_{N{l}_N^{\prime }}$, to the outer sum, we must have

resulting in all ${\Delta }_l=0$, i.e., $m_l^{\prime }=m_l$.

Similar as in Glynn’s formula [8], a reduction of the number of summations is still possible in Eq. (2). Assuming that $m_n$ is the minimum of the non-zero $m_l$ (i.e., $m_n\le m_l$ for $l=1,\dots ,n$) one can set $x_n=1$ and omit the summation over $x_n$ in Eq. (2). Indeed, let us show that

In this case, due to an arbitrary number of factors $U_{kn}$ with different $k$ in the product on the r.h.s. of Eq. (4), the equality condition in Eq. (3) for $m_1^{\prime },\dots ,m_n^{\prime }$ becomes an inequality for $m_1^{\prime },\dots ,m_{n-1}^{\prime }$

satisfied for $m_l\ge m_n$ only if ${\Delta }_l=0$ ($m_l^{\prime }=m_l$) for all $l=1,\dots ,n-1$.

Now let us estimate the number of computations in Eq. (4). There are $N$ multiplications of sums of matrix elements in the rows $1,\dots ,N$. The number of additions in the inner sum (over $l=1,\dots ,n$) can be reduced, similarly as in Refs. [8, 40], by ordering the factors $x_{l_1}\dots x_{l_{n-1}}$ in the outer sum in such a way that neighbour factors have just one element $x_{l_j}$ different. Then, for each such factor, the computation of the inner sum requires only one addition and one multiplication (to change one factor $x_{l_j}$) for each term in the outer sum in Eq. (4). Therefore, the total number of computations $C_m$ in Eq. (4) is defined solely by the outer sum and the product and reads

The reduction of the number of computations in Eq. (4) as compared to Eq. (2), where for the latter by the same counting we get $C_m^{\prime }\equiv O(N{\prod }_{l=1}^n(m_j+1))$ ($C_m^{\prime }$ is the estimate found in Ref. [42]), is significant only for large lower bound $min\left(m_l\right)$ of the non-zero occupations $m_l$ of output ports. For instance, Eq. (5) correctly estimates the number of computations $C_m=O\left(N\right)$ in the case of maximally bunched output $m=(N,0,\dots ,0)$. However, on average, in a Haar-random network $U$, the reduction in the number of computations according to Eq. (4) as compared to Eq. (2) is at most by $O\left(\ln N\right)$. To show that, consider the probability of maximal bunching of bosons at the output of a Haar-random network, given as [2]

Setting in Eq. (6) the probability to be $1-ϵ$ we obtain

One can conclude that, with probability arbitrarily close to $1$, on average, the denominator in Eq. (5) is at most $O\left(\ln N\right)$.

To quantify the number of classical computations for general $N\ge N$, when an arbitrary output configuration $m=(m_1,\dots ,m_M)$ occurs with a non-vanishing probability (with the uniform average probability in a Haar-random network), let us study how the total number of output ports occupied by bosons is distributed. Though the exact distribution is hard to figure out in an arbitrary given network, in a Haar-random network it has a narrow bell-shaped form, with the tails bounded by a binomial distribution for $N\gg 1$. This allows one to get an almost sure lower and upper bounds on the computational complexity for any network, i.e., the bounds apply to a randomly chosen network according to the Haar measure with arbitrarily close to $1$ probability.

Given the number $1\le n\le N$ of output ports occupied by bosons, the number of classical computations necessary to compute an output probability $p\left(m\right)$, according to the algorithm in Eq. (4), satisfies

where the lower bound is obtained by assuming the maximal bunching in a single output port, e.g., $m_1=N-n+1$ and $m_l=1$ for $l=2,\dots ,n-1$, whereas the upper bound by uniformly distributing bosons over the occupied output ports, i.e., $m_l=N/n$ for $l=1,\dots ,n$.

For $N\gg 1$ the two distributions of Eqs. (8) and (10) have bell-shaped form. Indeed, whereas the binomial is evidently bell-shaped about $n=⟨n⟩=Mx=N/(1+\rho )$, the distribution of the number of output ports occupied by bosons satisfies

where the two factors in each of the two relations come from the respective binomial coefficients. Therefore, the only extremum (maximum) probability for $N\gg 1$ is attained when $n^2\approx (M-n)(N-n)$, i.e., at the average value of Eq. (9). Moreover, by using Stirling’s approximation for the factorials [45], one can easily check that for $n=N/(1+\rho )$

i.e., the maximum of the binomial is below that of the distribution of the number of output ports occupied by bosons, which fact suggests that at some points $n_{\mp }$ from the left and from the right of the maximum $n=⟨n⟩$ the binomial distribution dominates that of the total number of occupied ports.

Let us formally define the lower $1\le n\le n_{-}$ and upper $n_{+}\le n\le N$ tails by the points $n_{\mp }$ of equal probability $P\left(n_{\pm }\right)=B_{n_{\pm }}\left(x\right)$. Now, let us show that $n_{\pm }=\frac{1\pm {\delta }_{\pm }}{1+\rho }N$, where there are such ${\delta }_{\pm }>0$ when $N\gg 1$ (or no upper tail of Eq. (8) exists if ${\delta }_{+}>\rho$). To find the points $n_{\pm }$ consider the equation

Applying the standard approximation to factorials [45] in the binomial $(\frac{N-1}{n-1})$, after simple algebra we obtain the following asymptotic equation for ${\delta }_{\pm }$

where $H\left(z\right)=-z\ln z-(1-z)\ln (1-z)$, which is asymptotically equivalent (for ${\delta }_{-}<1$ and ${\delta }_{+}<\rho$) to

Since the r.h.s of Eq. (11) varies between $0$ and $\ln 2$, i.e., the minimum and maximum values of $H\left(z\right)$ for $0\le z\le 1$, there is always a solution for ${\delta }_{\pm }$ for any $0\le \rho \le 1$ (no upper tail of Eq. (8) exists for ${\delta }_{+}>\rho$). Therefore, we have shown that there such $n_{\pm }=\frac{1\pm \delta }{1+\rho }N$, with some $\delta =max\left({\delta }_{\pm }\right)\in (0,\rho )$, that for $1\le n\le n_{-}$ and for $n_{+}\le n\le N$ the binomial dominates the distribution of Eq. (8).

Now we can use the Hoeffding-Chernoff bound [46], which states that the tails of the binomial distribution of Eq. (10), i.e., the sum of probabilities for $0\le n\le n_{-}$ and $n_{+}\le n\le M$, with $n_{\pm }=(1\pm \delta )Mx=\frac{1\pm \delta }{1+\rho }N$, are bounded by $e^{-{\delta }^{2}Mx/4}$. We get for $\delta$ as above

Now, setting $\exp \left(-\frac{{\delta }^{2}N}{4(1+\rho )}\right)=ϵ/2$ we get $\delta =O\left(1/\sqrt{N}\right)$ thus for sufficiently large $N\gg 1$ such $\delta$ will satisfy Eq. (11). Therefore, using the lower and upper bounds $n=n_{\pm }=\frac{1\pm \delta }{1+\rho }N$ in Eq. (7) we have the following theorem.

In theorem 1 we have taken into account also the vanishing density case, $\rho \to 0$ as $N\to \infty$: if $\delta >\rho$ for some $N\gg 1$ the upper bound is $O\left(N{2}^N\right)$ (the distribution in Eq. (8) does not possess the upper tail).

In Ref. [40] the estimate $O\left(N{2}^N\right)$ (the leading order) on the number of classical computations required to produce a single sample from the output distribution of quantum interference of $N$ single bosons on a unitary $M$-port was given. Note that the crucial fact is that the algorithm of Ref. [40] applies uniformly over all possible output configurations, since multi-sets of output ports $l_1\le \dots \le l_N$ are obtained by consecutive sampling from a conditional probability chain rule used there, each time sampling for one output port. The estimates of Ref. [40] corresponds to the classical hardness to compute a single permanent by Glynn’s formula [8], which disregards the speed up achieved by the modification of section 2, leading to the lower and upper bounds given by theorem 1. The modified Glynn method and theorem 1 are sufficient prerequisites for applying the algorithm of Ref. [40]. Namely, lemma 2 of Ref. [40] is extendable to the modified Glynn’s method, whereas a uniform estimate on the computational hardness of computing an output probability allows one to estimate the computations required for sampling the output ports using the conditional probability chain rule. Since the estimate of Ref. [40] in the leading order is essentially the classical hardness of computing a single matrix permanent by Glynn’s method, the number of computations (in the leading order) to produce a single sample of Boson Sampling in any density regime $0\le \rho \le 1$ satisfies the bounds in Eq. (13).

This work was supported by the National Council for Scientific and Technological Development (CNPq) of Brazil, grant 304129/2015-1 and by the São Paulo Research Foundation (FAPESP), grant 2018/24664-9.