On the classical complexity of sampling from quantum interference of indistinguishable bosons

TL;DR

This paper analyzes the classical computational complexity of simulating boson sampling experiments, providing bounds that relate the difficulty of sampling with smaller interferometers to larger ones, crucial for assessing quantum advantage.

Contribution

It offers a new complexity estimate for boson sampling on smaller interferometers, extending previous results to near-collision regimes with high probability.

Findings

Classical sampling complexity for N bosons on M ports is at least as hard as sampling with N/(1+N/M) bosons in the no-collision regime.

Sampling complexity scales with the ratio N/(1+N/M), linking small and large interferometer regimes.

Results are applicable to current experimental setups with fewer photons and ports, informing the pursuit of quantum advantage.

Abstract

Experimental demonstration of the quantum advantage over classical simulations with Boson Sampling is currently under intensive investigation. There seems to be a scalability issue to the necessary number of bosons on the linear optical platforms and the experiments, such as the recent Boson Sampling with $20$ photons on $60$-port interferometer by H.~Wang~\textit{et al}, \textit{Phys. Rev. Lett.} \textbf{123,} 250503 (2019), are usually carried out on a small interferometer, much smaller than the size necessary for the no-collision regime. Before demonstration of quantum advantage, it is urgent to estimate exactly how the classical computations necessary for sampling from the output distribution of Boson Sampling are reduced when a smaller-size interferometer is used. The present work supplies such a result, valid with arbitrarily close to $1$ probability, which reduces in the no-collision regime to the previous estimate by P.~Clifford and R.~Clifford. One of the results with immediate application to current experiments with Boson Sampling is that classically sampling from the interference of $N$ single bosons on an $M$-port interferometer is at least as hard as that with $\mathcal{N}= \frac{N}{1+N/M}$ single bosons in the no-collision regime, i.e., on a much larger interferometer with at least $\mathcal{M}\gg N^2$ ports.

Equations

$$\hat{p}(\mathbf{m})=\frac{|\mathrm{per}(U[1,\dots ,N|l_1,\dots ,l_N]){|}^2}{m_1!\dots m_l!}$$

$${\mathcal{C}}_{\mathbf{m}}=O\left(N^2\prod _{l=1}^M(m_l+1)\right).$$

$$x_l\in \{1,e^{\frac{2i\pi }{m_l+1}},\dots ,e^{\frac{2i\pi m_l}{m_l+1}}\},$$

$$\mathrm{per}(U[1,\dots N|l_1,\dots ,l_N])=\frac{{\prod }_{l=1}^n{m}_l!}{{\prod }_{l=1}^n(m_l+1)}\sum _{x_1}\dots \sum _{x_n}{x}_1\dots x_n\prod _{k=1}^N\left(\sum _{l=1}^n{x}_l{U}_{kl}\right).$$

