On the classical complexity of quantum interference of indistinguishable bosons
The classical complexity of sampling from the probability distribution of quantum interference of indistinguishable single bosons on unitary network with input and output ports is studied with the focus on how boson density for affects the number of computations required to produce a single sample. Glynn’s formula is modified for computation of probabilities of output configurations with output ports occupied by bosons, requiring only computations. It is found that in a unitary network chosen according to the Haar probability measure the tails of the distribution of the total number of output ports occupied by bosons are bounded by those of a binomial distribution. This fact allows to prove that for any with probability the number of computations in the classical sampling algorithm of P. Clifford and R. Clifford scales as at least and at most , where and .
Boson Sampling idea of Aaronson & Arkhipov  links sampling from the output distribution of the many-body quantum interference of indistinguishable single bosons on a unitary linear -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 , when an output port receives at most a singe boson , 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 and is impossible . The amplitudes of -boson interferences are given as matrix permanents of -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 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 .
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 , 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 .
Boson Sampling can be easier than the universal quantum computation  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 , superconducting qubits , neutral atoms in optical lattices  and dynamic Casimir effect . 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 .
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., . A significant reduction in the number of classical computations was demonstrated by using a Markov Chain Monte Carlo classical simulation algorithm  and subsequent analytical estimate  with the threshold for the quantum advantage elevated from the original  to bosons. This result is essentially the threshold of classical computability of a single matrix permanent by Glynn’s method , which allows simultaneous computation of matrix permanents appearing in the chain rule of conditional probabilities (lemma 2 of Ref. ) 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.  is general, the estimate on the classical complexity can be further reduced beyond the no-collision regime.
A finite boson density regime can be experimentally advantageous, since using a much smaller network, with only ports instead of at least , 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 .
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  is proposed, featuring the same speed-up for the matrix permanents of matrices with repeated columns/rows as found previously in Ryser’s method . The classical hardness of a probability at the output of Boson Sampling operating at a non-vanishing boson density regime 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 in the Haar measure (theorem 1 below). Moreover, since a modified Glynn method is employed for computation of matrix permanents, the algorithm of Ref.  (applicable uniformly over all regimes of boson density ) can be used for estimating the number of computations for classical sampling from the output distribution of Boson Sampling for any boson density . Thus, our lower and upper bounds also apply to the classical complexity of sampling from the output distribution of Boson Sampling for arbitrary .
2 Modified Glynn’s formula for matrix permanent of a matrix with repeated columns or rows
We consider quantum interference of perfectly indistinguishable single bosons on a unitary network with input and output ports, below fixing the input ports to be . We use the notations: for the density of bosons,
For a general multi-set with coinciding elements, the corresponding submatrix is rank-deficient, which reduces the computational complexity of its matrix permanent . 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  by analysing Ryser’s algorithm . Below, however, we will need the number of computations according to Glynn’s algorithm , since this algorithm is used for analysing the sampling complexity in Ref. . The algorithm itself has to be modified to obtain the same speed up as in Ref. .
Let us assume that only output ports are occupied by bosons and set (without loosing the generality) . Introducing for each an auxiliary complex variable taking values
we can rewrite the matrix permanent in the output probability Eq. (1) as follows
Indeed, the r.h.s. of Eq. (2) is a sum over all permutations in the product of elements of with the multi-set of columns corresponding to a configuration , where for a whole number , since there is at least one factor for all on the r.h.s of Eq. (2) and each auxiliary variable satisfies
where is the remainder of division . On the other hand, since the inner sum on the r.h.s. of Eq. (2) contributes factors consisting of products of exactly elements, i.e., , to the outer sum, we must have
resulting in all , i.e., .
Similar as in Glynn’s formula , a reduction of the number of summations is still possible in Eq. (2). Assuming that is the minimum of the non-zero (i.e., for ) one can set and omit the summation over in Eq. (2). Indeed, let us show that
satisfied for only if () for all .
Now let us estimate the number of computations in Eq. (4). There are multiplications of sums of matrix elements in the rows . The number of additions in the inner sum (over ) can be reduced, similarly as in Refs. [8, 40], by ordering the factors in the outer sum in such a way that neighbour factors have just one element different. Then, for each such factor, the computation of the inner sum requires only one addition and one multiplication (to change one factor ) for each term in the outer sum in Eq. (4). Therefore, the total number of computations 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 ( is the estimate found in Ref. ), is significant only for large lower bound of the non-zero occupations of output ports. For instance, Eq. (5) correctly estimates the number of computations in the case of maximally bunched output . However, on average, in a Haar-random network , the reduction in the number of computations according to Eq. (4) as compared to Eq. (2) is at most by . To show that, consider the probability of maximal bunching of bosons at the output of a Haar-random network, given as 
Setting in Eq. (6) the probability to be we obtain
One can conclude that, with probability arbitrarily close to , on average, the denominator in Eq. (5) is at most .
3 Estimate on classical computations for producing samples from the output distribution of many-boson interference
To quantify the number of classical computations for general , when an arbitrary output configuration 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 . 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 probability.
Given the number of output ports occupied by bosons, the number of classical computations necessary to compute an output probability , 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., and for , whereas the upper bound by uniformly distributing bosons over the occupied output ports, i.e., for .
Consider the probability distribution of the number of output ports occupied by bosons in a Haar-random network. Thanks to the uniform average probability the probability of output ports occupied by bosons is obviously
Indeed, we choose uniformly randomly out of output ports and distribute in them bosons with the uniform probability over the output configurations. The average number of the output ports occupied by bosons reads
As (for a fixed ) the tails of the distribution (8) (precise definition of the tails is given below) are below by the tails of the following binomial distribution
For the two distributions of Eqs. (8) and (10) have bell-shaped form. Indeed, whereas the binomial is evidently bell-shaped about , 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 is attained when , i.e., at the average value of Eq. (9). Moreover, by using Stirling’s approximation for the factorials , one can easily check that for
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 from the left and from the right of the maximum the binomial distribution dominates that of the total number of occupied ports.
Let us formally define the lower and upper tails by the points of equal probability . Now, let us show that , where there are such when (or no upper tail of Eq. (8) exists if ). To find the points consider the equation
Applying the standard approximation to factorials  in the binomial , after simple algebra we obtain the following asymptotic equation for
where , which is asymptotically equivalent (for and ) to
Since the r.h.s of Eq. (11) varies between and , i.e., the minimum and maximum values of for , there is always a solution for for any (no upper tail of Eq. (8) exists for ). Therefore, we have shown that there such , with some , that for and for the binomial dominates the distribution of Eq. (8).
Now we can use the Hoeffding-Chernoff bound , which states that the tails of the binomial distribution of Eq. (10), i.e., the sum of probabilities for and , with , are bounded by . We get for as above
For any with the probability at least in the Haar measure over unitary networks the number of classical computations required to compute (by a modified Glynn algorithm of section 2) a probability of the output distribution of quantum interference of indistinguishable bosons on a randomly chosen unitary network with input and output ports satisfies as
In theorem 1 we have taken into account also the vanishing density case, as : if for some the upper bound is (the distribution in Eq. (8) does not possess the upper tail).
In Ref.  the estimate (the leading order) on the number of classical computations required to produce a single sample from the output distribution of quantum interference of single bosons on a unitary -port was given. Note that the crucial fact is that the algorithm of Ref.  applies uniformly over all possible output configurations, since multi-sets of output ports are obtained by consecutive sampling from a conditional probability chain rule used there, each time sampling for one output port. The estimates of Ref.  corresponds to the classical hardness to compute a single permanent by Glynn’s formula , 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. . Namely, lemma 2 of Ref.  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.  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 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.
-  S. Aaronson and A. Arkhipov, Theory of Computing 9, 143 (2013).
-  A. Arkhipov and G. Kuperberg, Geometry & Topology Monographs 18, 1 (2012).
-  E. R. Caianiello, Nuovo Cimento, 10, 1634 (1953); Combinatorics and Renormalization in Quantum Field Theory, Frontiers in Physics, Lecture Note Series (W. A. Benjamin, Reading, MA, 1973).
-  S. Scheel, arXiv:quant-ph/0406127 (2004).
-  L. G. Valiant, Theoretical Coput. Sci., 8, 189 (1979).
-  S. Aaronson, Proc. Roy. Soc. London A, 467, 3393 (2011).
-  H. Ryser, Combinatorial Mathematics, Carus Mathematical Monograph No. 14. (Wiley, 1963).
-  D. Glynn, Eur. J. of Combinatorics, 31, 1887 (2010).
-  M. Jerrum, A. Sinclair, and E. Vigoda, Journal of the ACM 51, 671 (2004).
-  M. J. Bremner, A. Montanaro, and D. J. Shepherd, Quantum 1, 8 (2017).
-  J. Bermejo-Vega, D. Hangleiter, M. Schwarz, R. Raussendorf, and J. Eisert, Phys. Rev. X 8, 021010 (2018).
-  S. Boixo et al, Nature Physics, 14, 595 (2018).
-  X. Gao, S.-T. Wang, and L.-M. Duan, Phys. Rev. Lett. 118, 040502 (2017).
-  J. Preskill, arXiv:1801.00862 (2018).
-  R. Feynman, Int. J. Theoret. Phys. 21, 467 (1982).
-  P. W. Shor, Proceedings of the 35th Annual Symposium Foundations of Computer Science (IEEE, New York, 1994), p. 124.
-  G. Kalai, Notices of the AMS, 63, 508 (2016).
-  E. Knill, R. Laflamme, and G. J. Millburn, Nature 409, 46 (2001).
-  M. A. Broome, A. Fedrizzi, S. Rahimi-Keshari, J. Dove, S. Aaronson, T. C. Ralph, and A. G. White, Science 339, 794 (2013);
-  J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Kolthammer, X.-M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K. Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R. Smith, and I. A. Walmsley, Science, 339, 798 (2013).
-  M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit, and P. Walther, Nature Photonics, 7, 540 (2013).
-  A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F. Galvão, N. Spagnolo, C. Vitelli, E. Maiorino, P. Mataloni, and F. Sciarrino, Nature Photonics, 7, 545 (2013).
-  A. P. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L. O’Brien, and T. C. Ralph, Phys. Rev. Lett. 113, 100502 (2014).
-  M. Bentivegna et al, Sci. Adv. 1, e1400255 (2015).
-  K. R. Motes, A. Gilchrist, J. P. Dowling, and P. P. Rohde, Phys. Rev. Lett. 113, 120501 (2014).
-  Y. He et al, Phys. Rev. Lett. 118, 190501 (2017).
-  J. C. Loredo et al, Phys. Rev. Lett. 118, 130503 (2017).
-  H. Wang et al, Phys. Rev. Lett. 120, 230502 (2018).
-  H.-S. Zhong et al, Phys. Rev. Lett. 121, 250505 (2018).
-  C. Shen, Z. Zhang, and L.-M. Duan, Phys. Rev. Lett. 112, 050504 (2014).
-  B. Peropadre, G. G. Guerreschi, J. Huh, and A. Aspuru-Guzik Phys. Rev. Lett. 117, 140505 (2016).
-  A. Deshpande, B. Fefferman, M. C. Tran, M. Foss-Feig and A. V. Gorshkov, Phys. Rev. Lett. 121, 030501 (2018).
-  B. Peropadre, J. Huk and C. Sabín, Sci. Reps 8, 3751 (2018).
-  S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Phys. Rev. X 6, 021039 (2016).
-  J. J. Renema, A. Menssen, W. R. Clements, G. Triginer, W. S. Kolthammer, and I. A. Walmsley, Phys. Rev. Lett. 120, 220502 (2018).
-  M. Oszmaniec and D. J Brod, New J. Phys. 20, 092002 (2018).
-  R. García-Patrón, J. J. Renema, and V. S. Shchesnovich, arXiv:1712.10037 [quant-ph].
-  J. J. Renema, V. S. Shchesnovich, and R. García-Patrón, arXiv:1809.01953 [quant-ph].
-  A. Neville, C. Sparrow, R. Clifford, E. Johnston, P. M. Birchall, A. Montanaro, A. Laing, Nature Physics 13, 1153 (2017).
-  P. Clifford, and R. Clifford, arXiv:1706.01260 (2017).
-  M. C. Tichy, Entanglement and Interference of Identical Particles. PhD thesis, Freiburg University (2011); see Appendix B.
-  V. S. Shchesnovich, Int. J. Quantum Inform. 11, 1350045 (2013); see appendix D.
-  A. I. Barvinok, Math. Oper. Res., 21 65 (1996); see theorem (3.3).
-  H. Minc, Permanents, Encyclopedia of Mathematics and Its Applications, Vol. 6 (Addison-Wesley Publ. Co., Reading, Mass., 1978).
-  C. Mortici, J. Math. Ineqs. 5, 611 (2011).
-  W. Hoeffding, J. of the Amer. Stat. Assoc. 58, 13 (1963).