Classically simulating near-term partially-distinguishable and lossy boson sampling

TL;DR

This paper presents a new classical simulation method for Boson Sampling with partial distinguishability and loss, achieving polynomial runtime improvements over existing approaches in realistic experimental scenarios.

Contribution

The authors develop a classical simulation scheme for Boson Sampling with partial distinguishability and loss, offering significant runtime improvements in practical regimes.

Findings

Polynomial runtime improvement over ideal Boson Sampling simulation.

Enhanced efficiency in simulating realistic, lossy, and partially distinguishable Boson Sampling.

Applicable to near-term quantum experiments, challenging claims of quantum supremacy.

Abstract

Boson Sampling is the problem of sampling from the same distribution as indistinguishable single photons at the output of a linear optical interferometer. It is an example of a non-universal quantum computation which is believed to be feasible in the near term and cannot be simulated on a classical machine. Like all purported demonstrations of "quantum supremacy", this motivates optimizing classical simulation schemes for a realistic model of the problem, in this case Boson Sampling when the implementations experience lost or distinguishable photons. Although current simulation schemes for sufficiently imperfect boson sampling are classically efficient, in principle the polynomial runtime can be infeasibly large. In this work, we develop a scheme for classical simulation of Boson Sampling under uniform distinguishability and loss, based on the idea of sampling from distributions where at most k photons are indistinguishable. We show that asymptotically this scheme can provide a polynomial improvement in the runtime compared to classically simulating idealised Boson Sampling. More significantly, we show that in the regime considered experimentally relevant, our approach gives an substantial improvement in runtime over other classical simulation approaches.

Equations

$$\mathrm{Pr}[S^{\prime }]=\frac{1}{{\prod }_{i=1}^m{S}_i!S_i^{\prime }!}\sum _{\tau ,{\tau }^{\prime }\in {\text{S}}_n}\prod _{k=1}^n{U}_{s_k^{\prime },s_{\tau (k)}}{U}_{s_k^{\prime },s_{{\tau }^{\prime }(k)}}^{\ast }{\mathcal{S}}_{{\tau }^{\prime }(k),\tau (k)},$$

$$\mathrm{Pr}[S^{\prime }]$$

$$=\frac{1}{{\prod }_{i=1}^m{S}_i!S_i^{\prime }!}\sum _{\tau \in {\text{S}}_n}\prod _{k=1}^n{U}_{s_k^{\prime },s_{\tau (k)}}\sum _{{\tau }^{\prime }\in {\text{S}}_n}\prod _{k=1}^n{U}_{s_k^{\prime },s_{{\tau }^{\prime }(k)}}^{\ast }$$

$$=\frac{1}{{\prod }_{i=1}^m{S}_i!S_i^{\prime }!}\mathrm{Per}(U_{S,S^{\prime }})\mathrm{Per}(U_{S,S^{\prime }}^{\ast })$$

$$=\frac{1}{{\prod }_{i=1}^m{S}_i!S_i^{\prime }!}|\mathrm{Per}(U_{S,S^{\prime }}){|}^2,$$

$$\mathrm{Pr}[S^{\prime }]$$

$$=\frac{1}{{\prod }_{i=1}^m{S}_i!S_i^{\prime }!}\sum _{\tau \in {\text{S}}_n}\prod _{k=1}^n|U_{s_k^{\prime },s_{\tau (k)}}{|}^2$$

$$=\frac{1}{{\prod }_{i=1}^m{S}_i!S_i^{\prime }!}\mathrm{Per}(|U_{S,S^{\prime }}{|}^2),$$

$$\rho =\frac{1}{n!{\prod }_{i=1}^m{S}_i!}\sum _{\sigma ,{\sigma }^{\prime }\in {\text{S}}_n}\sigma |s⟩⟨s|{\sigma }^{\prime †}\prod _{k=1}^n{\mathcal{S}}_{{\sigma }^{\prime -1}(k),{\sigma }^{-1}(k)},$$

$$\mathrm{Pr}[S^{\prime }]=\sum _{\sigma \in S_n}\prod _{i=1}^n{\mathcal{S}}_{i,\sigma (i)}\mathrm{Per}(M\ast M_{1,\sigma }^{\ast }),$$

$$\mathrm{Pr}[S^{\prime }]=\sum _{j=0}^n\sum _{{\sigma }^j}{x}^j\mathrm{Per}(M\ast M_{1,\sigma }^{\ast }),$$

$$\mathrm{Pr}[S^{\prime }]=\sum _{j=0}^n\sum _{{\sigma }^j}{x}^j\sum _{\begin{array}{c}{J}^{\prime }\le S^{\prime }\\ |J^{\prime }|=j\end{array}}\mathrm{Per}(M_{J^{\prime },1}\ast M_{J^{\prime },{\sigma }_p}^{\ast })\mathrm{Per}(|M_{\overline{J^{\prime }},{\sigma }_u}{|}^2).$$

$$c_j=\sum _{\begin{array}{c}{J}^{\prime }\le S^{\prime }\\ |J^{\prime }|=j\end{array}}\mathrm{Per}(M_{J^{\prime },1}\ast M_{J^{\prime },{\sigma }_p}^{\ast })\mathrm{Per}(|M_{\overline{J^{\prime }},{\sigma }_u}{|}^2).$$

$$\mathrm{Var}[\mathrm{Per}(M_{J^{\prime },1}\ast M_{J^{\prime },{\sigma }_p}^{\ast })]=\frac{j!}{m^{2j}},$$

$$\mathrm{Var}[\mathrm{Per}(|M_{\overline{J^{\prime }},{\sigma }_u}{|}^2)]<\frac{(n-j)!}{m^{2(n-j)}}\sum _{l=0}^{n-j}\frac{1}{l!}.$$

$$\mathrm{Var}[c_j]$$

$$\to {\left(\frac{n!}{m^n}\right)}^2\text{ as }n\to \infty ,$$

$$\mathrm{Cov}[c_j,c_j^{\prime }]=0\forall j\ne j^{\prime }.$$

$$\mathrm{Var}[\Delta \mathrm{Pr}[S^{\prime }]]$$

$$=\mathrm{Var}\left[\sum _{j=k+1}^n{x}^j\mathrm{Var}[c_j]\right]$$

$$<{\left(\frac{n!}{m^n}\right)}^2\left(\frac{x^{2(k+1)}}{1-x^2}\right),$$

$${\rho }_{n,x}=\frac{1}{n!}\left(\sum _{\sigma ,{\sigma }^{\prime }\in S_n}\sigma |s⟩⟨s|{\sigma }^{\prime }{x}^{\sigma \cdot {\sigma }^{\prime }}\right),$$