$$\sum_{x_{l}}x^{m^{\prime}+1}_{l}=\sum_{j=0}^{m_{l}}e^{\frac{2i\pi(m^{\prime}+1)}{m_{l}+1}j}=\left\{\begin{array}[]{cc}0,&\mathrm{rem}(m^{\prime}+1,m_{l}+1)\neq 0\\ m_{l}+1,&\mathrm{rem}(m^{\prime}+1,m+1)=0\end{array}\right.,$$

$$N=\sum _{l=1}^n{m}_l^{\prime }=\sum _{l=1}^n{\Delta }_l(m_l+1)+N,$$

$$\mathrm{per}(U[1,\dots N|l_1,\dots ,l_N])=\frac{{\prod }_{l=1}^n{m}_l!}{{\prod }_{l=1}^{n-1}(m_l+1)}\sum _{x_1}\dots \sum _{x_{n-1}}{x}_1\dots x_{n-1}\prod _{k=1}^N\left(\sum _{l=1}^{n-1}{x}_l{U}_{kl}+U_{kn}\right).$$

$$N\ge \sum _{l=1}^{n-1}{m}_l^{\prime }=\sum _{l=1}^{n-1}{\Delta }_l(m_l+1)+N-m_n,$$

$${\mathcal{C}}_{\mathbf{m}}\equiv O\left(N\frac{{\prod }_{l=1}^n(m_l+1)}{\min (m_l+1)}\right).$$

$$\mathrm{Prob}(\mathrm{max}(m_{l})\leq m)\approx\left[1-\Bigl{(}\frac{\rho}{1+\rho}\Bigr{)}^{m+1}\right]^{M}.$$

$$m\lesssim \frac{\ln \left(\frac{N}{\rho ϵ}\right)}{\ln \left(\frac{1+\rho }{\rho }\right)}.$$

$$\Omega \left(N{2}^n\right)\le {\mathcal{C}}_{\mathbf{m}}\le O\left(N{\left(1+N/n\right)}^n\right),$$

$$P(n)\equiv \mathrm{Prob}\left(\sum _{l=1}^M\theta (m_l)=n\right)=\frac{\left(\genfrac{}{}{0ex}{}{M}{n}\right)\left(\genfrac{}{}{0ex}{}{N-1}{n-1}\right)}{\left(\genfrac{}{}{0ex}{}{M+N-1}{M-1}\right)},1\le n\le N.$$

$$⟨n⟩=\frac{MN}{M+N-1}=\frac{N}{1+\rho }\left[1+O\left(\frac{1}{M+N}\right)\right].$$

$$B_n(x)=\left(\genfrac{}{}{0ex}{}{M}{n}\right)x^n(1-x{)}^{M-n},x\equiv \frac{\rho }{1+\rho }.$$

$$P(n+1)=\frac{M-n}{n+1}\frac{N-n}{n}P(n),P(n-1)=\frac{n}{M-n+1}\frac{n-1}{N-n+1}P(n),$$

$$\frac{\left(\genfrac{}{}{0ex}{}{N-1}{n-1}\right)}{\left(\genfrac{}{}{0ex}{}{M+N-1}{M-1}\right)}=\sqrt{\frac{1+\rho }{\rho }}\frac{{\rho }^n}{(1+\rho {)}^M}\left[1+O\left(\frac{1}{N}\right)\right],x^n(1-x{)}^{M-n}=\frac{{\rho }^n}{(1+\rho {)}^M},$$

$${\left(\genfrac{}{}{0ex}{}{M}{n}\right)}^{-1}P(n)=x^n(1-x{)}^{M-n}.$$

$$\sqrt{\frac{(1\pm {\delta }_{\pm })(1+\rho )}{\rho \mp {\delta }_{\pm }}}\exp \left\{NH\left(\frac{1\pm {\delta }_{\pm }}{1+\rho }\right)\right\}=(1+\rho {)}^N\left[1+O\left(\frac{1}{N}\right)\right],$$

$$H\left(\frac{1\pm {\delta }_{\pm }}{1+\rho }\right)=\ln (1+\rho ).$$

$$\sum _{n=n_{-}}^{n_{+}}P(n)\ge \sum _{n=n_{-}}^{n_{+}}{B}_n(x)>1-2\exp \left(-\frac{{\delta }^{2}N}{4(1+\rho )}\right).$$

$$\Omega \left(N{2}^{\frac{1-\delta }{1+\rho }N}\right)\le {\mathcal{C}}_{\mathbf{m}}\le O\left(N{\left(1+r\right)}^{\frac{N}{r}}\right),\delta =2\sqrt{\frac{1+\rho }{N}\ln \left(\frac{2}{ϵ}\right)},r=\max \left(1,\frac{1+\rho }{1+\delta }\right).$$

$$\sum _{\mathbf{m}}f(l_1,\dots ,l_N)=\sum _{l_1=1}^M\dots \sum _{l_N=1}^M\frac{m_1!\dots m_M!}{N!}f(l_1,\dots ,l_N),$$

$$p(l_1,\dots ,l_N)=\frac{m_1!\dots m_M!}{N!}\hat{p}(\mathbf{m})=\frac{1}{N!}{\left|\mathrm{per}(U[1,\dots ,N|l_1,\dots ,l_N])\right|}^2.$$

$$p(l_1,\dots ,l_K)=\sum _{l_{K+1}=1}^M\dots \sum _{l_N=1}^{M}p(l_1,\dots ,l_N)=\frac{1}{N!}\sum _{\sigma ,\tau }\left[\prod _{k=1}^K{U}_{\sigma (k),l_k}{U}_{\tau (k),l_k}^{\ast }\right]\prod _{k=K+1}^N{\delta }_{\sigma (k),\tau (k)}.$$

$$p(l_1,\dots ,l_K)=\frac{1}{N!}\sum _{\pi }\frac{{\left|\mathrm{per}(U[\pi (1),\dots ,\pi (K)|l_1,\dots ,l_K])\right|}^2}{K!}.$$

$$p(l_1,\dots ,l_K)=\sum _{\pi }p(\pi )p(l_1,\dots ,l_K|\pi ),p(l_1,\dots ,l_K|\pi )\equiv \frac{{\left|\mathrm{per}(U[\pi (1),\dots ,\pi (K)|l_1,\dots ,l_K])\right|}^2}{K!}.$$

$$p(l_1,\dots ,l_N)=\sum _{\pi }p(\pi )p(l_1,\dots ,l_N|\pi ),$$

$$p(l_1,\dots ,l_K|\pi )=p(l_1|\pi )p(l_2|l_1;\pi )\dots p(l_K|l_1,\dots ,l_{K-1};\pi ),p(l_r|l_1,\dots ,l_{r-1};\pi )\equiv \frac{p(l_1,\dots ,l_r|\pi )}{p(l_1,\dots ,l_{r-1}|\pi )}.$$

$$p(l_1,\dots ,l_K|\pi )=\frac{1}{K!}{\left|\sum _{\alpha =1}^K{U}_{\pi (\alpha ),l_K}\mathrm{per}(U[\pi (1),\dots ,\pi (\alpha -1),\pi (\alpha +1),\dots ,\pi (K)|l_1,\dots ,l_{K-1}])\right|}^2$$

Full text

1 Introduction

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 each output port receives at most a singe boson [2], matrix elements of the unitary matrix describing a multiport 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 $\epsilon$ with the computations polynomial in $N$ and $1/\epsilon$ is impossible [1]. The classical hardness originates from the fact, that the amplitudes of $N$-boson interferences are given as matrix permanents of $N$-dimensional submatrices of a unitary matrix [3, 4], whose computation is believed to be classically hard [5, 6] (see also a review [7]). The fastest known algorithms [8, 9] compute a matrix permanent of an arbitrary matrix in $O(N2^{N})$ computations. Even a relative error approximation to the absolute value of a matrix permanent is classically hard [1] (given two conjectures are true; see also Ref. [10]). Approximation algorithms to an additive error are known: a probabilistic approximation to matrix permanent of an arbitrary matrix [11], generalised for matrices with repeated rows or columns [12] and with repeated rows and columns [13]. On the other hand, deterministic approximation algorithm to a relative error was found only for diagonally dominant matrices [14]. In contrast, distinguishable bosons (classical particles) result in an output probability distribution expressed as the matrix permanents of positive matrices, which can be estimated polynomially in $N,1/\epsilon$ [15] and approximated by the deterministic algorithm of Ref. [12].

Boson Sampling, is among several proposals [16, 17, 18, 19] considered for the quantum supremacy demonstration, i.e., demonstration of a provable computational advantage of a quantum system over digital computers [20], the first step on the way of using the computational advantage promised by quantum mechanics [21, 22]. The importance of such a demonstration cannot be underestimated in view of the opposite hypothesis [23].

Boson Sampling can be easier than the universal quantum computation [24] implemented with linear optics, since it requires neither interaction between photons nor error correction schemes. The proof of principle experiments [25, 26, 27, 28] have followed the initial idea, moreover, improvements and advances are constantly reported [29, 30, 31, 32, 33, 34, 35]. Recently, an experimental Boson Sampling with $20$ photons on $60$-port interferometer have been demonstrated [36]. Moreover, alternative platforms include ion traps [37], superconducting qubits [38], neutral atoms in optical lattices [39] and dynamic Casimir effect [40].

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 [41] and subsequent analytical estimate [42] 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 [9], which allows simultaneous computation of matrix permanents appearing in the chain rule of conditional probabilities with the largest size of a matrix equal to the total number of bosons.

On the other hand, inevitable experimental imperfections can allow for further speed up, moreover, even efficient classical simulations algorithms are possible [43, 44, 45, 46, 47]. Currently used planar optical platforms have an exponential scaling of boson losses with the optical depth of a multiport (i.e., the total number of layers of beamsplitters and phase shifters), leaving only a finite-size window in optical depth for the quantum supremacy demonstration with Boson Sampling [46] and making a large planar optical multiport $M\gg N^{2}$ for $N\sim 50$ [41, 42] unsuitable for the quantum advantage demonstration due to strong losses.

One is therefore forced to either search for other platforms for experimental Boson Sampling, which do not have exponentially scaling losses of bosons, or consider sending more photons to an interferometer up to linear scaling of the interferometer size $M\sim N$, not discarded as a possibility [1] for the quantum advantage demonstration. Current experiments use smaller interferometers than necessary for the no-collision regime (due to losses or other reasons) [34, 35, 36]. In such regimes, however, bunching of bosons at output of a unitary multiport reduces the computational hardness of the corresponding output probabilities [48, 49, 50], due to boson bunching at the output ports and the fact that the computational hardness of a general matrix permanent depends on the matrix rank [51]. Thus though the previous algorithm [41, 42] can be applicable to such a regime, the number of computations can be further reduced due to boson bunching and an appropriate estimates of this effect is in order.

To derive an estimate on the computational hardness valid for arbitrary $M\geq N$, or arbitrary density of bosons $\rho={N}/{M}$, is the goal of the present work. In contrast to the no-collision regime, where the number of computations for an output probability does not depend on the output distribution of bosons, thanks to vanishing probability of boson bunching at the output ports, for a finite-density regime $\rho=\Omega(1)$ it depends strongly on the occupations of the output ports [48, 49, 50]. Averaging of the number of computations over the $\binom{M+N-1}{N}$ output distributions of bosons seems to be an unfeasible task, since the output probabilities are hard to estimate. In an experiment a fixed multiport must be chosen uniformly randomly over the unitary group (i.e., according to the Haar measure) [1], but averaging over the Haar-random unitary is not equivalent to averaging over the output distributions of bosons for a fixed multiport. Thus, estimating the number of computations for a fixed multiport for general $M\geq N$ is a difficult problem. This problem is solved below by observing that the distribution of the total number of output ports occupied by bosons, averaged over the Haar-random unitary, is a narrow one of the bell-shaped form. Since the tails of such a distribution are exponentially small, we cut off two small regions corresponding to unitary matrices on the low and high ends of the computational complexity, and estimate the number of computations from above and below for the rest of multiports. In this way, the resulting estimates are valid for a fixed multiport with a small probability of failure. This approach allows us to obtain an experimentally relevant result that applies to a fixed multiport (except a small fraction of multiports at the low end of the computational complexity, such as the multiports described by diagonally dominant unitary matrices [14]): simulating Boson Sampling with $N$ bosons on $M$-port interferometer is at least as hard as with $\mathcal{N}=N/(1+N/M)$ bosons in the no-collision regime, i.e., on a $\mathcal{M}$-port with at least $\mathcal{M}\gg N^{2}$. Though such an estimate on the number of computations does not substitute an actual complexity-theoretical proof of hardness of Boson Sampling in a finite-density regime, it is an important result for the experimental efforts to scale up Boson Sampling, while simultaneously combat losses of bosons.

The rest of the text is organized as follows. In section 2 a modification of Glynn’s method [9] is proposed, with the same speed-up for the matrix permanents of matrices with repeated columns/rows as found previously in Ryser’s method [49]. The classical hardness of a probability at the output of Boson Sampling operating in a finite-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 multiport 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 of section 3. Moreover, since a variant of Glynn’s method is employed for computation of matrix permanents, the algorithm of Ref. [42] (applicable uniformly over all regimes of density of bosons) is used to give a reduced estimate on the number of computations required to produce a single sample from the output distribution of Boson Sampling in any density regime $0<\rho\leq 1$, section 4, where theorem 2 gives the main result. In the last section 5 open problems related to the presented results are discussed.

2 Modified Glynn’s formula for matrix permanent of a matrix with repeated columns or rows

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

[TABLE]

for the probability [1] of detecting bosons in a multi-set of output ports $l_{1}\leq\ldots\leq l_{N}$ corresponding to output configuration $\mathbf{m}=(m_{1},\ldots,m_{M})$, $m_{1}+\ldots+m_{M}=N$, $U[1,\ldots,N|l_{1},\ldots,l_{N}]$ for the submatrix of $U$ on the rows $1,\ldots,N$ and columns $l_{1},\ldots,l_{N}$, $\mathrm{per}(\ldots)$ for the matrix permanent [52].

For a general multi-set $l_{1},\ldots,l_{N}$ with coinciding elements, the corresponding submatrix $U[1,\ldots,N|l_{1},\ldots,l_{N}]$ is rank-deficient, which reduces the computational complexity of its matrix permanent [51]. An estimate of the number of classical computations $\mathcal{C}_{\mathbf{m}}$ for evaluation of such a matrix permanent with account of the speed up due to reduced matrix rank was found before [49] by analysing Ryser’s algorithm [8]. For $N\leq M$ we have

[TABLE]

The essential result of Eq. (2) was reproduced also in Ref. [50]. Below, however, we will need the number of computations according to Glynn’s algorithm [9], since this algorithm is used for analysing the sampling complexity in Ref. [42]. The algorithm itself is modified below, in a similar way as in Ref. [12], to reduce the number of computations to the bare necessary.

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

[TABLE]

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

[TABLE]

Indeed, the r.h.s. of Eq. (3) is a sum over all permutations in the product of $N$ elements of $U[1,\ldots N|l^{\prime}_{1},\ldots,l^{\prime}_{N}]$ with the multi-set of columns $l^{\prime}_{1},\ldots,l^{\prime}_{N}$ corresponding to a configuration $\mathbf{m}^{\prime}=(m^{\prime}_{1},\ldots,m^{\prime}_{n},0,\ldots,0)$, divided by $m^{\prime}_{1}!\ldots m^{\prime}_{n}!$, where $m^{\prime}_{l}=m_{l}+\Delta_{l}(m_{l}+1)$ for a whole number $\Delta_{l}\geq 0$, since there is at least one factor $x_{l}$ for all $l=1,\ldots,n$ on the r.h.s of Eq. (3) and each auxiliary variable $x_{l}$ satisfies

[TABLE]

where $\mathrm{rem}(s,q)$ is the remainder of division $s/q$. Since the product of summations over columns of $U$ on the r.h.s. of Eq. (3) contributes factors consisting of products of exactly $N$ elements, i.e., $x_{l^{\prime}_{1}}\ldots x_{l^{\prime}_{N}}U_{1l^{\prime}_{1}}\ldots U_{Nl^{\prime}_{N}}$, to the summations over $x_{1},\ldots,x_{n}$, we must have

[TABLE]

resulting in all $\Delta_{l}=0$, i.e., $m^{\prime}_{l}=m_{l}$. Eq. (3) is an alternative generalised Glynn estimator to that found in Ref. [12].

Similar as in Glynn’s formula [9], a reduction of the number of summations is still possible in Eq. (3). Assuming that $m_{n}$ is the minimum of the non-zero $m_{l}$ (i.e., $m_{n}\leq m_{l}$ for $l=1,\ldots,n$) one can set $x_{n}=1$ and omit the summation over $x_{n}$ in Eq. (3). Indeed, let us show that

[TABLE]

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. (5), the equality condition in Eq. (4) for $m^{\prime}_{1},\ldots,m^{\prime}_{n}$ becomes an inequality for $m^{\prime}_{1},\ldots,m^{\prime}_{n-1}$

[TABLE]

satisfied for $m_{l}\geq m_{n}$ only if $\Delta_{l}=0$ ($m^{\prime}_{l}=m_{l}$) for all $l=1,\ldots,n-1$.

Now let us estimate the number of computations in Eq. (5). There are $N$ multiplications of sums of matrix elements in the rows $1,\ldots,N$. The number of additions in the inner sum (over $l=1,\ldots,n$) can be reduced, similarly as in Refs. [9, 42], by ordering the factors $x_{l_{1}}\ldots 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. (5). Therefore, the total number of computations $\mathcal{C}_{\mathbf{m}}$ in Eq. (5) is defined solely by the outer sum and the product and reads

[TABLE]

Let us note that Eq. (6) correctly estimates the number of computations $C_{\mathbf{m}}=O(N)$ in the case of maximally bunched output $\mathbf{m}=(N,0,\ldots,0)$ as well as the correct result for the non-bunched outputs $O(N2^{N})$ by the Ryser-Glynn algorithm [8, 9], a feature not present in any of the previous estimates [49, 50]. Nevertheless, an interesting observation is that the reduced estimate still obeys the majorization pattern pointed out in Ref. [53].

On average, in a Haar-random multiport $U$, the speedup in the number of computations due to Eq. (5) (due to the denominator in Eq. (6)) as compared to Eq. (3) is at most by $O(\ln N)$. To show that, consider the probability of maximal bunching of bosons at the output of a Haar-random network, given as [2]

[TABLE]

Setting in Eq. (7) the probability to be $1-\epsilon$ we obtain

[TABLE]

3 Classical computations for calculating a single probability of the output distribution of many-boson interference

To quantify the number of classical computations for general $M\geq N$, when an arbitrary output configuration $\mathbf{m}=(m_{1},\ldots,m_{M})$ occurs with a non-vanishing probability (with the uniform average probability in a Haar-random network), we study how the total number of output ports occupied by bosons is distributed. It is shown below that, for $N\gg 1$ in a Haar-random multiport the total number of output ports occupied by bosons has a bell-shaped distribution, with the tails bounded by a binomial distribution. This allows one to get an almost sure lower and upper bounds on the computational complexity for any network, i.e., the bounds apply with arbitrarily close to $1$ probability to a randomly chosen network according to the Haar measure.

Given the number $1\leq n\leq N$ of output ports occupied by bosons, the number of classical computations necessary to compute an output probability $\hat{p}(\mathbf{m})$, according to the algorithm in Eq. (5), satisfies

[TABLE]

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,\ldots,n-1$, whereas the upper bound by uniformly distributing bosons over the occupied output ports in Eq. (6), i.e., $m_{l}=N/n$ for $l=1,\ldots,n$ (an upper bound for $m_{l}=[N/n]$ for $\leq n-1$ and $m_{n}=N-(n-1)[N/n]$ in Eq. (6)).

Consider now the probability distribution of the number of output ports occupied by bosons in a Haar-random network. Thanks to the uniform average probability $\langle\hat{p}(\mathbf{m})\rangle=\binom{M+N-1}{M-1}^{-1}$ the probability of $n$ output ports occupied by bosons is obviously (here $\theta(x)=1$ for $x>0$ and $\theta(x)=0$ otherwise)

[TABLE]

Indeed, we choose uniformly randomly $n$ out of $M$ output ports and distribute in them $N-n$ bosons with the uniform probability over the output configurations. The average number $\langle n\rangle$ of the output ports occupied by bosons reads

[TABLE]

As $N\gg 1$ (for a fixed $\rho=N/M$) the tails of the distribution (10) (the precise definition of the tails is given below) lie below the tails of the following binomial distribution

[TABLE]

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

[TABLE]

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. (11). Moreover, by using Stirling’s approximation for the factorials [54], one can easily check that for $n=N/(1+\rho)$

[TABLE]

i.e., the maximum of the binomial distribution lies 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 point $n=\overline{n}$ the binomial distribution dominates that of the total number of occupied ports.

Let us formally define the left $1\leq n\leq n_{-}$ and right $n_{+}\leq n\leq N$ tails by the points $n_{\mp}$ of equal probability $P(n_{\pm})=B_{n_{\pm}}(x)$. Now, let us show that such points do exist and are given by the expression $n_{\pm}=\frac{1\pm\delta_{\pm}}{1+\rho}N$, for some $\delta_{\pm}>0$ when $N\gg 1$. To find the points $n_{\pm}$ consider the equation

[TABLE]

Applying the standard approximation to factorials [54] in the binomial $\binom{N-1}{n-1}$, after simple algebra we obtain the following asymptotic equation for $\delta_{\pm}$

[TABLE]

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

[TABLE]

For $0<\rho<1$ the r.h.s of Eq. (13) varies between [math] and $\ln 2$, i.e., the minimum and maximum values of $H(z)$ for $0\leq z\leq 1$, thus there is always a solution for $\delta_{\pm}>0$ (satisfying also $\delta_{-}<1$ and $\delta_{+}<\rho$). Therefore, we have shown that there are such $n_{\pm}=\frac{1\pm\delta}{1+\rho}N$, with some $\delta=\mathrm{max}(\delta_{\pm})\in(0,\rho)$, that for $1\leq n\leq n_{-}$ and for $n_{+}\leq n\leq N$ the binomial distribution of Eq. (12) dominates the distribution of Eq. (10).

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

[TABLE]

Now, setting $\exp\left(-\frac{\delta^{2}N}{4(1+\rho)}\right)=\epsilon/2$ we get $\delta=2\sqrt{\frac{1+\rho}{N}\ln(\frac{2}{\epsilon})}$, thus for sufficiently large $N\gg 1$ such $\delta$ will satisfy Eq. (13) (and $\delta<\rho$, for a finite-density regime $\rho=\Omega(1)$). Therefore, using the lower and upper bounds $n=n_{\pm}=\frac{1\pm\delta}{1+\rho}N$ in Eq. (9) we have the following theorem.

Theorem 1

For any $\epsilon>0$, with the probability at least $1-\epsilon$ in the Haar measure over the unitary multiports, 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 $N$ indistinguishable bosons on a unitary multiport with $M\geq N$ input and output ports satisfies for $N\gg 1$ the following bounds:

[TABLE]

In formulating theorem 1, by using the max$(\ldots)$ in the definition of $r$ we have taken into account also the vanishing density case, and more generally, the case of $\delta>\rho$, as $N\to\infty$, thus recovering the previous estimate $O(N2^{N})$ (according to Eq. (13) the distribution in Eq. (10) does not possess the right tail in this case). Moreover, by continuity the validity is extended to the case of $\rho=1$.

4 Classical complexity of sampling from $N$-boson quantum interference on a unitary $M$-port

In Ref. [42] the estimate $\mathcal{C}=O(N2^{N}+N^{2}M)$ on the number of classical computations $\mathcal{C}$ 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. [42] applies uniformly over all possible output configurations, since multi-sets of output ports $l_{1}\leq\ldots\leq l_{N}$ are obtained by consecutive sampling from a conditional probability chain rule used there, each time sampling for one output port. The leading order in the above estimate is just the number of computations for a single output probability in the distribution. The estimate of Ref. [42], however, is obtained by using Glynn’s formula [9], which disregards the speedup found in section 2, leading to the lower and upper bounds given by theorem 1. Below we recall the main steps in the algorithm of Ref. [42] and apply it to derive the sampling complexity with account of the speed up of section 3.

The algorithm of sampling of Ref. [42] uses symmetry of the output probability Eq. (1) with respect to permutations of input and of output ports to enlarge the space of events, given by occupations of the output ports $\mathbf{m}=(m_{1},\ldots,m_{M})$, $m_{1}+\ldots+m_{M}=N$, to the space of independent output ports $l_{1},\ldots,l_{N}$, with $1\leq l_{j}\leq M$. This is achieved by using the following summation identity valid for any symmetric function $f(l_{1},\ldots,l_{N})$

[TABLE]

where the summation on the left hand side runs over all $\mathbf{m}$ such that $m_{1}+\ldots+m_{M}=N$. Therefore, by Eq. (16) the probability $p(l_{1},\ldots,l_{N})$ of an ordered sequence of output ports $l_{1},\ldots,l_{N}$ becomes

[TABLE]

Let us now introduce the marginal probability of the first $K$ output ports to be $l_{1},\ldots,l_{K}$ by performing the summation over arbitrary $l_{K+1},\ldots,l_{N}$:

[TABLE]

Observe that, due to the delta-symbols, the permutations on the right hand side of Eq. (18) must act as follows $\sigma(1,\ldots,N)=(k_{1},\ldots,k_{K},k_{K+1},\ldots,k_{N})$ and $\tau(1,\ldots,N)=(k^{\prime}_{1},\ldots,k^{\prime}_{K},k_{K+1},\ldots,k_{N})$, where $(k^{\prime}_{1},\ldots,k^{\prime}_{K})$ is a permutation of $(k_{1},\ldots,k_{K})$. Introducing a permutation $\pi\in S_{N}$, which reorders the input ports in such a way that the first $K$ coincide with $k_{1},\ldots,k_{K}$, expanding the permutations $\sigma=(\sigma^{(I)}\otimes\sigma^{(II)})\pi$ and $\tau=(\tau^{(I)}\otimes\sigma^{(II)})\pi$, with $\sigma^{(I)},\tau^{(I)}$ acting on the first $K$ and $\sigma^{(II)}$ on the last $N-K$ input ports, and performing the summation over $\sigma^{(II)}$ in Eq. (18) we obtain the following expression

[TABLE]

In deriving Eq. (19) we have taken into account that for each subset $k_{1},\ldots,k_{K}$ of $K$ input ports there are $K!(N-K)!$ permutations $\pi$ which select the first $K$ input ports from this subset, therefore summation over $\pi$ is accompanied by the factor $K!(N-K)!$ in the denominator with $(N-K)!$ cancelling the same factor in the numerator due to summation over $\sigma^{(II)}$. Now, by comparing the right hand side of Eqs. (17) and (19) one can see that they have a similar form, except the summation over $\pi$ in Eq. (19). The latter summation can be understood as averaging over all possible permutations $\pi\in S_{N}$, if we assume the uniform probability $p(\pi)=1/N!$. Under this condition, the marginal probability becomes

[TABLE]

One can easily verify that our definition of the conditional probability $p(l_{1},\ldots,l_{K}|\pi)$ in Eq. (20) is self-consistent, i.e., taking a marginal probability of another (larger) marginal probability results in the marginal probability in the form in Eq. (20). The probability in Eq. (18) can be also put in this form due to the symmetry of the output probability, mentioned above. We have

[TABLE]

where there are exactly $N!$ coinciding terms $p(l_{1},\ldots,l_{N}|\pi)=\frac{1}{N!}\left|\mathrm{per}(U[\pi(1),\ldots,\pi(N)|l_{1},\ldots,l_{N}])\right|^{2}$.

The sampling algorithm of Ref. [42] is based on the chain rule for the conditional probability in Eq. (20):

[TABLE]

Sampling for output ports $l_{1},\ldots,l_{N}$, where $N$ input bosons end up, is performed as follows.

Sample for permutation $\pi\in S_{N}$ with the uniform probability $p(\pi)=1/N!$. This requires $O(N^{2})$ computations (e.g., sample uniformly from $1,\ldots,N$ without replacement, which is equivalent to uniformly sampling the permutations by the probability chain rule). 2. 2.

Given $\pi$ as above, sample for the first output port $l_{1}$ with the conditional probability $p(l_{1}|\pi)$ as in Eq. (20). This requires $O(N)$ computations. 3. 3.

Given the output ports $l_{1},\ldots,l_{K-1}$, sample for the next output port $l_{K}$ by using the conditional probability $p(l_{K}|l_{1},\ldots l_{K-1};\pi)$ of Eq. (22). Since the previously obtained output ports are known at this step, the sampling can be achieved by using only the numerator in Eq. (22), i.e., the probability $p(l_{1},\ldots,l_{K}|\pi)$ for a given set of $l_{1},\ldots,l_{K-1}$. Such probabilities (for $K=2,\ldots,N$) are computed sequentially by using the modified Glynn formula of section 2 and the Laplace expansion of the matrix permanent over unknown output port $l_{K}$:

[TABLE]

Now, let us show that the bounds of theorem 1 apply to the leading order of the number of computations required to produce a single sample in the above algorithm, if we take into account Eq. (5). As above, let $K$ stand for the number of sequentially sampled output ports, whereas let $s$ be the corresponding number of different output ports occupied by $K$ bosons (i.e., $1\leq s\leq n$, where $n$ is the total number of output ports occupied by $N$ bosons, as in section 3). The conditional probabilities are computed by employing the Laplace expansion of Eq. (23) in the modified Glynn’s formula of Eq. (5), by following similar steps as in lemma 2 of Ref. [42].

Consider first the lower bound on computations given by Eq. (9). By Eq. (23) the total number of computations in steps 1 – 3 above for cumulative sequences of output ports $l_{1},\ldots,l_{K}$, where $1\leq K\leq N$, becomes

[TABLE]

where $\sum\limits_{K\geq s}K$ is the sum over the sequences with $s$ different occupied ports. The r.h.s. of Eq. (24) is obviously bounded from below by the sequence of sampled ports where the first $N-n+1$ output ports are the same, whereas after this the sampled sequence gets each time a new output port. The first $(N-n+1)M$ computations of permanents of rank $1$ for different output ports $1\leq l\leq M$ are dominated by the term $\sum_{K}MK=MN^{2}$ in Eq. (24). We have therefore $\mathcal{C}\geq\mathcal{C}_{\mathrm{lb}}$ with

[TABLE]

which reduces to the estimate found in Ref. [42] in the case of vanishing density $\rho\ll 1/N$ (when $n=N$).

Now consider the upper bound on the number of computations $\mathcal{C}^{(K)}_{\mathbf{m}}$ at each step $K=m_{1}+\ldots+m_{M}$. To derive an upper bound we can drop the denominator in Eq. (6) (see also the discussion at the end of section 2) to get a simpler upper bound on the number of computations $\mathcal{C}^{(K)}_{\mathbf{m}}=O\left(N\prod_{l=1}^{M}(m_{l}+1)\right)$. Then the number of computations $\mathcal{C}^{(K)}_{\mathbf{m}}$ for the steps $K$ and $K+1$ are related as follows (except a constant factor, omitted for simplicity)

[TABLE]

where $l$ is the output port to which a boson is added in step $K+1$ and $m$ is the maximal bunching at an output port, see Eq. (8) of section 2. Therefore, with probability $1-\epsilon$, the total number of computations satisfies $\mathcal{C}\leq\mathcal{C}_{\mathrm{ub}}$, where

[TABLE]

and $m$ is given by Eq. (8). Recalling theorem 1 of section 3, we have the following result.

Theorem 2

For any $\epsilon>0$, with the probability at least $1-\epsilon$ in the Haar measure over the unitary multiports, the number of classical computations $\mathcal{C}$ required to produce a single sample from the output distribution of quantum interference of $N$ indistinguishable bosons on a unitary multiport with $M\geq N$ input and output ports satisfies for $N\gg 1$ the following bounds:

[TABLE]

Note that, for $N\gg 1$ and a fixed $\epsilon$, $r$ is a growing function of $\rho$, whereas $(1+r)^{\frac{1}{r}}$ is monotonously decreasing with the lower bound for $r=2$ (for $\rho=1$ and $\delta=0$) being $\sqrt{3}\approx 1.73$. In its turn, $2^{\frac{1-\delta}{1+\rho}}$ is also decreasing with $\rho$ with the minimum (for $\rho=1$ and $\delta=0$) $\sqrt{2}\approx 1.41$.

The upper and lower bounds of theorem 2 reduce to the estimate of Ref. [42], $\mathcal{C}=O(N2^{N}+N^{2}M)$, in the no-collision (vanishing density) regime, i.e., when $\rho\ll 1/N$.

From the lower bound in Eq. (2) one can conclude that for $N\gg 1$ the $N$-boson quantum interference in the finite-density regime is at least as hard to simulate classically as that in the no-collision regime with $\mathcal{N}=N/(1+\rho)$ bosons.

5 Open problems

In the present work we have focused on single bosons at the input of a unitary multiport, however, one can adopt the results for arbitrary Fock states at the input. We are bound by the modified Glynn formula Eq. (5) which can use either repeated columns or rows to obtain a speed up over the usual Glynn formula [9]. For an input Fock state $\mathbf{s}=(s_{1},\ldots,s_{k})$, $s_{1}+\ldots+s_{k}=N$, one can swap the input ports and the output ports (i.e., rows and columns) in Eq. (3) to obtain a similar expansion where now multiple occupations of the input ports$,s_{1},\ldots,s_{k}$ are used instead $m_{1},\ldots,m_{n}$ of the output ports. In this way, one obtains an estimate on the number of computations similar to that of Eq. (6) with $s_{1},\ldots,s_{k}$ replacing $m_{1},\ldots,m_{n}$. Hence, the minimum runtime for either column-based or row-based expansion reads

[TABLE]

It is left as an open question if there is a more advanced generalisation of Glynn method, or that of Eq. (5), that requires smaller number of computations in the case of Fock states at the input than the estimate given by Eq. (29). A related problem is to find the speed up of the number of computations to produce a single sample due to Fock state at the input (multiple occupations of the input ports). Such a speedup can be possible already with the estimate of Eq. (29). Indeed, in the algorithm of section 4 the permanents of sizes $1\leq K\leq N$ are sequentially computed. For $K\geq 2$ in each such computation one could either employ the row or column based expansion as in Eq. (3). Therefore, at each size $K$ the computation could take advantage over the best speedup due to either repeated rows or columns. This problem is left for future work.

Another open question is direct estimate of the tails of the probability distribution of output ports occupied by bosons, Eq. (10), instead of using the binomial distribution which results in sufficient but not necessary bounds on the tails of the distribution. With the direct estimate of the tails, tighter bounds on the number of computations (i.e., smaller gaps between lower and upper bound), than those in theorems 1 and 2, could be obtained. Moreover, this would allow to generalise the estimate on the complexity of theorems 1 and 2 of sections 3 and 4 to arbitrary (non-fractional) density of bosons $\rho=N/M$, i.e., for the case of Fock states input with $M<N$. Indeed, to prove the theorems we have used the auxiliary binomial distribution in section 3 which bounds the tails of the distribution of the number of ports occupied by bosons, but the important step given by Eq. (13), which defines the tails, has solutions only for $\rho\leq 1$.

6 Acknowledgements

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.