$${\rho }_{n,x}=\frac{1}{n!}\left(\sum _{i=0}^n{x}^i\sum _{\begin{array}{c}\sigma ,{\sigma }^{\prime }\in S_n\\ \exists I\subseteq [n]\text{ s.t. }|I|=i\\ {\sigma }^{-1}(j)\ne {\sigma }^{\prime -1}(j)\forall j\in I\\ {\sigma }^{-1}(j)={\sigma }^{\prime -1}(j)\forall j\notin I\end{array}}\sigma |s⟩⟨s|{\sigma }^{\prime †}\right).$$

$${\rho }_{n,x}$$

$$=\sum _{i=0}^n{p}_i\sum _{\begin{array}{c}I\subseteq [n]\\ |I|=i\end{array}}\left(\frac{1}{n!}\sum _{\begin{array}{c}\sigma ,{\sigma }^{\prime }\in S_n\\ {\sigma }^{-1}(j)={\sigma }^{\prime -1}(j)\forall j\notin I\end{array}}\sigma |s⟩⟨s|{\sigma }^{\prime †}\right)$$

$$=\sum _{i=0}^n{p}_i\sum _{\begin{array}{c}I\subseteq [n]\\ |I|=i\end{array}}{\rho }_I,$$

$$p_i=x^i(1-x{)}^{n-i}.$$

$$\frac{1}{n!}\left(\sum _{i^{\prime }=i}^n{x}^{i^{\prime }}(1-x{)}^{n-i^{\prime }}(\genfrac{}{}{0ex}{}{n-i}{i^{\prime }-i})\right)$$

$$=\frac{1}{n!}\left(x^i\sum _{i^{\prime }=i}^n{x}^{i^{\prime }-i}(1-x{)}^{n-i^{\prime }}(\genfrac{}{}{0ex}{}{n-i}{i^{\prime }-i})\right)$$

$$=\frac{1}{n!}\left(x^i\sum _{i^{\prime }=0}^{n-i}{x}^{i^{\prime }}(1-x{)}^{n-i-i^{\prime }}(\genfrac{}{}{0ex}{}{n-i}{i^{\prime }})\right)$$

$$=\frac{x^i(x+1-x{)}^{n-i}}{n!}$$

$$=\frac{x^i}{n!}.$$

$$\mathrm{Tr}[{\rho }_{n,x}]$$

$$=\sum _{i=0}^n{x}^i(1-x{)}^{n-i}\sum _{\begin{array}{c}I\subseteq [n]\\ |I|=i\end{array}}1$$

$$=\sum _{i=0}^n{x}^i(1-x{)}^{n-i}(\genfrac{}{}{0ex}{}{n}{i})$$

$$=(x+1-x{)}^n$$

$$=1.$$

$$E[\Delta P]=\sum _{S^{\prime }}E|\mathrm{Pr}[S^{\prime }]-{\mathrm{Pr}}_k[S^{\prime }]|,$$

$$E\left[|\mathrm{Pr}[S^{\prime }]-{\mathrm{Pr}}_k[S^{\prime }]|\right]=$$

Full text

I Introduction

Since the idea of quantum computers was first proposed, speedups have been theoretically demonstrated for a variety of problems Montanaro (2016). More recently, this has turned to interest in near term quantum speedups on realistic hardware Harrow and Montanaro (2017). Boson Sampling is an example of such a problem Aaronson and Arkhipov (2011), which considers the complexity of sampling from the same probability distribution as $n$ collision-free indistinguishable photons output from an $m$-mode linear optical interferometer drawn uniformly at random. Aaronson and Arkhipov showed that the ability to exactly simulate Boson Sampling in polynomial time would imply that $P^{\#P}=BPP^{NP}$, and the Polynomial Hierarchy would collapse to the third level. This is because the output probabilities of a linear optical interferometer are proportional to the permanent of an $n\times n$ complex matrix Scheel and Buhmann (2008), which is $\#P$-Hard to compute. The same result was also shown for approximately sampling from the same distribution, modulo two conjectures related to the permanents of Gaussian matrices Aaronson and Arkhipov (2011). This led to significant interest in Boson Sampling, with experimental demonstrations ranging from two to six photons Broome et al. (2013); Spring et al. (2013); Tillmann et al. (2013); Crespi et al. (2013); Carolan et al. (2015); Wang et al. (2017); Zhong et al. (2018); Paesani et al. (2018).

This has driven interest in improving classical simulations of Boson Sampling, in order to verify any speedup gained from an experimental implementation. Neville et al. Neville et al. (2017) showed that Boson Sampling can be classically performed for 30 bosons across 900 modes at a rate of 250 samples in under five hours on a standard laptop. This was further improved by Clifford and Clifford Clifford and Clifford (2018), who showed that sampling from $n$ photons across $m$ modes can be classically performed in $O(n2^{n}+\operatorname{\textrm{Poly}}(n,m))$ time and $O(m)$ space.

Recently classical simulations have been expanded to consider practical issues such as photon distinguishability, based on a rich collection of theoretical work in this area de Guise et al. (2014); Tamma (2014); Shchesnovich (2014); Rohde (2015); Shchesnovich (2015); Tamma and Laibacher (2015); Laibacher and Tamma (2015); Tichy (2015); Tillmann et al. (2015); Menssen et al. (2017); Laibacher and Tamma (2017); Spagnolo et al. (2014); Giordani et al. (2018). Renema et al. Renema et al. (2018a) demonstrated that Boson Sampling with partially-distinguishable photons can be simulated in time which grows polynomially with $n$. This was later expanded to consider loss as well Renema et al. (2018b). However, the runtime might still not be efficient in practice, as the polynomial can be large. There is also a further disadvantage in that the error bounds are the average case for a random linear optical interferometer, meaning that there could be interferometers for which the algorithm performs significantly worse. A significant improvement could be achieved through adapting the method of Clifford & Clifford to this algorithm, but there are challenges with this approach.

Another motivation for adapting Clifford & Clifford to photonic imperfection is the task of classically simulating other photonic regimes, such as Gaussian Boson Sampling Hamilton et al. (2017). The probability distribution of Gaussian Boson Sampling is that of $n$ indistinguishable squeezed states at the output of an $m$-mode linear optical interferometer, and depends on the Hafnian of a matrix. Unlike Boson Sampling, there is no known polynomial time classical algorithm for computing the outcome of $n$ fully distinguishable squeezed states. On the other hand, it is classically efficient to sample distinguishable squeezed states, via a similar approach to that used for classically sampling distinguishable single photons Aaronson and Arkhipov (2014). As a result, adapting the Clifford & Clifford algorithm to non-ideal Boson Sampling models provides a first step towards being able to classically simulate imperfections in Gaussian Boson Sampling.

Here we consider the cost of classically simulating Boson Sampling when the photons are partially distinguishable or lossy. We look at the same model of distinguishability as considered in Renema et al. (2018a, b), and use techniques for modelling photon distinguishability in first quantization Moylett and Turner (2018); Stanisic and Turner (2018) to show that this is akin to choosing the indistinguishable photons of a Boson Sampling experiment via the binomial distribution. We combine this with the well-studied model of uniform loss, which also follows a binomial distribution. This gives rise to a method which is able to naturally apply the Clifford & Clifford algrithm and take advantage of its efficiency. This algorithm also offers a worst-case error bound for any linear optical interferometer, rather than a random one. Although this approach only offers a polynomial improvement compared to the runtime for ideal Boson Sampling (unlike the exponential improvement shown in Renema et al. (2018a, b)) we use analytical bounds to show that for photon numbers of experimental interest our algorithm can make a significant improvement over alternative approaches.

This article is laid out as follows. In Sec. II, we give an overview of distinguishability in first quantization, and summarise the work of classical simulation algorithms including Renema et al. (2018a, b). In Sec. III, we show what the model of distinguishability looks like in first quantization, and provide an alternative classical simulation. In Sec. IV, we consider average error bounds for a Haar-random unitary interferometer, via the methods explained in Sec. II.3. In Sec. V, we improve this bound to a worst-case error bound, by computing an upper bound for the trace distance between our approximation and the model. In Sec. VI, we expand these results to consider uniform loss, and show how distinguishability and loss relate to each other. In Sec. VII, we explore these error bounds for experimentally interesting numbers of photons, and show that there are some cases where our algorithm offers an improvement. Finally, we briefly consider non-uniform loss, where loss is is a function of the number of optical components, and use the methods of García-Patrón et al. (2018); Oszmaniec and Brod (2018) to show that classical simulations with non-uniform loss also become easier when distinguishability is introduced. We conclude with some open research questions in Sec. IX.

II Preliminaries

II.1 Partial distinguishability in first quantization

Studying photonics in first quantization allows an arbitrary linear optical experiment to be expressed in the quantum circuit model  Moylett and Turner (2018); Stanisic and Turner (2018); Oszmaniec and Brod (2018). Bosonic symmetrisation can be implemented efficiently, e.g. through use of the quantum Schur transform Bacon et al. (2007), and the main sources of errors in linear optics can be modelled: loss as erasure of qudits and distinguishability as correlation to extra photonic degrees of freedom, both of which lead to decoherence of the quantum circuit.

We start by considering arbitrary distinguishability of photonic scattering, models of which have been proposed in Tamma and Laibacher (2015); Laibacher and Tamma (2015); Tichy (2015). Let $m,n$ be integers such that $m\in O(n^{2})$. Let $U\in\operatorname{\textrm{U}}(m)$ be an $m$-mode linear optical interferometer, and let $S,S^{\prime}$ be $n$-photon $m$-mode Fock vectors, meaning that $S=(S_{0},\dots,S_{m})$, $S^{\prime}=(S^{\prime}_{0},\dots,S^{\prime}_{m})$, and $\sum_{i=0}^{m}S_{i}=\sum_{i=0}^{m}S^{\prime}_{i}=n$. The probability of seeing output occupation $S^{\prime}$ from input $S$ with interferometer $U$ is

[TABLE]

where $\mathcal{S}_{k,l}$ is an $n\times n$ Gram matrix describing the (in)distinguishability of pairs of photons, $S$ is a Fock array describing the arrangement of input photons, $s$ is the corresponding array of single particle states in first quantisation, likewise for $S^{\prime}$ and $s^{\prime}$, $U$ is a $m\times m$ unitary matrix describing the interferometer, and Sn is the group of permutations of $n$ particles. Of particular note is that when $\mathcal{S}$ is the all-1 matrix, we have

[TABLE]

where $U_{S,S^{\prime}}$ is a matrix defined by taking rows and columns of $U$ according to $S$ and $S^{\prime}$ Scheel and Buhmann (2008). This is the output probability distribution for Boson Sampling with indistinguishable photons. Likewise, we can see that when $\mathcal{S}$ is the identity matrix, we have

[TABLE]

which is the distribution for fully distinguishable photons. In first quantisation, the input state can be written as Moylett and Turner (2018)

[TABLE]

For convenience, we will assume here and throughout the rest of this paper that our photons start in the Fock state $\ket{S}=\ket{1^{n}0^{m-n}}$. We will use the corresponding first quantized states of this form in the following.

II.2 Classical simulation of fully (in)-distinguishable lossless Boson Sampling

Boson Sampling under ideal conditions (lossless indistinguishable single photons) is intractable for sufficiently large $n$. Until recently the only classical simulation method explicitly known was to compute the entire probability distribution before taking a sample, though it was widely believed that more efficient, albeit still exponential time, approaches existed. A brute force method cannot scale, due to both the number of possible outcomes and the complexity of computing even one $n\times n$ complex matrix permanent.

Two major results gave the first explicit classical simulation strategies which were faster than brute-force sampling. The first, by Neville et al. Neville et al. (2017), demonstrated that Boson Sampling experiments with up to 30 photons could be simulated on a single laptop, and suggests that a supercomputer could handle up to 50 photons. This was achieved by starting with the classical distribution of $n$ distinguishable photons, and then using Metropolised Independence Sampling to adapt the distribution to that of ideal Boson Sampling. The second result, by Clifford & Clifford Clifford and Clifford (2018), gave a classical algorithm for exact Boson Sampling and runs in time equivalent to computing two $n\times n$ matrix permanents per sample with a polynomial overhead. This is through a combination of optimizations, particularly computing marginal probabilities and sampling via the chain rule. Our approach here is to make these more efficient techniques applicable to realistic situations with distinguishability and loss.

It is also worth discussing at this point how to classically simulate boson sampling with completely distinguishable photons, for which a polynomial time exists Aaronson and Arkhipov (2014). In this case, there is no photon interference, so photons can be sampled individually. This is done by taking a photon which starts in mode $i$, and sampling output mode $j$ with probability $|U_{j,i}|^{2}$. Repeating for all photons gives us the complete sample in $O(mn)$ time.

II.3 Expanding in terms of fixed points

In Renema et al. (2018a, b), Renema et al. consider a model where inter-photon distinguishability is measured by an inner product of pure states Tamma (2014); de Guise et al. (2014); Shchesnovich (2015); Rohde (2015); Tamma and Laibacher (2015); Tichy (2015). The probability distribution of arbitrarily distinguishable bosons is modelled as

[TABLE]

where $\mathcal{S}$ is the same matrix describing the distinguishability as in the previous section, $M$ is a matrix defined by the rows and columns of our interferometer $U$ selected based on our photon output $S^{\prime}$ and input $S$, $M^{*}_{1,\sigma}$ is the conjugate matrix with the identity permutation applied to rows and permutation $\sigma$ applied to columns, and $*$ denotes element-wise multiplication. They further restrict to a model where the indistinguishability overlap is defined by a single parameter $\mathcal{S}_{i,j}=x+(1-x)\delta_{i,j},x\in[0,1]$. The sum over permutations can be ordered based on how many fixed points a permutation has, giving

[TABLE]

where $\sigma^{j}$ denotes permutations which have $n-j$ elements as fixed points. Each permanent can be broken down via the Laplace expansion into a sum of a complex matrix permanent multiplied by a positive matrix permanent:

[TABLE]

Here we are now choosing submatrices of $M$, with $J^{\prime}$ representing the $\binom{n}{j}$ possible combinations of rows from $M$, $\bar{J^{\prime}}$ representing the remaining rows, and $\sigma_{p}$ and $\sigma_{u}$ representing permuted and unpermuted elements of $\sigma$ respectively. The $J^{\prime}\leq S^{\prime}$ notation is used to indicate that $J^{\prime}$ is a Fock state such that $J_{i}\leq S^{\prime}_{i}\,\forall i\in[m]$.

The classical simulation method used truncates the number of fixed points in a permutation as being at most $k$, with the remainder of the probability treated as an error margin. It is important to note while these approximations are real, they are not necessarily positive. This is due to the truncation, where positive higher order terms which would have corrected the probability to be positive are now missing from the approximation. To correct this, any negative approximations are rounded up to [math]. These probabilities are then used to train a Metropolised Independence Sampler, akin to the technique of Neville et al. (2017). Training this sampler requires approximating a number of probabilities dependent on the underlying distribution, each of which involves computing $O(n^{2k})$ permanents of $k\times k$ complex matrices, and the same number of permanents of $(n-k)\times(n-k)$ matrices with non-negative entries. The permanents of $k\times k$ complex matrices can be computed classically in $O(k2^{k})$ time via Ryser’s algorithm, and the permanents of matrices with non-negative entries can be approximated up to multiplicative error in polynomial time Jerrum et al. (2004); Huber and Law (2008). As long as $k$ is independent of $n$, this means that there is a polynomial runtime.

To work out a suitable value $k$, define coefficients $c_{j}$ as

[TABLE]

Assuming the matrices are Gaussian, the variance of each permanent can be bounded as

[TABLE]

and

[TABLE]

This leads to two key results. The first is that the variance of $c_{j}$ tends towards a constant value:

[TABLE]

and the second is that the covariance for different values of $j$ is zero:

[TABLE]

From this one can approximate the variance of the error as a geometric series, which as $n\rightarrow\infty$ tends towards the inequality

[TABLE]

where $P_{k}$ is the probability distribution when truncated at $j\leq k$.

Finally one can use a Markov inequality to show that if the variance of the error is of the form $(n!/m^{n})^{2}\epsilon^{2}$, the average error of the simulation is at most $\epsilon$ Renema et al. (2018b). Crucially, this value of $\epsilon$ is only dependent on $x$ and $k$ and no longer dependent on $n$. This means that for any value of $x$, one can choose a suitable value of $k$ to achieve a required error $\epsilon$, and run a classical simulation in time polynomial in $n$.

Although this runtime is polynomial in terms of $n$ and can therefore be considered asymptotically efficient, it might not be classically simulable in practice. There are three main contributions to this: First, the algorithm is reliant on Metropolised Independence Sampling, which potentially requires many probabilities to be approximated per sample. Second, approximating each probability requires $O(n^{2k})$ permanents of $k\times k$ matrices, which even for small $k$ could be a large number of permanents. And third, approximating each probability requires $O(n^{2k})$ permanents of $(n-k)\times(n-k)$ matrices with non-negative elements. Although approximating permanents of matrices with non-negative elements can be achieved in polynomial time, classical algorithms still have a runtime ranging from $O((n-k)^{4}\log(n-k))$ to $O((n-k)^{7}\log^{4}(n-k))$, depending on the sparsity of the matrix Huber and Law (2008). These issues are the main points to address in order to achieve a practical classical algorithm for Boson Sampling. Clifford & Clifford could help to alleviate these issues, but there is a challenge due to the fact that the approximation used in Renema et al. does not correspond to a bosonic state. This in turn leads to negative probabilities, which are not clear how to correct for the Clifford & Clifford algorithm.

II.4 Classical simulation of lossy Boson Sampling

Another common source of imperfections in linear optics is that of photon loss, which arises through a number of different means. Indeed, any large-scale demonstration of Boson Sampling is bound to face photon loss, and therefore needs to take such issues into account. Some results have already shown instances where hardness is still retained, such as when only a constant number of photons are lost Aaronson and Brod (2016); Wang et al. (2018).

Neville et al. compared the classical simulation of their approach to a Boson Sampling experiment where any photon loss was considered a rejected experiment Neville et al. (2017). In Renema et al. (2018b), the method described in Sec. II.3 was adapted to consider uniform loss, showing that the same result can be found, with the only difference being that $x$ is now replaced by $\alpha=\sqrt{\eta}x$, where $\eta$ is the probability of each individual photon surviving. Crucially, this result demonstrated that Boson Sampling where a constant fraction of photons were lost can be simulated in $O(\ell^{2k}k2^{k})$, where $\ell$ is the number of photons which survive and $k$ is only dependent on the constant $\ell/n$, distinguishability $x$, and the desired accuracy of the simulation. This can be expanded to classically simulating Boson Sampling under uniform loss by sampling $\ell$ from the binomial distribution before sampling output photons, which offers a runtime of $O(n^{2k}k2^{k})$. Novel classical simulations for Boson Sampling under loss have also been considered by use of classically simulable states such as thermal García-Patrón et al. (2018) or separable Oszmaniec and Brod (2018) states.

There has also been some consideration of how classical simulations can be generalised to non-uniform loss. This usually means photon loss that is dependent on the number of optical components, with each component having transmission probability $\tau$. Classical simulation methods can be generalised to this model by identifying a layer of uniform losses from the circuit, followed a non-uniform lossy circuit which can be simulated classically through the use of additional modes for lost photons García-Patrón et al. (2018); Oszmaniec and Brod (2018). These results showed that Boson Sampling under non-uniform loss can be classically simulated as long as the smallest number of components a photon encounters is logarithmic in $n$.

III Expanding in terms of states

We will now introduce a new expression for partially distinguishable particles. We start by writing the input state of $n$ photons, with pairwise distinguishability parameter $x$ as in the previous section, in first quantization

[TABLE]

where we have used $\sigma\cdot\sigma^{\prime}$ to denote the number of places where permutations $\sigma$ and $\sigma^{\prime}$ match. For reference, the expansion of Renema et al. (2018a, b) is carried out by identifying $\sigma$ and $\sigma^{\prime}$ that match for a fixed set of $i$ points:

[TABLE]

Note here that the sums over permutations do not correspond to physical states. This can be seen by the fact that for $i\neq 0$ this summation has no elements along the diagonal of the density matrix, as $\sigma$ and $\sigma^{\prime}$ need to differ in exactly $i$ places.

We instead look at an alternative expansion, in order to decompose the model into a linear combination of physical states:

[TABLE]

where $\rho_{I}$ is the state where photons in modes $j\in I$ are fully indistinguishable from each other, all other photons are fully distinguishable, and $0\leq p_{i}\leq 1$ is a coefficient dependent on $x$ and $n$ determining the probability of a state with $i$ indistinguishable single photons.

Note that unlike Eq. (23), where permutations must differ in exactly $i$ points, in Eq. (24) we allow permutations to differ in at most $i$ points. This means that elements closer to and along the diagonal of the density matrix are also part of this summation, and this means that each sum over $\sigma,\sigma^{\prime}\in\operatorname{\textrm{S}}_{n}$ forms a valid density matrix.

Already we can see how a classical simulation might work – if we are able to sample $p_{i}$ efficiently, then we can choose $\rho_{I}$ by selecting $i$ photons uniformly at random to be indistinguishable. These $i$ photons can be classically simulated using Clifford & Clifford Clifford and Clifford (2018), while the remaining $n-i$ photons are treated as fully distinguishable photons, each of which can be simulated individually in polynomial time Aaronson and Arkhipov (2014); Neville et al. (2017).

III.1 The $p_{i}$ are binomially distributed

Here, we will show that the coefficients $p_{i}$ follow the binomial distribution

[TABLE]

To see that the matrix elements of Eq. (26) with $p_{i}$ binomially distributed equal those of Eq. (23), consider $\sigma,\sigma^{\prime}$ which differ at points in the set $I$, where $|I|=i$; the coefficient here should be $x^{i}$. Contributing to this element of the density matrix will be the state $\rho_{I}$, as well as other states $\rho_{I^{\prime}}$, where $I\subseteq I^{\prime}$. The number of such sets $I^{\prime}$ is $\binom{n-i}{i^{\prime}-i}$, as it is equivalent to choosing $i^{\prime}$ from $n$ elements when $i$ elements have already been chosen. The corresponding matrix element is

[TABLE]

It is not hard to see that the state is normalised

[TABLE]

Thus this model of fixed pairwise distinguishability can be written as an expansion in terms of valid states, where indistinguishable photons are drawn from a binomial distribution.

III.2 Classical simulation

We can now see explicitly how a simulation for Boson Sampling with distinguishable photons would work. First, we sample an integer $i\in[n]$ according to the Binomial distribution with coefficients $n$ and $x$. Next, we sample a subset $I$ of the photons uniformly at random from the $\binom{n}{i}$ possible subsets of size $i$. These are the indistinguishable photons of our simulation, which we simulate using Clifford and Clifford in $O(i2^{i}+\operatorname{\textrm{Poly}}(i,m))$ time. The remaining $n-i$ photons are considered to be distinguishable. Rather than needing to compute the output probabilities of these photons colletively, which could take between $O(n-i)^{4}\log(n-i)$ and $O(n-i)^{7}\log^{4}(n-i)$ time via permanents of matrices with non-negative entries Huber and Law (2008), we can instead sample each distinguishable photon individually. To do so, we take a distinguishable photon in mode $a$, and compute the probability of this photon being measured in mode $b$ as $|U_{b,a}|^{2}$. Thus we can compute all output probabilities and obtain a sample for a single distinguishable photon in $O(m)$ time, meaning that we can obtain a sample for all $n-i$ distinguishable photons in $O(m(n-i))$ time Aaronson and Arkhipov (2014); Neville et al. (2017).

The run time is dominated by the time taken to sample our indistinguishable photons, which can be as large as $O(n2^{n}+\operatorname{\textrm{Poly}}(n,m))$ if we are unlucky. By truncating our Binomial sampling up to some level $k$, we can simulate Boson Sampling up to some level of error. The extent of this error will be the focus of Secs. IV & V.

IV Average case error

We can use the same strategies used in Renema et al. (2018a, b) to derive an error bound for Boson Sampling via state truncation for a Haar-random interferometer. We shall do this by considering the total variation distance between our approximation and the model for partial distinguishability for a Gaussian matrix. This is given by

[TABLE]

where $\operatorname{\textrm{Pr}}_{k}$ is the probability distribution truncated at $k$ indistinguishable photons via our approximation. For a specific outcome $S^{\prime}$, we can expand the right hand side to

[TABLE]

where $P_{I}$ is now the probability distribution with indistinguishable photons defined by set $I$, and $p_{i}^{\prime}=p_{i}/(\sum_{i=0}^{k}x^{i}(1-x)^{n-i}\binom{n}{i})$ is the normalised version of the $p_{i}$ coefficients defined in Sec. III. Note that $P_{I}$ is akin to the distribution arising from state $\rho_{I}$ in Sec. III. We can use the triangle inequality to bound this value to

[TABLE]

where we have introduced variables $\Delta P_{\leq k}$ and $\Delta P_{>k}$ for convenience. We shall consider the expected values of these terms for a Gaussian matrix separately, starting with the latter. Using the Laplace expansion, we find that

[TABLE]

where $U_{I,J}$ is a matrix defined from our interferometer $U$ by selecting columns according to $I$ and rows according to $J$, and $|U_{\bar{I}\bar{J}}|^{2}$ is a matrix whose elements are the absolute values squared of $U_{\bar{I},\bar{J}}$.

We next need to consider the expected values of the matrix permanents in Eq. (43) for a Haar random unitary. To do this, we shall assume the matrix describing our interferometer is Gaussian. This allows us to assume that each entry of $U_{I,J}$ and $|U_{\bar{I},\bar{J}}|^{2}$ is independent, and that the two matrices are independent of each other. Starting with $\operatorname{\textrm{Per}}(|U_{\bar{I},\bar{J}}|^{2})$, we note that this is a Gaussian matrix of size $(n-i)\times(n-i)$, and each entry is the square of two independent Gaussians, meaning that each entry of $|U_{\bar{I},\bar{J}}|^{2}$ has expected value $1/m$. From this, we can calculate the expected value as

[TABLE]

For $U_{I,J}$, we note that each element of $U_{I,J}$ is an independent Gaussian entry, with mean value [math] due to symmetry, and second order moment $\operatorname{\mathbb{E}}[|U_{i,j}|^{2}]=1/m$. The second order moment for the permanent can then be calculated using the same methods as in Aaronson and Arkhipov (2011):

[TABLE]

Because the two matrices are independent Renema et al. (2018a), we can express the expected value of their product as

[TABLE]

Plugging this into Eq. (43), we find that

[TABLE]

where first we use the fact that since our input and output are both collision-free, there are $\binom{n}{i}$ ways of choosing $I$ and $J$, then apply cancellation, and finally substitute the values of $p_{i}$. This gives our error bound for terms not in our approximation.

Next we shall consider the term $\Delta P_{\leq k}$. First we can note that $p_{i}^{\prime}\geq p_{i}\forall i\leq k$, due to the normalisation of $p_{i}^{\prime}$. Thus we can rewrite this term as

[TABLE]

Using the same techniques for computing the Laplace expansion and calculating the expected value of permanents of Gaussian matrices, we can show that

[TABLE]

Next we expand $p_{i}^{\prime}$:

[TABLE]

and use this expansion as well as the value of $p_{i}$ to calculate

[TABLE]

Finally, we use these values and sum over all collision-free $S^{\prime}$, of which there are $\binom{m}{n}\approx m^{n}/n!$ to bound our total variation distance for a Haar-random unitary as

[TABLE]

It seems like there should be some room for improvement in this bound. In particular, the use of the triangle inequality suggests that there might be more precise approximations of the expected distance.

V Worst case error

We shall now give an improved error bound, using a different technique. This bound will be an improvement in two ways: first by improving the error bound by a factor of two, and second by being a bound for the worst-case error of any linear-optical interferometer, rather than a bound for a Haar-random unitary interferometer. We do so by finding an upper bound for the trace distance between our ideal partially distinguishable state and the approximation that results from truncating at some $k$, the size of the largest indistinguishable set of particles. As the trace distance is an upper bound for any POVM measurement, we know that this will provide an upper bound for the difference in distribution produced by any interferometer.

Denoting our truncated state

[TABLE]

where the denominator is a normalising factor, and similarly

[TABLE]

We can now estimate the trace distance by rewriting $\rho_{n,x}$ as

[TABLE]

which, by convexity, gives

[TABLE]

If we want this bound to be small, we can use techniques like those used for working out the tails of the binomial distribution. For example, it is known that for $nx<k<n$ that

[TABLE]

where $D(k/n||x)$ is the relative entropy between coins with bias $k/n$ and $x$, respectively Arratia and Gordon (1989). Choosing a value of $k=n\alpha$ for $\alpha>x$ will give an error bound of $\textrm{exp}(-nD(\alpha||x))$. Thus, choosing such a value for $k$ would give an error that decreases as $n$ increases, albeit at the cost of needing to increase $k$ linearly with $n$. While this is still an exponential time algorithm, it would offer a polynomial speedup over the Clifford & Clifford method for Boson Sampling with indistinguishable photons.

We can also note that the trace distance is only dependent on the initial states and not the measurement outcomes. As a result, this error bound also applies in the case where the output is not collision free.

VI Incorporating loss

We now consider how to adapt this simulation to Boson Sampling under uniform loss. We shall assume that each photon survives with probability $\eta$.

In Oszmaniec and Brod (2018); Moylett and Turner (2018), it was shown that the initial state for Boson Sampling with a fixed number of lost photons can be represented in the first quantisation as the initial state

[TABLE]

where $s_{L}$ is the state where photons in the subset $L$ of the original input photons have survived. In order to generalise this to uniform loss, we append $n-\ell$ ”lost” photons in an additional spatial mode (single particle state [math]) which isn’t affected by the interferometer:

[TABLE]

Note that in the same way that it doesn’t matter which particles are traced out when initially applying the loss, it similarly doesn’t matter which particles are replaced with the $\ket{0}\bra{0}$ state. Uniform loss matches that of choosing which subset of photons survive according to the binomial distribution Oszmaniec and Brod (2018); Renema et al. (2018b). We can combine this model with the distinguishability model of Sec. III, giving

[TABLE]

We can now see how our classical simulation for Boson Sampling under distinguishability can be adapted to accommodate loss as well. First, we choose a subset of photons $L$ to indicate the photons that were not lost. From this subset, we choose another subset of photons $I\subseteq L$ to indicate the indistinguishable photons, which are simulated via the Clifford & Clifford algorithm. The photons in $L\setminus I$ are all distinguishable photons, and can be simulated classically as before. The classical complexity of this algorithm depends on the number of indistinguishable photons we choose. As in Section III, by truncating this to be some maximum size $k$ we can get an algorithm that runs in $O(k2^{k}+\operatorname{\textrm{Poly}}(k,m))$ time.

To understand the precision of this algorithm, we first note that if $|L|\leq k$, then we can classically simulate any number of indistinguishable photons within our desired runtime. As a result, we only need to truncate when $|L|>k$, and only need to do so up to $|I|\leq k$.

As with Sec. IV, we start by considering the error bound for a random interferometer. We note that in cases where at most $k$ photons survive our approximation is exact, so these outcomes do not contribute to our total variation distance. For the remainder, we see that

[TABLE]

where $\Delta P_{L}$ denotes the error of our simulation with photons in input modes denoted by $L$. Using our bound in Sec. IV as well as the rule of conditional binomial distributions (see below), we can bound this as

[TABLE]

Again, an improvement over the use of the triangle inequality can lead to an improvement in this bound.

For the worst-case error, we construct analogous states to those in Sec. V:

[TABLE]

[TABLE]

and note that $\rho_{n,\eta,x}$ is a linear combination of these states. As a result, the worst-case error of this simulation, using the convex properties of the trace distance, can be bounded as

[TABLE]

where in the second line we have used the rule of conditional binomial distributions. Using the same result as used for Eq. (68), we can bound the error to a constant value if $k=n\alpha$ for $\alpha>\eta x$ Arratia and Gordon (1989). This shows a relationship between distinguishability and loss similarly to, but not exactly the same, as the one found in Renema et al. (2018b): the more distinguishable photons are, the more we can classically simulate photon loss, and vice versa. It is remarkable to see that these two algorithms have a different dependence on $x$ and $\eta$: while state truncation depends on $\eta x$, point truncation depends on $\eta x^{2}$. It is not immediately clear where this difference comes from, and we leave it as an open question.

VII Empirical Errors

A natural question at this point is how to assess the performance of this new approach over that of Renema et al. (2018a, b). It is not immediately clear how to find a fair comparison, as each approach has its own strengths and weaknesses. Truncating based on fixed points has the benefit of the error asymptotically tending towards a constant as $n$ increases, which means that $k$ can be chosen independently of $n$ and does not need to increase. But this comes at the cost of a potentially large, albeit polynomial, runtime of at least $O(n^{2k}k2^{k}(n-k)^{4}\log(n-k))$ Renema et al. (2018a); Huber and Law (2008). Truncating based on states, on the other hand, provides a significant improvement in runtime based on $k$, and is able to run in $O(2^{k}+\operatorname{\textrm{Poly}}(k,m,n))$ time, but at the cost of $k$ increasing linearly with $n$ for constant error.

We therefore consider a variety of comparisons. In Sec. VII.1, we start by considering the highest value of $x$ and $\eta$ simulable by each approach when given the same values of $n$ and $k$. We then introduce the runtime for each algorithm in Sec. VII.2, by comparing how fast they can simulate particular values of $x$ and $\eta$ for increasing $n$. Finally in Sec. VII.3, we compare the highest value of $x$ and $\eta$ simulable by both algorithms for a 90-photon Boson Sampling experiment, where $k$ is varying but under the condition that the algorithms have similar run times. The motivation for this is that $90$ photons has been suggested as strict upper bound for what is achievable using classical computation Dalzell et al. (2018).

Before we go further, we make a few observations on the calculations of error bounds and runtimes used in this section. Rather than using the asymptotic error bounds for fixed point truncation, which assume $n\rightarrow\infty$, we have used bounds for finite $n$. This provides an improvement in the error of up to $1/\sqrt{e}$. For the runtime of fixed point truncation, we explicitly calculate $\sum_{i=0}^{k}\binom{n}{i}R(n,n-i)i2^{i}(n-i)^{4}\log(n-i)$, where $R(n,n-i)=\binom{n}{i}\lceil i!/e\rfloor$ is the number of permutations with $n-i$ fixed points, and we have assumed that computing every permanent of the $(n-i)\times(n-i)$ matrices with non-negative entries requires $O((n-i)^{4}\log(n-i))$ time 111Note that this computation could in the worst case take $O(n-i)^{7}\log^{4}(n-i)$ time, depending on the matrix sparsity Huber and Law (2008).. For Metropolised Independence Sampling, we choose the number of probabilities to approximate via state truncation as 100, which matches currently used “burn in” and “thinning” timesNeville et al. (2017), though it is worth noting that this number could be different depending on the distribution of partially-distinguishable and lossy bosons. For the runtime of state truncation, we use $2k2^{k}+mk(k-1)/2+m(n-k)$, where we have used the fact that the first term is approximately equivalent to computing two matrix permanents, the second term is the polynomial overhead of the Clifford & Clifford approach, and the final term is the polynomial overhead of sampling the fully distinguishable photons Clifford and Clifford (2018). For these calculations, we have also assumed that $m=n^{2}$ 222Note that $m\in O(n^{2})$ is only sufficient to ensure that the probability of seeing collisions from a Haar-random interferometer are small Arkhipov and Kuperberg (2011). The classical hardness of Boson Sampling is also dependent on entries of $U$ being drawn independently with high probability. To ensure this, $m$ could be required to be as large as $n^{6}$ Aaronson and Arkhipov (2011). However, it is widely believed, and often referenced in Boson Sampling experiments, that $m=n^{2}$ should be sufficient Aaronson and Arkhipov (2011); Broome et al. (2013); Spring et al. (2013); Tillmann et al. (2013); Crespi et al. (2013); Carolan et al. (2015); Wang et al. (2017); Zhong et al. (2018)..

Finally, we note that in Secs. VII.1 and VII.3, we only consider distinguishability and loss separately, by comparing the highest value of $x$ simulable in a lossless system, and the highest value of $\eta$ simulated in a full indistinguishable system. Ideally one would compare highest combinations of $\eta$ and $x$ which are classically simulable. However, doing so is complicated by the fact that both methods handle combinations of noise differently: point truncation handles them as the parameter $\eta x^{2}$, whereas state truncation handles them as $\eta x$. With this in mind, we plot values for distinguishability and loss separately, and note that, for the same performance, a reduction in one implies an increase in the other.

VII.1 Comparison at same level of truncation

We start by comparing the performance of the two algorithms when truncated at the same level $k$. This is of interest as in both approaches $k$ is considered to be a parameter defining the interference between photons. To do so, we consider the error bounds of classically simulating $n$-photon Boson Sampling for $n$ ranging between 2 and 100. The values chosen for $k$ depend on $n$: we consider $k=n-1$ as the upper limit of what the two algorithms can achieve without simulating the full distribution, and also $k=n/2$ as a more feasible, though still exponential time, value.

The result is plotted in FIG. 1, where in (1a) we show the highest value of $x$ simulable assuming no loss ($\eta=1$) and in (1b) we show the highest value of $\eta$ simulable assuming the photons are fully indistinguishable ($x=1$). For all cases, we are considering simulations up to 10% error.

There are a number of things we can note from FIG. 1. First is that when $k=n-1$, we can see that both algorithms tend to the same maximum values of distinguishability and loss. In the case of distinguishability, we can easily see why by considering the error bounds of both algorithms. One can see from Eq. (67) that state truncation will have a simple error bound in this case of $\epsilon\leq x^{n}$, meaning that for constant error the largest value of $x$ simulable is $x=\epsilon^{1/n}$. For point truncation, Eq. (20) shows that the error is similarly bounded as $\epsilon\leq x^{n}/\sqrt{e}$, leading to a largest value of $x=(\epsilon\sqrt{e})^{1/n}$. Thus, although the highest value of $x$ simulable via point truncation is higher than that via state truncation, the difference will decrease in the limit of large $n$. Curiously we see the same effect as well in the case of loss, but now the highest value of $\eta$ simulable via state truncation is higher than that of point truncation. Again, this can be shown to hold theoretically: For state truncation the error scales as $\epsilon\leq\eta^{n}$ according to Eq. (78), corresponding to $\eta=\epsilon^{1/n}$; whereas for point truncation we see from Eq. (20) and substituting $x=\sqrt{\eta}$ that the error scales as $\epsilon\leq\eta^{n/2}/\sqrt{e}$, meaning a maximum value of $\eta$ is $\eta=(e\epsilon^{2})^{1/n}$. In the limit of large $n$ these differences will also tail off.

For $k=n/2$, we see that for both distinguishability and loss point truncation is more powerful than state truncation. Although this is harder to formally prove, there is intuition to see why this is the case. For state truncation, we know that for a small error to be achievable we need $k\geq n\eta x$, as this is the mean of the binomial distribution. Thus for $k=n/2$, we have that $\eta x\leq 1/2$, and in both cases we see the highest value of $x$ and $\eta$ tending to a value below $1/2$. For point truncation on the other hand, we know that the error tends to a constant value only dependent on $k$ and $\eta x^{2}$ in the limit of large $n$. As a result, it is unsurprising that for $k$ increasing linearly with $n$ the highest values of $x$ and $\eta$ will increase.

VII.2 Comparison of runtimes

We next consider the runtime required to simulate $n$-photon Boson Sampling up to 10% error via either method. The motivation for this comparison is that the runtime of the two algorithms at the same value of $k$ are significantly different. In particular, the runtime of state truncation is only dependent on $k$ and not $n$, whereas the runtime for point truncation depends on a scaling of approximately $O(n^{2k})$.

To understand how the runtimes scale, in FIG. 2 we plot the runtimes of classically simulating $n$-photon Boson Sampling experiments via the two approaches for fixed values of $\eta$ and $x$. The values of $k$ chosen for each algorithm are the smallest values for an error of at most 10%. For choosing $\eta$ and $x$, we give two example cases. The first (FIG. 2, red), where $\eta=0.755$ and $x=0.975$, is an example of a hypothetical best experiment we could build with current technology, with the most lossless sources (82%) Slussarenko et al. (2017), interferometers (99%) Wang et al. (2018) and detectors (93%) Marsili et al. (2013), and the highest level of photon indistinguishability (97.6%) He et al. (2018). The second case (FIG. 2, green), where $\eta=x=0.5$, is an example of how the two algorithms perform in what would be considered a poor experiment for both distinguishability and loss. Actual Boson Sampling experiments are likely to fall between these two extremes.

In both cases, state truncation appears to outperform point truncation for near-term photon experiments, with point truncation eventually being able to perform faster for larger values of $n$. When $\eta=x=0.5$, point truncation performs better when $n$ is approximately larger than 230. In the case of $\eta=0.755,x=0.976$, point truncation only performs better when $n$ is approximately larger than $390$ photons. This gives an idea of the regions in which the polynomial runtime of point truncation can be better or worse than the exponential runtime of state truncation.

It is also worth noting that just because point truncation is faster than state truncation for large enough $n$ does not necessarily mean that either algorithm is efficient in these cases. When $\eta=x=0.5$, point truncation only becomes more efficient at instances where both algorithms already require the order of $10^{22}$ operations. And in the case where $\eta=0.755,x=0.976$, both algorithms have runtimes on the order of $10^{92}$ operations before point truncation outperforms state truncation.

VII.3 Comparison at same runtime

Now we consider both truncation level and runtime, and compare the algorithms when restricted to comparable runtimes. To do this, we shall consider the challenge of simulating a 90-photon Boson Sampling experiment, and the largest values of distinguishability and loss that can be simulated at level $k$ with 10% error. The motivation for this is that $90$ photons has been suggested as strict upper bound for what is achievable using classical computation Dalzell et al. (2018).

The results are shown in Fig. 3, detailing for each algorithm the highest value of $x$ simulable when $\eta=1$ (3a) and the highest value of $\eta$ simulable when $x=1$ (3b). In both figures, the blue line indicates state truncation at level $k$, and the green line indicates point truncation at level $k$. However, the runtime of state truncation at level $k$ and point truncation at level $k$ are likely to be drastically different. To take runtime into consideration as well, we consider the orange line which indicates point truncation at level $k^{\prime}$, where $k^{\prime}$ is the smallest integer such that the approximated runtime of point truncation at level $k^{\prime}$ is longer than that of state truncation at level $k$. This allows us to compare the performance of the two algorithms when restricted to similar runtimes.

Considering distinguishability in FIG. 3a, we can note that point truncation with comparable runtime performs better up to $k\leq 45$, after which the methods are roughly comparable with state truncation performing marginally better, before becoming more dominant for $k\geq 60$. It has been suggested that boson sampling with 50 indistinguishable photons is roughly the limit of what can be classically simulated on a supercomputer Neville et al. (2017); Clifford and Clifford (2018); Zhang et al. (2018), so it appears that when considering distinguishability, the algorithms are roughly comparable in this case.

When considering loss in FIG. 3b, we see a noticeable improvement for state truncation. Now point truncation under the same runtime only performs better up to $k\leq 22$, with state truncation performing considerably better for higher values of $k$. Boson Sampling with up to 30 indistinguishable photons is already known to be classically simulable on a standard laptop Neville et al. (2017), so this appears to offer a noticeable improvement even for fast classical simulations.

VIII Non-uniform loss

We finish by briefly considering non-uniform loss, where each photon survives a lossy optical component with probability $\tau$. This model of loss has been considered before García-Patrón et al. (2018); Oszmaniec and Brod (2018), but without the incorporation of distinguishability. We can do this using the same methods as other non-uniform loss results, by extracting non-uniform losses into a layer of uniform losses followed by a lossy interferometer. The uniform loss layer means that each photon has probability $\eta=\tau^{s}$ of surviving, where $\tau$ is the loss of each optical component and $s$ is the smallest number of lossy optical components a photon interacts with. If we take the total number of lossy components to be $d$, the remaining lossy circuit can be modelled as an $(m+d)$-mode interferometer, with lost photons ending up in the additional $d$ modes. Thus we can achieve the same error as Eq. (78) in $O(k2^{k}+\operatorname{\textrm{Poly}}(k,m,d))$ time. In typical schemes for linear interferometers, $d$ is at most polynomial in $m$ Reck et al. (1994); Clements et al. (2016), so the overhead from these additional modes is small. We can bound the error to a constant if $k>nx\tau^{s}$. Taking the logarithm on both sides and rearranging for $s$, we find that this holds if

[TABLE]

This matches results in García-Patrón et al. (2018); Oszmaniec and Brod (2018), showing that boson sampling can be classically simulated if each photon encounters at least a logarithmic number of lossy components. It also shows how distinguishability can affect the simulability of lossy components in Boson Sampling: if our photons are more distinguishable, corresponding to a smaller value of $x$, then we can simulate shallower (i.e. less total loss) optical circuits.

IX Conclusion

In recent years significant improvements have been made in the ability of classical computers to simulate Boson Sampling under various imperfections. However, while it is of theoretical interest to demonstrate asymptotic improvements in classical simulation, the whole reason for proposals such as Boson Sampling is to offer speedups for near-term devices. Although our algorithm will not scale polynomially as the number of photons increases, we find that a substantial improvement over current classical algorithms can be achieved for the numbers of photons that experimentalists are currently aiming for. In doing so, we have effectively set a benchmark for what is required of a 50-90 photon Boson Sampling device.

There are a number of ways one could improve this classical simulation. In particular, the approach of Ref. Renema et al. (2018a) for truncation when looking at near-term devices is dependent on Metropolised Independence Sampling. A direct adaptation of the Clifford & Clifford algorithm to this approach would almost certainly offer an improvement over our algorithm. However, such an adaptation is non-trivial, due to the fact that the term in the expansion are not states, something that motivated our work here.

There are other open questions we would like to consider as well. The first would be to improve the average-case error bounds to match, if not improve over, our worst-case error bound. This would most likely involve an alternative to using the triangle inequality. The second would be to find a way of explaining the difference in dependence on $\eta$ and $x$ between point and state truncation, and ideally improving either algorithm in the process.

IX.1 Note Added

During this work we were made aware of independent work by V. Shchesnovich, which also shows that the model of distinguishability considered by Renema et al. corresponds to that of selecting indistinguishable photons via the binomial distribution Shchesnovich (2019a). This is derived using significantly different methods from those used in this manuscript, and does not consider classical simulation of distinguishability via the above method (though this has been anticipated Shchesnovich (2019b)).

No underlying data was produced during this study.

Acknowledgements.

A.E.M. was supported by the Bristol Quantum Engineering Centre for Doctoral Training, EPSRC grant EP/L015730/1. R.G.P. Acknowledges the support of F.R.S.- FNRS and the Fondation Wiener Anspach. J.J.R. acknowledges support from NWO Vici through Pepijn Pinkse. We would like to thank Nicolás Quesada for bringing the point about Gaussian Boson Sampling in Sec. I to our attention.