No imminent quantum supremacy by boson sampling

TL;DR

This paper argues that achieving quantum supremacy through boson sampling is unlikely in the near term, as classical algorithms can simulate larger photon systems than current experiments, and technological improvements are needed.

Contribution

The paper introduces advanced classical algorithms for boson sampling and provides analyses showing the difficulty of demonstrating quantum supremacy with current and near-term technology.

Findings

Classical algorithms can simulate up to 30 photons in boson sampling.

Current experiments with 5 photons are far from demonstrating supremacy.

Photon losses significantly hinder the realization of quantum supremacy.

Abstract

It is predicted that quantum computers will dramatically outperform their conventional counterparts. However, large-scale universal quantum computers are yet to be built. Boson sampling is a rudimentary quantum algorithm tailored to the platform of photons in linear optics, which has sparked interest as a rapid way to demonstrate this quantum supremacy. Photon statistics are governed by intractable matrix functions known as permanents, which suggests that sampling from the distribution obtained by injecting photons into a linear-optical network could be solved more quickly by a photonic experiment than by a classical computer. The contrast between the apparently awesome challenge faced by any classical sampling algorithm and the apparently near-term experimental resources required for a large boson sampling experiment has raised expectations that quantum supremacy by boson sampling is on the horizon. Here we present classical boson sampling algorithms and theoretical analyses of prospects for scaling boson sampling experiments, showing that near-term quantum supremacy via boson sampling is unlikely. While the largest boson sampling experiments reported so far are with 5 photons, our classical algorithm, based on Metropolised independence sampling (MIS), allowed the boson sampling problem to be solved for 30 photons with standard computing hardware. We argue that the impact of experimental photon losses means that demonstrating quantum supremacy by boson sampling would require a step change in technology.

Tables

Table 1. Table 1: p 𝑝 p -value data for various problem sizes and sampler combinations. Here, τ 𝜏 \tau is shorthand for both τ burn subscript 𝜏 burn \tau_{\mathrm{burn}} and τ thin subscript 𝜏 thin \tau_{\mathrm{thin}} .

n	Sampler 1	Sampler 2	2-sample KS test $p$-values
3	MIS	naïve brute-force	0.651	0.34	0.653	0.368	0.089
	MIS	rejection	0.209	0.415	0.679	0.053	0.167
	rejection	naïve brute-force	0.03	0.588	0.39	0.277	0.331
7	MIS	naïve brute-force	0.448	0.606	0.858	0.734	0.468
	MIS	rejection	0.628	0.615	0.246	0.581	0.775
	rejection	naïve brute-force	0.27	0.252	0.287	0.954	0.263
12	MIS	rejection	0.628	0.919	0.998	0.77	0.519
20	MIS ($\tau =100$)	MIS ($\tau =1000$)	0.548	0.258	0.832	0.483	0.585

Equations

$$T\left(x^{\prime }|x\right)=\min \left(1,\frac{f(x^{\prime })}{f(x)}\frac{g(x)}{g(x^{\prime })}\right).$$

$$c_t(n)=a100n{2}^n$$

$$q_t(n,\eta )=\frac{e}{R{\eta }^n}$$

$$\text{QA}(n,\eta)=\max\Big{[}0,\log_{10}\Big{(}\frac{c_{t}}{q_{t}}\Big{)}\Big{]}.$$

$${\text{QS}}_1:\text{QA}>10$$

$${\text{QS}}_2:q_t<1\text{week},c_t>100\text{yrs}.$$

$$\mathrm{Pr}(S)=\frac{{\left|\mathrm{Per}\left(A_S\right)\right|}^2}{s_1!\dots s_n!}$$

$$\mathrm{Per}(X)=\sum _{\sigma \in {\mathcal{S}}_n}\prod _{i=1}^n{x}_{i,\sigma \left(i\right)}$$

$$|{\Psi }_{in}⟩=|\underset{n}{\underbrace{1,\dots ,1}},\underset{m-n}{\underbrace{0,\dots ,0}}⟩$$

$$\sum _{x}g(x)\frac{f(x)}{\lambda g(x)}=\frac{1}{\lambda },$$

$$\tilde{f}(x)=\{\begin{array}{ll}\frac{f(x)}{{\sum }_{x,f(x)\le \tilde{\mu }}f(x)}& \text{if f(x)≤μ}\\ 0& \text{if f(x)>μ}\end{array}.$$

$$\frac{1}{2}\sum _x|f(x)-\tilde{f}(x)|=\sum _{x,f(x)>\tilde{\mu }}f(x).$$

$$T\left(x^{\prime }|x\right)=\min \left(1,\frac{f(x^{\prime })}{f(x)}\frac{g(x)}{g(x^{\prime })}\right)$$

$$\underset{P_{\mathrm{D}}}{\mathrm{Pr}}(S)=\frac{\mathrm{Per}\left(|A_S{|}^2\right)}{s_1!\dots s_n!},$$

$$\frac{P\left(\mathcal{Q}|N_{\mathrm{events}}\right)}{P\left(\mathcal{R}|N_{\mathrm{events}}\right)}=\prod _{x=1}^{N_{\mathrm{events}}}\left(\frac{q_x}{r_x}\right)=\mathcal{X},$$

$$P_{\mathrm{ind}}\equiv P\left(\mathcal{Q}|N_{\mathrm{events}}\right)=\frac{1}{\mathcal{X}+1}\prod _{x=1}^{N_{\mathrm{events}}}\left(\frac{q_x}{r_x}\right)$$

$$P_{\mathrm{CFS}}\approx \left(\genfrac{}{}{0ex}{}{m}{n})/(\genfrac{}{}{0ex}{}{m+n-1}{n}\right).$$

$$q_t\left(n\right)={\left(R(n)P_{\mathrm{CFS}}{\eta }^n\right)}^{-1}$$

$$q_t^{\mathrm{lin}}\left(n\right)$$

$$\sim {\left(R(n){\left(\frac{4{\eta }_{\mathrm{f}}}{5}\right)}^n{{\eta }_0}^{4n^2}\right)}^{-1}$$

$$q_t^{\mathrm{quad}}\left(n\right)$$

$$\sim {\left(R(n){{\eta }_{\mathrm{f}}}^n\left(\frac{1}{e}\right){{\eta }_0}^{n^3}\right)}^{-1}$$

$${\eta }_{\mathrm{f}}={\eta }_{\mathrm{QD}}{\eta }_{\mathrm{de}}{\eta }_{\det },$$

$$c_t(n)=3\times 10^{-13}{n}^2{2}^n$$

$$\text{QA}(n,\eta)=\max\Big{[}0,\log_{10}\Big{(}\frac{c_{t}}{q_{t}}\Big{)}\Big{]}.$$

$${\text{QS}}_1:\text{QA}>10$$

$${\text{QS}}_2:q_t<1\text{week},c_t>100\text{yrs}.$$

$$\mathrm{Pr}(S)=\frac{1}{\left|{\Lambda }_{n,k}\right|}\sum _{T\in {\Lambda }_{n,k}}{\left|\mathrm{Per}\left(A_{S,T}\right)\right|}^2$$

$$L=\left(\begin{array}{cc}\sqrt{\eta }& \sqrt{1-\eta }\\ \sqrt{1-\eta }& -\sqrt{\eta }\end{array}\right)\otimes 1_m$$

$$M_I$$

$$M_O$$

$${\rho }_m=\int \lambda (\mathbf{\alpha },\mathbf{\beta })|\mathbf{\alpha }⟩⟨\mathbf{\beta }|d^m\mathbf{\alpha }{d}^m\mathbf{\beta }$$

$$\mathcal{U}(T)|\mathbf{\alpha }⟩=\underset{i=1}{\overset{m}{⨂}}|\sum _j{T}_{ij}{\alpha }_j⟩=|T\mathbf{\alpha }⟩$$

$$\begin{split}\rho_{O}&=\operatorname{Tr}_{\text{an}}\left[\mathcal{U}(M_{O})\rho\mathcal{U}^{\dagger}(M_{O})\right]\\ &=\operatorname{Tr}_{\text{an}}\left[\int d^{2m}\bm{\alpha}d^{2m}\bm{\beta}\Big{|}\sqrt{\eta}U\bm{\alpha}\Big{\rangle}\!\Big{\langle}\sqrt{\eta}U\bm{\beta}\Big{|}\otimes\big{|}\sqrt{1-\eta}U\bm{\alpha}\big{\rangle}\!\big{\langle}\sqrt{1-\eta}U\bm{\beta}\big{|}\right]\\ &=\int d^{2m}\bm{\alpha}d^{2m}\bm{\beta}\Big{|}\sqrt{\eta}U\bm{\alpha}\Big{\rangle}\!\Big{\langle}\sqrt{\eta}U\bm{\beta}\Big{|}\operatorname{Tr}\left[\big{|}\sqrt{1-\eta}U\bm{\alpha}\big{\rangle}\!\big{\langle}\sqrt{1-\eta}U\bm{\beta}\big{|}\right]\\ &=\int d^{2m}\bm{\alpha}d^{2m}\bm{\beta}\Big{|}\sqrt{\eta}U\bm{\alpha}\Big{\rangle}\!\Big{\langle}\sqrt{\eta}U\bm{\beta}\Big{|}\operatorname{Tr}\left[\mathcal{U}(U)\big{|}\sqrt{1-\eta}\bm{\alpha}\big{\rangle}\!\big{\langle}\sqrt{1-\eta}\bm{\beta}\big{|}\mathcal{U}^{\dagger}(U)\right]\\ &=\int d^{2m}\bm{\alpha}d^{2m}\bm{\beta}\Big{|}\sqrt{\eta}U\bm{\alpha}\Big{\rangle}\!\Big{\langle}\sqrt{\eta}U\bm{\beta}\Big{|}\operatorname{Tr}\left[\big{|}\sqrt{1-\eta}\bm{\alpha}\big{\rangle}\!\big{\langle}\sqrt{1-\eta}\bm{\beta}\big{|}\right]\\ &=\rho_{I}.\end{split}$$

Full text

It is believed that new types of computing machines will be constructed to exploit quantum mechanics for an exponential speed advantage in solving certain problems compared with classical computers Aaronson and Arkhipov (2011); Broome et al. (2013); Tillmann et al. (2013); Spring et al. (2013); Crespi et al. (2013); Lund et al. (2017); Latmiral et al. (2016); Wang et al. (2016a); Shor (1997); Lloyd (1996). Recent large state and private investments in developing quantum technologies have increased interest in this challenge. However, it is not yet experimentally proven that a large computationally useful quantum system can be assembled, and such a task is highly non-trivial given the challenge of overcoming the effects of errors in these systems.

Boson sampling is a simple task which is native to linear optics and has captured the imagination of quantum scientists because it seems possible that the anticipated supremacy of quantum machines could be demonstrated by a near-term experiment. The advent of integrated quantum photonics Politi et al. (2008); O’Brien et al. (2009) has enabled large, complex, stable and programmable optical circuitry Shadbolt et al. (2012); Carolan et al. (2015); Harris et al. (2015), while recent advances in photon generation Ding et al. (2016); Wang et al. (2016b); Somaschi et al. (2016); Spring et al. (2017) and detection Lita et al. (2008); Hadfield (2009) have also been impressive. The possibility to generate many photons, evolve them under a large linear optical unitary transformation, then detect them, seems feasible, so the role of a boson sampling machine as a rudimentary but legitimate computing device is particularly appealing. Compared to a universal digital quantum computer, the resources required for experimental boson sampling appear much less demanding. This approach of designing quantum algorithms to demonstrate computational supremacy with near-term experimental capabilities has inspired a raft of proposals suited to different hardware platforms Bremner et al. (2016a, b); Boixo et al. (2016); Bermejo-Vega et al. (2017).

Based on a simple architecture, the boson sampling problem is similarly straightforward to state. A number $n$ of indistinguishable noninteracting bosons (e.g. photons) should be injected into $n$ input ports of a circuit comprised of a number $m$ of linearly coupled bosonic modes. The circuit should be configured so that the transformation between input and output ports is described by a uniformly (“Haar”) random unitary matrix. The probability for the $n$ bosons to be detected at given set of $n$ output ports is equal to the square of the absolute value of the permanent of the transfer matrix that describes the transformation.

While choosing a number of modes $m\sim n^{5}\log^{2}n$ guarantees that the distribution of any $n\times n$ sub-matrix is approximately equal to that of a matrix of elements drawn independently from the complex normal distribution Aaronson and Arkhipov (2011), the less impractical scaling of $m\sim n^{2}$ is typically targeted. This polynomial relation between $n$ and $m$ is also important because it ensures a not too large probability that two or more of the bosons arrive at the same output port, i.e. bunch; the conjectured hardness only applies to collision-free events, i.e. no bunching. Because approximating the permanent of a random matrix is conjectured to be computationally hard Aaronson and Arkhipov (2011), calculating any transition probability is intractable; the collection of all of the possible collision-free transition probabilities ($m$ choose $n$) constitutes an exponentially large probability distribution, where each element is exponentially hard to calculate. Running an ideal boson sampler would solve the problem of producing samples from this distribution.

Importantly, a strong case for the classical hardness of boson sampling can be made even when the distribution being sampled from is only approximately correct Aaronson and Arkhipov (2011): assuming certain conjectures from computational complexity theory, there can be no efficient classical algorithm to sample from any distribution within a small distance from the boson sampling distribution.

Current estimates for the regime in which photonic experiments could achieve quantum supremacy have been between 20 and 30 photons Aaronson and Arkhipov (2011); Preskill (2012); Bentivegna et al. (2015) and, recently, as low as 7 Latmiral et al. (2016). However, our classical algorithm, based on Metropolised independence sampling (MIS), while necessarily inefficient for large $n$, was able to output a sample for $n=30$ bosons in half an hour on a standard laptop, and would enable an $n=50$ sample to be produced in under 10 days on a supercomputer Wu et al. (2016). MIS Liu (1996, 2004) is a specific Markov Chain Monte Carlo (MCMC) method. For an instance of the problem, our aim is to construct a Markov chain where each state in the chain signifies a boson sampling detection event. New states in the chain are proposed from a classical mockup of the boson sampling distribution: the distribution of distinguishable particles, where probabilities of detection events are equal to permanents of real and positive transition matrices, and sampling becomes efficient Aaronson and Arkhipov (2013). Proposed states are then accepted or rejected by comparison with the actual probability to observe that event for indistinguishable bosons. This means that, at each step, only a relatively small number of permanents must be calculated; a calculation of the full and exponentially large boson sampling distribution is not required.

More precisely, let $P_{\mathrm{D}}$ be the distinguishable particle distribution with probability mass function $g(x)$, over the set of tuples of length $m$ with elements in $\{0,1\}$ which sum to $n$. And let $P_{\mathrm{BS}}$ be the boson sampling distribution over these tuples, with probability mass function $f(x)$. Then starting at a random (according to $P_{\mathrm{D}}$) tuple $x$, propose a new random tuple $x^{\prime}$. The transition from $x$ to $x^{\prime}$ is accepted with probability

[TABLE]

Repeating this procedure generates a Markov chain, which will converge such that thereafter, the chain is sampling from $P_{\mathrm{BS}}$.

Not all states in the Markov chain are retained as detection events. The time taken for the Markov chain to converge means that a number of tuples at the beginning of the chain must be discarded, known as the ‘burn in’ period, $\tau_{\mathrm{burn}}$. For the size of simulations covered here, empirical tests find that a burn in period of 100 is sufficient for convergence to have occurred. In addition, autocorrelation between states in the chain can occur, for example because two consecutive states in the chain will be identical whenever a proposed new state $x^{\prime}$ is not accepted. We empirically find for the size of problem we tackle that autocorrelation is suppressed with a thinning procedure that retains only every $100$th state (see appendix Sec. IV.4). Generally, the burn in period and thinning interval are reduced by a greater overlap between target and proposal distributions, as measured by (1), the transition probability. We find that a proposal distribution of the distinguishable particle transition probabilities has a high acceptance rate of $\sim 40\%$, a sign that the distributions overlap well. In Fig. 1 we present a schematic of the MIS based approach to boson sampling, alongside schematics of a quantum photonic approach and the boson sampling problem itself.

The classical resources used to produce the thinned Markov chain are far fewer than those required by the brute force approach of calculating all permanents in the full distribution. The asymptotically fastest known algorithm for exactly computing the permanent of an $n\times n$ matrix is Ryser’s algorithm Ryser (1963), whose runtime when implemented efficiently is $\mathcal{O}(n2^{n})$. Generating the first tuple in a sample requires the computation of $\tau_{\mathrm{burn}}$ $n\times n$ real valued matrix permanents and $\tau_{\mathrm{burn}}$ $n\times n$ complex valued matrix permanents. Each subsequent sample requires $\tau_{\mathrm{thin}}$ $n\times n$ real valued and complex valued matrix permanents to be computed. The relative scaling of the approaches to classical boson sampling using a standard laptop is shown in Fig. 2(b), setting both $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$ equal to $100$. The MIS sampler is over 50 orders of magnitude faster for the $n=30$ case than the brute force computation of the entire distribution.

We used this algorithm on a standard laptop to produce samples of size 20,000 for up to 20 bosons, and used a local server, which allowed around 30 times more chains to be run in parallel, to produce 250 samples for 30 bosons in 900 modes in less than five hours. As in the experimental case, a central challenge is to provide evidence for sampling from the correct distribution. Here we addressed this using standard statistical techniques. The likelihood ratio tests Bentivegna et al. (2014) in Fig. 2(a) show a rapid growth in confidence in the hypothesis that these samples are from the indistinguishable boson distribution rather than the distinguishable particle distribution for $n=7,12,20$ and $30$.

Further verification results are shown in Fig. 2(c) for the case of 7 bosons in 49 modes. For each tuple in a sample size of 20,000 produced by our classical algorithm, we calculate $-\log(|\text{Per}(A_{S})|^{2})$, where $A_{S}$ is the matrix associated to each tuple, and produce a probability mass histogram. The same function and associated histogram is plotted for a sample of 20,000 tuples chosen from the distinguishable particle distribution (note that $\text{Per}(|A_{S}|^{2})$ gives the probability to observe a transition of distinguishable particles). Fig. 2(c) shows the clear difference between the two distributions, which we analyse with a (bootstrapped) 2-sample Kolmogorov-Smirnov (KS) test Praestgaard (1995); Sekhon (2011). We are able to reject the null hypothesis, that the two samples are chosen from the same distribution, at a significance level of $0.001$. Conversely, a similar KS test between distributions from our MIS algorithm and a rejection sampling algorithm, and between our MIS algorithm and the brute force approach of calculating all permanents in the full distribution, both found large $p$-values (see appendix Sec. IV.4).

We next compare our classical approach with plausible experimental parameters. It is worth noting that asymptotically, experimental boson sampling will have a slower runtime than our algorithm. This is because photon losses scale exponentially badly with $n$ Aaronson and Arkhipov (2011); Rohde et al. (2014). The runtime for an experiment with a transmission probability resulting from fixed loss (generation, coupling, detection) $\eta_{\mathrm{f}}$ and a transmission probability $\eta_{0}$, resulting from loss per unit optical circuit depth scales as \mathcal{O}\big{(}(\frac{1}{\eta_{\mathrm{f}}})^{n}(\frac{1}{\eta_{0}})^{dn}\big{)} for an optical circuit depth of $d$, which is worse than Ryser’s algorithm if $d$ grows with $n$ for any $\eta_{\mathrm{f}},\eta_{0}<1$. However, the region of interest for quantum supremacy is likely to be restricted to $n<100$, where low-loss experiments still have the potential to produce large speedups. Assuming that our MIS sampler continues to perform equally well for larger instance sizes, we can compare its runtime with current and future experiments. The classical and quantum runtimes for an instance of size $n$ bosons in $m=n^{2}$ modes can be estimated as

[TABLE]

where $a$ is the time scaling of the classical computer (for computing one real and one complex permanent), the factor of $e$ is an approximation to the probability of obtaining a collision-free event Arkhipov and Kuperberg (2012), $R$ is the experimental source repetition rate, and $\eta=\eta_{\mathrm{f}}\eta_{0}^{d}$ is the experimental transmission probability of a single photon including the efficiencies of photon generation, coupling, circuit propagation and detection (note that $R$ and $\eta$ will generally be a function of $n$). We define the quantum advantage (QA) as the improvement in quantum runtime versus classical runtime measured in orders of magnitude,

[TABLE]

We now consider two plausible notions of quantum supremacy. First, we can define supremacy as a speedup so large that it is unlikely to be overcome by algorithmic or hardware improvements to the classical sampler, for which we choose a speedup of ten orders of magnitude. Secondly, we may wish to define supremacy as the point at which a computational task is performed in a practical runtime on a quantum device, for which we choose under a week, but in an impractical runtime on a classical device, for which we choose over a century.

These criteria can be summarised as

[TABLE]

In order to make concrete estimates of future runtimes, we need to fix $a$ and $R$. Choosing $a=3n\times 10^{-15}s$ as the time scaling for computing one real and one complex matrix permanent recently reported for the supercomputer Tianhe 2 Wu et al. (2016) and $R=10$GHz, which is faster than any experimentally demonstrated photon source to our knowledge, we can plot QA against $n$ and $\eta$.

We first note that current approaches using spontaneous parametric down conversion (SPDC) photon pairs are generally inefficient with $\eta<0.002$ Broome et al. (2013); Tillmann et al. (2013); Spring et al. (2013); Crespi et al. (2013); Spagnolo et al. (2014); Carolan et al. (2014, 2015); Bentivegna et al. (2015). Recently, improved rates have been demonstrated with quantum dot photon sources Loredo et al. (2017); He et al. (2016); Wang et al. (2016a). The current leading experimental demonstration however is still restricted to $\eta\approx 0.08$ for $n=5$ Wang et al. (2016a) where $q_{t}\approx 10^{9}c_{t}$. This calculation includes the rate used in the experiment ($76n^{-1}\text{MHz}$) and includes a suppression factor caused by a lower collision-free event rate using a linear instead of quadratic mode scaling. In Wang et al. Wang et al. (2016a) a number of realistic, near-term experimental improvements are suggested to reach 20 photon boson sampling. Using these projections we find that $\eta$ is increased to $\approx 0.35$, which would be a major experimental breakthrough. However, as shown in Fig. 2(d), in this case we predict that the classical runtime would still be over six orders of magnitude faster. Fig. 2(e) shows the regions of quantum advantage and quantum supremacy with current and projected experiments.

In Ref. Aaronson and Brod (2016), the authors showed that the boson sampling problem can be modified to allow for a fixed number of lost photons at the input of the circuit whilst retaining computational hardness. In appendix Sec. VI, we show that if the overall losses in the experiment are path-independent then this is equivalent to loss at the input. The MIS sampler can be readily adapted to this scenario by adding an initial step which generates a uniformly random input subset, followed by the usual MIS method for this input state. The dashed contours and lines in Fig. 2(e) and (d) take into account the adjusted classical and quantum runtimes when up to two lost photons are allowed. Although allowing loss helps the experiments to compete, the complexity of realistic experimental regimes such as losing a constant fraction of photons remains unknown and it is easy to see that losing too many photons eventually allows the problem to become efficiently solvable classically.

We conclude that a near term experimental demonstration of quantum supremacy by boson sampling would require a technological step change, reaching photon numbers of over 50 and ultra-low loss interferometers with thousands of modes. A convincing demonstration would additionally need to solve general instances of the problem to a high degree of precision which will require full programmability and increased control, both of which are likely to add loss. Although the boson sampling algorithm could be run on a fault-tolerant quantum computer, this approach would lose the appealing simplicity of the original proposal. However, it may be the case that a more limited, partially error-corrected device could demonstrate quantum supremacy, and theoretical and experimental schemes which can overcome photon loss should be investigated. We finally note that our algorithm is unoptimised and we expect significant improvements can be made, pushing the supremacy threshold further still from the current experimental reality.

I Introduction

The boson sampling problem has, since its introduction by Aaronson and Arkhipov in 2011 Aaronson and Arkhipov (2011), been the subject of great interest within theoretical and experimental quantum physics. Informally, boson sampling is the problem of sampling from the distribution obtained by transmitting photons through a linear-optical network. No efficient classical algorithm for this task is known, but by its nature, boson sampling is suited to being carried out on a specialised linear optical quantum device which could be substantially simpler than a fully universal linear-optical quantum computer. As such, boson sampling is considered a leading candidate for a problem which could see enhanced performance when solved by a quantum device compared with a classical device – so-called quantum supremacy Preskill (2012); Lund et al. (2017) – in the near future.

Estimates for the size of experiment required to demonstrate quantum supremacy via boson sampling have evolved over time. Aaronson and Arkhipov initially predicted this to be in the region of 20-30 photons in a 400-900 mode linear optical network Aaronson and Arkhipov (2011); Preskill suggested “about 30” photons Preskill (2012); and Bentivegna et al. Bentivegna et al. (2015) suggested that, by making some minor modifications to the problem, 30 photons and 100 modes would suffice. Recently Latmiral et al. Latmiral et al. (2016) have reported that the number of photons and modes required to achieve quantum supremacy could be as small as 7 and 50 respectively.

This estimate is based on a comparison with the simplest possible classical technique for simulating boson sampling: evaluating the entire probability distribution. However, there exist well-known classical techniques for exact or approximate sampling that do not require the entire distribution of interest to be determined. Here we assess the requirements on a linear-optical implementation in order to achieve quantum supremacy in light of these.

We find that, using projected loss parameter estimates with a promising experimental boson sampling setup Wang et al. (2016a), quantum supremacy (via the standard boson sampling problem or in a modified version with loss Aaronson and Brod (2016) ) is not achieved.

II The boson sampling problem

The original boson sampling problem, as described by Aaronson and Arkhipov Aaronson and Arkhipov (2011), is based on sampling from the probability mass function defined by the linear scattering of multiple bosons (in practice, photons) prepared in a Fock state and measured in the Fock basis.

Given $U\in SU(m)$, let $A$ be the column-orthonormal, $m\times n$ matrix formed by taking the first $n$ columns of $U$. Also let $\Phi_{m,n}$ be the set of all possible tuples of length $m$ of non-negative integers whose elements sum to $n$.

For some tuple $S=\left(s_{1},\dots,s_{m}\right)\in\Phi_{m,n}$, let $A_{S}$ be the sub-matrix of $A$ with $s_{i}$ copies of row $i$. Boson sampling is the problem of sampling from the probability distribution $P_{\mathrm{BS}}$ over $\Phi_{m,n}$, for a given input $U$, with probabilities defined as:

[TABLE]

where $\mathrm{Per}(X)$ is the permanent of an $n\times n$ matrix $X=\left(x\right)_{ij}$, defined by:

[TABLE]

where $\mathcal{S}_{n}$ is the group of permutations of integers 1 to $n$.

The matrix permanent is similar to the more commonly encountered matrix determinant, but the definition lacks an alternating sign. Although this difference is seemingly minor, computing the permanent of a complex-valued matrix falls in the #P-hard computational complexity class Valiant (1979); Aaronson (2011), and is in general vastly more demanding than computing the determinant. Indeed, the fastest algorithms known for computing the permanent run in time exponential in $n$.

Using the fact that computing the permanent is a #P-hard problem, Aaronson and Arkhipov Aaronson and Arkhipov (2011) have shown that the existence of an efficient, exact classical algorithm for boson sampling would collapse the polynomial hierarchy to its third level, a complexity-theoretic consequence considered extremely unlikely. However, the boson sampling distribution could (in principle) be sampled from using a simple linear-optical circuit with $n$ photons in $m$ modes, where the circuit corresponds to the desired unitary $U$. To resolve the apparent contradiction with #P-hardness, note that their result only shows that a classical sampler from the boson sampling distribution could be used to compute the permanent within the polynomial hierarchy, not that a quantum sampler could.

In addition to these results on the exact form of the problem, Aaronson and Arkhipov proposed two additional (yet plausible) conjectures which together would imply that an approximate form of boson sampling could not be carried out efficiently on a classical computer. In approximate boson sampling, the task is to output a sample from an arbitrary probability distribution $P$ within total variation distance $\epsilon$ of the real boson sampling distribution $P_{\mathrm{BS}}$, for some small constant $\epsilon$. This is a particularly important result for experimental implementations of boson sampling, as it allows for the inevitable imperfections in a real boson sampling device, up to a point Rohde and Ralph (2012); Arkhipov (2015a).

An important feature of the approximate boson sampling problem is that the sub-unitary matrices $A_{S}$ must look like Gaussian random matrices. In order for this to be true, it must hold that $m\gg n$ (the currently best proven bound is $m=\mathcal{O}(n^{5}\log^{2}n)$ Aaronson and Arkhipov (2011), although it is generally accepted that $m=\mathcal{O}(n^{2})$ should be sufficient). A side effect of this condition is that, due to the bosonic birthday paradox Arkhipov and Kuperberg (2012), $\mathrm{Pr}\left(s_{i}>1\right)$ is small, and hence that one can consider the problem restricted to the “collision free subspace” (CFS) of the overall Hilbert space of the bosonic states, where there is no bunching at the output ($s_{i}\leq 1$ for all $i$). So $\Pr(S)=\left|\mathrm{Per}\left(A_{S}\right)\right|^{2}$. Then our task becomes to sample from the restriction of the boson sampling distribution to this subspace, suitably renormalised by a constant close to 1. This is the problem on which we focus in this work.

In the setting of linear optics, an instance of the boson sampling problem for random $U$ can be implemented by injecting the $n$-photon Fock state

[TABLE]

into a linear optical network of beam-splitters and phase-shifters with transfer matrix chosen Haar-randomly, and a single-photon detector coupled to each output mode. These detectors need not be able to discern the number of output photons in each mode, as one can restrict the problem to the CFS. A sample is obtained by recording the $n$-fold coincident detection signals at the output.

Schemes for implementing arbitrary unitary linear optical transfer matrices using beam-splitters and phase-shifters are well known Reck et al. (1994); Clements et al. (2016) and have been realised successfully Carolan et al. (2015); Harris et al. (2015). The current record for implementing the “classic” version of boson sampling is 5 photons in 9 modes Wang et al. (2016a).

For large $n$, producing the input state $\left|\Psi_{in}\right\rangle$ is a major experimental challenge. A “scattershot” variant of boson sampling has therefore been proposed which is easier to implement experimentally, and yet satisfies similar complexity-theoretic properties to classic boson sampling Lund et al. (2014); Aaronson (2013). This approach is useful when using photon pair sources such as spontaneous parametric down conversion (SPDC) sources, where each source has a small probability of generating a photon pair per pulse of a pump laser. In the scattershot setup, $m$ sources are pumped concurrently, and some subset of these sources may generate a photon pair. One photon in each pair is sent to a detector in order to herald the presence of the other, which is sent to an input mode of the linear optical circuit (which may now be any of the $m$ modes, rather than one of the first $n$ modes).

In principle, and with perfectly efficient single photon detectors, this allows one to select the instances where $2n$ detection events occur ($n$ herald signals and $n$ post-circuit signals) and generate samples from a boson sampling distribution with a known input configuration. The benefit of this approach is that, as the number of pairs generated is distributed binomially, one can tune the probability that each source generates a photon pair so as to maximise the probability of an $n$ photon-pair outcome, and this is much greater than if one were limited to $n$ sources. Scattershot boson sampling has been implemented with 6 sources and 13 modes Bentivegna et al. (2015).

From the point of view of classical simulation, the scattershot boson sampling problem is not significantly harder than the classic boson sampling problem. The only difference is that, rather than sampling from a distribution defined by a fixed $m\times n$ submatrix $A$ of an $n\times n$ unitary matrix $U$, one first chooses $A$ itself at random, according to a probability distribution which is easy to sample from classically. We therefore focus on classic boson sampling here, but stress that all our results also apply to scattershot boson sampling.

III Classical sampling techniques for boson sampling

We now describe some techniques that can be used to implement boson sampling on a classical computer. All the classical algorithms we consider will ultimately be based on computing permanents, and our goal will be to minimise the number of permanents computed per sample obtained. The asymptotically fastest known algorithm for exactly computing the permanent of an $n\times n$ matrix is Ryser’s algorithm Ryser (1963), whose runtime when implemented efficiently is $\mathcal{O}(n2^{n})$. The Balasubramanian–Bax/Franklin–Glynn algorithm achieves the same asymptotic performance but may be preferred from the perspective of numerical stability Wu et al. (2016). For the values of $n$ we consider here, numerical stability is unlikely to be a significant issue Wu et al. (2016).

Although this exponential runtime may seem discouraging, note that this is substantially faster than the naïve approach of evaluating the sum (8) directly, which would take time $\mathcal{O}(n!)$. On a standard personal computer, the permanent of an $n\times n$ matrix can be computed in under 1 second for $n\approx 25$ (see Wu et al. (2016); Latmiral et al. (2016) and Sec. IV.1 below for more detailed numerical experiments), and it is predicted in Wu et al. (2016) that a supercomputer could solve the case $n=50$ in under 2 hours. In addition, some proposed verification techniques for boson sampling experiments (e.g. likelihood-ratio tests) rely themselves on computing permanents Aaronson and Arkhipov (2011, 2013).

III.1 Brute force exact sampling

Perhaps the most obvious way to tackle the boson sampling problem on a classical computer is to compute all of the probabilities in the CFS, and sample from the probability mass function associated with these probabilities Latmiral et al. (2016). It becomes apparent that this scheme is immensely computationally demanding when one considers how the dimension of the CFS scales with the size of the problem. Indeed, this approach requires computing $\binom{m}{n}$ permanents of $n\times n$ complex valued matrices before a single sample can be output. Using the lower bound $\binom{m}{n}\geq(m/n)^{n}$, for $m\geq n^{2}$ at least $n^{n}$ permanents must be computed; and even for $n=10$ we need to compute more than $17\times 10^{12}$ permanents. Thus our computation quickly becomes swamped by the number of permanents to compute, rather than the complexity of computing the permanent itself.

We therefore seek an approach which relies on computing a number of permanents which scales more favourably with $n$ and $m$ than $\binom{m}{n}$.

III.2 Rejection sampling with a uniform proposal

Rejection sampling is a general approach for exactly sampling from a desired distribution $P$ with probability mass function $f(x)$, given the ability to sample from a distribution $Q$ with probability mass function $g(x)$, where $f(x)\leq\lambda g(x)$ for some $\lambda$ and all $x$. The algorithm proceeds as follows:

Generate a sample $x$ from $Q$. 2. 2.

With probability $\frac{f(x)}{\lambda g(x)}$, output $x$. Otherwise, go to step 1.

It is easy to show that the sample eventually output by the rejection sampling algorithm is distributed precisely according to $P$. The probability that an accepted sample is generated is

[TABLE]

so the expected number of samples from $Q$ required to output a sample from $P$ is just $\lambda$. The simplest case in which we can apply the algorithm is where $Q$ is the uniform distribution on $N$ elements, and we have the upper bound $f(x)\leq\mu$ for some $\mu$ and all $x$. Then the expected number of uniform samples required to obtain a sample from $P$ is $\mu N$, which will be minimised when $\mu=\max_{x}f(x)$.

Here we take $P$ to be the boson sampling distribution restricted to the CFS, and $Q$ to be the uniform distribution on the CFS (so $N=\binom{m}{n}$). Note that $P$ is subnormalised, so is not quite a probability distribution. However, the rejection sampling algorithm is blind to this subnormalisation (as this is effectively the same as increasing $\lambda$), so will generate samples from the renormalised distribution.

Each iteration of rejection sampling requires the computation of one permanent, corresponding to $f(x)$. However, we seem to have a problem: in order to use rejection sampling most efficiently, it is required to know the maximum value of $f(x)$, which corresponds to the largest permanent of all $n\times n$ submatrices $A$ of a given $m\times m$ unitary matrix $U$. One might expect that computing this would be computationally hard. Without any bound on this quantity, we would be forced to use the trivial bound $\mu=1$, corresponding to $\binom{m}{n}$ permanent computations being required to obtain one sample from $P$.

Fortunately, as we are only attempting to perform approximate boson sampling, we only require a good estimate of $\mu=\max_{x}f(x)$. Recall that it was argued in Aaronson and Arkhipov (2011) that sampling from a distribution within total variation distance $\epsilon$ of the real boson sampling distribution $P_{\mathrm{BS}}$ should be computationally difficult, for some small constant $\epsilon$. Imagine that our guess $\widetilde{\mu}$ for $\mu$ is too small, such that $\sum_{x,f(x)>\widetilde{\mu}}f(x)=\epsilon>0$. Then if $x$ is sampled uniformly at random and $f(x)>\widetilde{\mu}$, step 2 of the rejection sampling algorithm will fail. If we modify the rejection sampling algorithm to simply produce a new uniform sample in this case and repeat, it is easy to see that we can view the modified algorithm as sampling from the truncated distribution $P_{\operatorname{low}}$ with probability mass function

[TABLE]

Then the total variation distance between $P_{\operatorname{low}}$ and $P$ is

[TABLE]

So if the probability mass of $P$ above $\mu$ is at most $\epsilon$, we have sampled from a distribution within distance $\epsilon$ of $P$. We have found that we are able to use a simple random restart hill climbing algorithm to provide a suitable estimate of $\mu$ with $\mathcal{O}(m^{2}n)$ computations of $n\times n$ matrix permanents.

Our random restart hill climbing algorithm works as follows. We start by randomly sampling a submatrix, represented by the tuple $S$, from the uniform distribution. For one pass, we greedily try replacing each row in the sampled submatrix by each row from $A$ in turn, accepting only if this increases $\Pr(S)$ while also making sure to avoid selecting the same row twice. We perform repeat passes until there is no improvement of the probability over a complete pass. At this point we randomly resample a new starting submatrix and repeat from the beginning. The total number of permanent calculations for one pass is $n(m-n)$.

This method is not guaranteed to find a global maximum. However, in our experiments for $n\leq 7$ where we can still compute the full probability mass function exactly, we found the estimates for the maximum probability to be exactly equal to the global maximum in the overwhelming majority of cases. In the range $8\leq n\leq 12$ where we no longer are able to compute the exact maximum probability, the bounds from our hill-climbing algorithm also allowed us to sample using rejection sampling efficiently and then compare our results with our Metropolised independence sampler (qv). This provided further evidence for both sampling techniques.

III.3 Metropolised independence sampling

Markov Chain Monte Carlo (MCMC) methods have been long been the standard tool within statistics for the task of sampling from complicated probability distributions. We employ for our problem of boson sampling a special case of the MCMC procedure known as Metropolised independence sampling Liu (1996, 2004) (MIS), which enables us to sample efficiently from the target distribution.

The general principle of MCMC methods is to construct a Markov chain which has the target distribution (in our case, the boson sampling distribution of tuples) as its stationary distribution. The overall MCMC method requires us only to specify a proposal probability distribution, in our case over tuples $(s_{1},\dots,s_{m})$. At each turn we sample from this probability distribution and accept the new state according to the prespecified rules of the MCMC method. In the general case the proposal distribution can depend on the current state of the chain but in MIS these proposals are entirely independent. By choosing a proposal distribution which is not too far from our target boson sampling distribution, this simplification turns out to be particularly useful for us. Empirically we show that it not only ensures fast convergence but a simple and efficient thinning procedure almost completely eliminates dependence between successive samples.

Let $Y$ be some proposal distribution over the tuples in $\Phi_{m,n}$ with probability mass function $g(x)$. As before, let $P$ be an instance of the CFS restricted (subnormalised) boson sampling distribution over these tuples, with probability mass function $f(x)$. Starting at some random (according to $Y$) tuple $x$, propose a new random tuple $x^{\prime}$. We accept the proposal and transition from $x$ to $x^{\prime}$ with probability

[TABLE]

where we note that any normalisation factor multiplying probabilities would cancel at this stage. Repeating this generates a Markov chain with stationary distribution equal to the target distribution, $P$. If one can argue that the chain has converged then each sample from the generated chain after this point will be from $P$. Each state in the Markov chain requires generating a single tuple sample from $Y$, one evaluation of $h(x^{\prime})$ and one evaluation of $g(x^{\prime})$, which entails one computation of an $n\times n$ complex valued matrix permanent.

It should be noted that it takes some time for the Markov chain to converge to its stationary distribution, and so some number of tuples at the beginning of the chain must be discarded (referred to as “burn in”) in order for samples to come from the target distribution. In addition to this, for example as there is a probability that the state of the Markov chain stays the same between iterations, there exists some autocorrelation in the chain. In order to generate independent samples from the Markov chain, one usually throws away a number of states (referred to as “thinning”) at least equal to the autocorrelation period of the chain before harvesting a sample state. Both burn period $\tau_{\mathrm{burn}}$ and thinning interval $\tau_{\mathrm{thin}}$ depend both on the target distribution and the choice of proposal distribution $Y$. In general, the greater the overlap between target and proposal distributions (as measured by the acceptance probability (10)), the smaller $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$ need be.

A candidate for the proposal distribution is, as above in the case of rejection sampling, the uniform distribution over the tuples. However, we have found that a more suitable proposal distribution is the distribution $P_{\mathrm{D}}$ of distinguishable particles for a given input $U$, as it has a greater overlap with the target distribution. Indeed, our experiments have found that for $n=20$, the acceptance rate when using $P_{\mathrm{D}}$ as the proposal distribution is roughly 40%.

This distribution is tantalisingly similar to the usual boson sampling distribution (7):

[TABLE]

where for a complex matrix $A$ with elements $A_{ij}$, $|A|^{2}$ denotes the matrix with elements $|A_{ij}|^{2}$. However, there is a classical algorithm which can sample from this distribution in time $O(mn)$ Aaronson and Arkhipov (2013).

Using $P_{\mathrm{D}}$ as the proposal distribution, one evaluation of $f(x^{\prime})$ requires computing one $n\times n$ real valued matrix permanent. So, in order to generate the first tuple in a sample, $\tau_{\mathrm{burn}}$ $n\times n$ real valued matrix permanents and $\tau_{\mathrm{burn}}$ $n\times n$ complex valued matrix permanents must be computed. Each subsequent sample requires $\tau_{\mathrm{thin}}$ $n\times n$ real valued and complex valued matrix permanents to be computed.

Testing (see Section IV.4) has suggested that, using the distinguishable particle proposal distribution, it is sufficient for $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$ to grow very slowly with $n$, and for neither to grow above 100 for $n\leq 25$.

Adapting our MIS-based method to carry out the scattershot boson sampling problem is simple. For each sample that we wish to output, we can first sample (efficiently) from the uniform distribution on $n$-fold input modes (which fixes the columns of $U$ contributing to submatrices), before running the algorithm in the way described above for a single sample. In this case, $\tau_{\mathrm{thin}}$ becomes meaningless and we are only interested in $\tau_{\mathrm{burn}}$, as we start a new chain for each sample.

Alternatively, our proposal distribution can be changed to include the uniform distribution over $n$-fold input modes, meaning that each state in a given Markov chain can correspond to a different input configuration.

Although in this work we do not examine the more general situation of boson sampling where there can be more than one boson in an output mode, we anticipate that relaxing the CFS restriction will not increase the run time of the MIS method. In fact, it is possible that the average run time could be decreased with this relaxation, as there exists an algorithm for computing the permanent which is exponential in matrix rank, rather than matrix size Barvinok (1996). However, due to there usually existing very large permanents of sub-unitary matrices with many repeated rows, relaxing the CFS restraint has an adverse effect on the average run time of our rejection sampling method.

IV Numerical results

For concreteness and convenience, we will examine the performance of our samplers for the case where the number of modes scales exactly as the square of the number of photons, i.e. $m=n^{2}$. This regime is generally expected to be hard to simulate classically Aaronson and Arkhipov (2011); Latmiral et al. (2016).

IV.1 Computing matrix permanents

In order to make scaling based arguments for the expected run time of our sampling methods, we need to both confirm that we are able to compute permanents of $n\times n$ matrices in $\mathcal{O}(n2^{n})$ time and compute the associated constant of proportionality.

Timing results are presented in fig. 3, where we see the expected scaling for large $n$. We attribute the difference between the fitting function $t(n)=cn2^{n}$ and the data points for small $n$ to extra operations in our permanent computation code (checks, memory allocation etc.), on top of an implementation of Ryser’s algorithm. From the fit we were able to calculate $c$ for the case of computing permanents of complex and real valued matrices on a personal computer111Specifically, numerical experiments in this paper were performed on a Dell Latitude E5450 laptop with 2.30GHz Intel Core i5-5300U CPU finding that $c=5.180\times 10^{-10}$ and $c=2.106\times 10^{-10}$ respectively. In the case of computing permanents of real valued matrices we note that all values are in fact nonnegative, and therefore it may be beneficial to use the fully-polynomial randomized approximation scheme of Jerrum, Sinclair and Vigoda Jerrum et al. (2004).

Importantly, these values (and approximate values for more powerful machines) allow us to predict the average run time for our sampling methods for any value of $n$.

IV.2 Rejection sampling and random restart hill climbing

In addition to how long it takes to compute a single matrix permanent, we are also interested in how many matrix permanent evaluations are required before a sample value is output.

In the case of rejection sampling with a uniform proposal distribution, we must compute the upper bound $\mu$ using our random restart hill climbing algorithm before implementing the rejection sampling algorithm. We emphasise here that $\mu$ need only be computed once per input $U$, and so for large enough sample sizes, the computational contribution of computing $\mu$ gets washed out. For concreteness, we focus here on outputting a sample consisting of 100 tuples, and present the number of required permanent computations averaged over these 100 tuples.

In fig. 4 we see that, by using rejection sampling instead of brute force sampling, one can reduce the number of computations required at a rate which increases rapidly with $n$. In fact, a significant portion of the computation for all values of $n$ investigated here is spent getting an estimate of $\mu$. It should be noted, however, that this is an artefact of our choice of number of restarts implemented in our random restart hill climbing algorithm. We chose the number of restarts to be $4m$, based on empirical evidence that this results in the probability mass above this estimate being small. For $600$ random instances of the problem with $2\leq n\leq 7$, we found that the largest submatrix permanent was found $598$ times. When the largest value was not found, the probability mass above the estimate was found to be small ($<2\times 10^{-5}$).

As an example for reference, on a personal computer the average time to produce one sample for $n=12$ (averaged over 100 samples) was 2.03s.

IV.3 Metropolised independence sampling

Determining the theoretical runtime of MIS is straightforward: $\tau_{\mathrm{thin}}$ and $\tau_{\mathrm{burn}}$, together with the time to compute a matrix permanent, completely determine the time to produce a sample, as discussed above. We include these bounds in fig. 4, taking $\tau_{\mathrm{thin}}=\tau_{\mathrm{burn}}=100$, and amortised over 100 samples. Observe that for $n=20$, each tuple can be output in 1.58s on a personal computer; even for $n=25$, only 62.5s are required.

IV.4 Assessment of our sampling methods

How can we confirm that our samplers are sampling from an approximately correct distribution? Boson sampling is believed to be computationally hard to verify Aaronson and Arkhipov (2011); Gogolin et al. (2013); Aaronson and Arkhipov (2013), so, to date, experimental implementations of the problem have relied upon circumstantial evidence to support the claim that they are, in fact, sampling from the target distribution. This, coupled with the fact that reliability of MCMC methods is notoriously difficult to certify, requires us to proceed with caution.

However, the fact that we have 3 completely different samplers, each of which we can ask to provide us with a sample for the same instance of the problem, allows us to gather mutually supportive evidence for them all sampling from the correct target distribution; if we cannot reject that all pairs of samples come from the same underlying distribution, we can have some confidence that they are all in fact sampling from the target distribution.

A natural quantity to use to compare the distributions on matrices that we consider is $|\mathrm{Per}(A_{S})|^{2}$ of the sampled matrix $A_{S}$, which is equal to the probability of the sampled matrix in the boson sampling distribution; for convenience we actually use $\log|\mathrm{Per}(A_{S})|^{2}$.

We generated samples of size 20000 of tuples (assigning each of them an index of an integer value between 0 and $\binom{m}{n}-1$), together with their corresponding $-\log(|\mathrm{Per}(A_{S})|^{2})$, using various samplers and for problems of size $n=3,7,12,20$. The data are plotted in fig. 5.

We expect that the distribution of $-\log(|\mathrm{Per}(A_{S})|^{2})$ differs between samples of indistinguishable bosons and distinguishable particles. This feature manifests itself clearly for each attempt and each method of attempted boson sampling; the distinguishable particle values (yellow) are distributed visibly differently to those of other sampling methods (all other colours). More rigorously, we performed 2-sample bootstrap Kolmogorov-Smirnov (KS) tests Praestgaard (1995); Sekhon (2011) comparing the distribution of $-\log(|\mathrm{Per}(A_{S})|^{2})$ values for tuples obtained by sampling from the distinguishable particle distribution with those obtained using our classical boson sampling methods. The result was that we could reject the null hypothesis that this feature of the samples was distributed identically at the 0.1% significance level (i.e. $p<0.001$ for all samplers and all problem sizes considered). This result was found to be consistent over 5 repetitions of taking a sample for each input unitary.

Not only do we see that our sampling methods give significantly differently distributed $-\log(|\mathrm{Per}(A_{S})|^{2})$ to the distinguishable particle sampler, we see no significant difference between these distributions for each of our sampling methods. Applying the same 2-sample bootstrap KS tests to compare these distributions returned that we cannot reject the null hypothesis that they are distributed identically at any reasonable significance level. We see this as evidence that all of our samplers are sampling from the target boson sampling distribution.

Due to some sampling methods being more efficient than others, we witness problem sizes beyond which certain methods become unreasonably slow. As a rough guide, on a personal computer these problem sizes are $n=7$ for brute force sampling and $n=12$ for rejection sampling. Thus, in the case of $n=20$, we compare one MIS sampler (with $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=100$) with another MIS sampler (with $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=1000$).

In table 1 we present the $p$-values obtained when performing 2-sample bootstrap Kolmogorov-Smirnov tests on the values $-\log(|\mathrm{Per}(A_{S})|^{2})$ for sets of $20000$ tuples obtained with various samplers. Independent samples were generated $5$ times in order to repeat the tests. The null hypothesis in this test is that the values are identically distributed. As $p$-values are uniformly distributed when the null hypothesis is true, we see no significance in any relatively small values in table 1; they are to be expected when quoting a large number of $p$-values.

Note that the assessment of MCMC algorithms is non-trivial: in particular, slow convergence to the target distribution and autocorrelation within the chain can result in an erroneous sample being output. We therefore also performed a modified version of the likelihood ratio test described by Bentivegna et al. Bentivegna et al. (2014). Following the notation of Bentivegna et al. (2014), the test works as follows. We define two hypotheses: $\mathcal{Q}$, the indistinguishable boson hypothesis, and $\mathcal{R}$, an alternative hypothesis. In our case, the alternative hypothesis $\mathcal{R}$ is the distinguishable particle hypothesis. Let $q_{x}$ be the probability of seeing the sampled event $x$ according to hypothesis $\mathcal{Q}$, and $r_{x}$ be the corresponding probability under hypothesis $\mathcal{R}$. Assigning equal priors to each hypothesis, Bayes’ theorem tells us that

[TABLE]

where, for example, $P\left(\mathcal{Q}|N_{\mathrm{events}}\right)$ is the probability of hypothesis $\mathcal{Q}$ being correct given that a sample of size $N_{\mathrm{events}}$ with events $\left\{k_{x}\right\}$ was obtained. A large value of $\mathcal{X}$ would correspond with us having a high degree of confidence that hypothesis $\mathcal{Q}$ provides a better description of the observed data than hypothesis $\mathcal{R}$.

It is convenient to rewrite equation (12) in terms of the probability assigned to hypothesis $\mathcal{Q}$,

[TABLE]

which is now normalised such that $P\left(\mathcal{Q}|N_{\mathrm{events}}\right)+P\left(\mathcal{R}|N_{\mathrm{events}}\right)=1$.

One, perhaps subtle, point here is that $q_{x}$ and $r_{x}$ do not simply correspond to $\left|\mathrm{Per}\left(A_{x}\right)\right|^{2}$ and $\mathrm{Per}\left(\left|A_{x}\right|^{2}\right)$ respectively. The reason for this is that we are restricted to the CFS, and the probability of a sample being collision-free differs between indistinguishable bosons and distinguishable particles. So $q_{x}$ and $r_{x}$ must be normalised independently, such that they independently sum to 1 over all events $x$. Unfortunately, doing this exactly would require summing all probabilities in the CFS for the specific instance of the problem being considered – a computationally daunting task, as previously discussed. To approximate the normalisation of $q_{x}$ efficiently, here we instead average the probability of the output being collision-free over the Haar measure, as reported in Aaronson and Arkhipov (2011); Arkhipov and Kuperberg (2012):

[TABLE]

For the hypothesis $\mathcal{R}$, we can efficiently sample output tuples, and so the ratio of collision-free tuples to tuples with collisions in a large sample provides an approximate normalisation for $r_{x}$.

For each problem size, we use this likelihood ratio test to assess the performance of MIS samplers with different $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$. The reasoning is that, as the proposal distribution for the sampler is the distinguishable particle distribution, we might expect that if the chain has not converged to the target distribution, this will manifest itself as samples looking more like they are from the distinguishable particle distribution than they should. Also, we expect that samplers with larger $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$ are more likely to sample from the target distribution. Because of this, we expect to be able to observe a point at which increasing $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$ has, on average, no effect on the outcome of a likelihood ratio test between hypotheses $\mathcal{Q}$ and $\mathcal{R}$.

The data corresponding to this for $n=3,7,12,20$ are presented in fig. 6. We can see that, indeed, for a given problem size, increasing $\tau_{\mathrm{burn}}$ and $\tau_{\mathrm{thin}}$ beyond a point does not improve the likelihood ratio test results. Not only this, but we observe that this point gradually increases in size with the problem size $n$. For example, with problem size $n=3$, $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=5$ provides results comparable with $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=1000$. However, with problem size $n=20$, $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=5$ provides markedly worse samples on average than $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=1000$. In fact, $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=100$ seems to be the point at which increasing these values provides no benefit.

Additionally, in fig. 7 we see that autocorrelation between $-\log(|\mathrm{Per}(A_{S})|^{2})$ for 100 lags in an un-thinned chain is negligible for random instances of the problem for $n=3,7,12,20,25$. This is further supported by the autocorrelation being negligible at all lags for chains generated by applying our thinning procedure to the un-thinned chains.

Based on this, we claim that $\tau_{\mathrm{burn}}=\tau_{\mathrm{thin}}=100$ is sufficient for MIS for problems up to size $n=20$, and is likely to be sufficient for values of $n$ greater than 20. The consequence of this is that a sample from the boson sampling distribution can be output, using a classical computer, in roughly the time that it takes to compute 100 $n\times n$ real valued matrix permanents and 100 $n\times n$ complex valued matrix permanents. For $n=20$ and $m=400$, this equates to a computational saving of a factor of $\sim 10^{31}$ when compared with brute force sampling.

IV.5 Moving to larger $n$

Although MIS is efficient in the number of matrix permanent computations required, it still requires the inefficient computation of matrix permanents by Ryser’s algorithm Ryser (1963). As such, gathering a $20000$ tuple sample as was obtained in fig. 5 is challenging for $n>20$ using readily available hardware. We were able to perform boson sampling using MIS for $n=25$ and $m=625$, with a sample size of $1000$ tuples for plotting the distributions of indexed tuples and $-\log(|\mathrm{Per}(A_{S})|^{2})$ of these sampled tuples, alongside a sample from the distinguishable particle distribution. We also performed a likelihood-ratio test to compute the probability of the indistinguishable boson hypothesis against the distinguishable particle hypothesis, averaged over 100 samples of 25 tuples for MIS with different values of $\tau$.

In addition to this, we performed boson sampling using MIS for $n=30$ and $m=900$, obtaining a $250$ tuple sample with $\tau=100$, and compared the convergence of $P_{\mathrm{ind}}$ with other values of $\tau$.

All of the data for $n>20$ was taken on a cluster of 4 servers222 2 dual-socket Intel Xeon CPU E5-2697A v4 @ 2.60GHz and 2 dual-socket Intel Xeon CPU E5-2680 v3 @ 2.50GHz at the University of Bristol, allowing for $122$ MIS instances (and therefore $122$ Markov chains) run in parallel (vs a maximum of $4$ on our personal computer). The data are presented in fig. 8. For $n=25,30$, we notice no change in the performance of MIS beyond having to compute permanents of larger matrices.

V Comparison with experimental implementations

Limiting ourselves to the standard boson sampling problem for now, we will consider how the time to obtain a single tuple sample scales with experimental parameters for a physical boson sampling device.

Parameters which have an effect on the rate at which samples can be generated experimentally are: the $n$-photon generation rate $R(n)$, the single photon transmission probability $\eta$, and the Haar-averaged probability of the $n$-photon state being in the CFS given by eqn. (14), which we take to be an approximation of the relative size of the CFS for any boson sampling input $U$.

It is important to notice here that in fact $\eta$ depends on $n$. This is because the loss is composed of fixed loss $1-\eta_{\mathrm{f}}$ (source efficiency, detection loss, insertion loss etc.) and loss within the circuit, $1-\eta_{\mathrm{c}}$. The circuit depth – and therefore the amount of material (or number of components) each photon must propagate through – necessarily grows with the size of the problem, as non-trivial interference between all photons and across all modes must occur.

Circuits generally considered for boson sampling have depth $d$ which scales at least linearly with $m$ Reck et al. (1994); Clements et al. (2016) and hence at least quadratically with $n$. However, Aaronson and Arkhipov have shown that boson sampling can be implemented using a different circuit construction with depth $d=\mathcal{O}(n\log m)$ Aaronson and Arkhipov (2011) if the input modes are fixed, and recently a circuit and predictions of performance with a linear mode scaling of $m=4n$ have been reported Wang et al. (2016a).

This might suggest that $d=\mathcal{O}(n)$ could also be sufficient, although note that when $m=\mathcal{O}(n)$ the probability of a sample being collision-free is exponentially small in $n$ (see eqn. (14)), and no theoretical argument for computational hardness is known.

Here we compare the rates with which we expect to obtain samples from a boson sampling experiment for the two regimes $d=n^{2}$, $d=4n$ where experimental attempts have been documented Spring et al. (2013); Crespi et al. (2013); Tillmann et al. (2013); Broome et al. (2013); Carolan et al. (2014); Spagnolo et al. (2014); Bentivegna et al. (2015); Carolan et al. (2015); He et al. (2016); Wang et al. (2016a); Loredo et al. (2017). The transmission probability of each photon through the circuit is $\eta_{\mathrm{c}}={\eta_{0}}^{d}$, where $\eta_{0}$ is the probability of a photon surviving per 2-mode coupling length (or component) in the interferometer. Using these parameters, we can write the time in which we expect to obtain a single sample from a quantum device as a function of the problem size:

[TABLE]

where $\eta=\eta_{\mathrm{f}}\eta_{0}^{d}$ Expanding $\eta$ for linear and quadratic mode scaling respectively, we see that

[TABLE]

where it should be emphasised that $\eta_{\mathrm{f}},\eta_{0}<1$.

Using these bounds, we see that the time required to obtain a sample scales exponentially in either $n^{2}$ or $n^{3}$. However, classical computation of the permanent using Ryser’s algorithm scales only exponentially in $n$. Therefore, for very large $n$ and realistic loss parameters, the runtime scaling of a classical sampling algorithm would be favourable compared with the boson sampling experiment, even if the algorithm needed to compute exponentially many permanents to obtain one sample.

In order to provide a meaningful comparison of our classical boson sampling methods to an experimental implementation, we focus on the current leading experimental approach by Wang et al. Wang et al. (2016a). This approach involves the use of a quantum dot based photon source, and has resulted in increased photon numbers and significantly higher rates compared to the current best SPDC based demonstrations (see Extended Data Table 1 in Wang et al. (2016a)). In fig. 9 (also see fig.2d in the main text), we consider the rate at which a similar experiment with projected near-future parameters Wang et al. (2016a) can generate sample values with both $m=4n$ and $m=n^{2}$, and compare this to the average sample value generation rate for a personal computer running the MIS algorithm for boson sampling.

Specifically, in Wang et al. (2016a) Wang et al. report a quantum dot based photon source with $R(n)=76n^{-1}$MHz. Also reported are loss parameter values predicted to be achievable in the near future:

[TABLE]

where $\eta_{\mathrm{QD}}$ is the single photon source end-user brightness ($=0.74$), $\eta_{\mathrm{de}}$ is the demultiplexing efficiency for each channel ($=0.845$) and $\eta_{\mathrm{det}}$ is the efficiency of each detector ($=0.95$). The authors of Wang et al. (2016a) report the efficiency of their 9 mode interferometer to be $\sim 99\%$, and so we approximate the efficiency of a similar interferometer scaled to $m$ modes as $\eta_{\mathrm{c}}=0.99^{\frac{m}{9}}$ and hence $\eta_{0}=0.99^{1/9}$ as we are assuming $d=m$.

We see that, although the personal computer is outperformed for a small problem size, beyond $n\approx 14$ the personal computer outperforms the experimental boson sampler. Clearly, the point at which the classical device becomes dominant would increase if the experiment had access to an $n$-photon state generator with greater repetition rate. On the other hand, if the classical device were more powerful (i.e. could perform more floating point operations per second), it would start to become dominant at even smaller problem sizes. More important is that the rate falls off noticeably more rapidly for the experimental sampler with these loss parameters, for both regimes of scaling the number of modes with the number of photons.

To emphasise this point, we perform the following analysis. Assuming that our MIS sampler continues to perform equally well for larger instance sizes, we can compare its runtime with current and future experiments. The classical runtime for an instance of size $n$ bosons in $m=n^{2}$ modes can be estimated as

[TABLE]

where $3\times 10^{-15}s$ is the time scaling for computing one real and one complex permanent recently reported for the supercomputer Tianhe 2 Wu et al. (2016). We note here that we do not assume that we can use the efficient (i.e Gray code ordered) version of Ryser’s algorithm here, as we would need to parallelise the computation of the permanent to deal with large $n$ problems. We define the quantum advantage (QA) as the order of magnitude improvement in quantum runtime versus classical runtime,

[TABLE]

We now consider two possible definitions of quantum supremacy. First, we can define supremacy as a speedup so large that it is unlikely to be overcome by algorithmic or hardware improvements to the classical sampler, for which we choose a speedup of ten orders of magnitude. Secondly, we may wish to define supremacy as the point at which a computational task is performed in a practical runtime on a quantum device, for which we choose under a week, but in an impractical runtime on a classical device, for which we choose over a century.

These can be summarised as

[TABLE]

For a quadratic mode scaling and a photon source with rate $R(n)=10$GHz (beyond the capabilities of the current best single photon sources), we plot QA against $n$ and $\eta$ in fig.1 of the main text, where $q_{t}^{\mathrm{quad}}$ is used. In order to achieve both quantum supremacy criteria in this regime, we see that it is required to have $\eta>0.6$ and $n>60$. We can perform the same analysis with $q_{t}^{\mathrm{lin}}$ and $R(n)=76n^{-1}$MHz as per Wang et al. (2016a). The data for this are plotted in fig. 9. Although using a linear mode scaling increases $\eta$, we see that the requirements for quantum supremacy become even more stringent. Indeed, in order to achieve both quantum supremacy criteria here we require $\eta>0.9$ and $n>70$.

VI Lossy boson sampling

We have seen that loss can cause significant difficulties in the setting of the standard boson sampling problem, leading to a classical sampler outperforming the quantum experiment. In this section, we explore whether modifying the problem to build loss in from the start could mitigate these issues. One might expect that losing photons would also make the boson sampling problem easier classically, because the matrices whose permanents are required to be computed would be smaller. However, it is conceivable that the quantum advantage would increase, enabling a demonstration of quantum supremacy.

A lossy variant of boson sampling was introduced by Aaronson and Brod Aaronson and Brod (2016). In this variant, we assume that $n-k$ photons are lost before entering the circuit enacting the linear optical transfer matrix, so $k$ photons remain. Probabilities in this setting are not determined directly by $\left|\mathrm{Per}\left(A_{S}\right)\right|^{2}$ for some $m$-tuple $S$, but by the average of this quantity over all possible ways of losing $n-k$ photons from $n$ photons:

[TABLE]

where $\Lambda_{n,k}$ is the set of $k$-subsets of $\{1,\dots,n\}$ and $A_{S,T}$ is the $k\times k$ submatrix of $A$ obtained by taking columns of $A$ according to $T$ and rows of $A$ according to $S$, which remains a subset of $\{1,\dots,m\}$. Note that once again we restrict to the collision-free subspace, making the assumption that the probability of a collision is low enough that this does not significantly affect the probabilities.

It was shown by Aaronson and Brod Aaronson and Brod (2016) that, if at most a constant number of photons in total are lost before entering the circuit, the lossy boson sampling problem remains hard, under the same assumptions as the original boson sampling problem. However, it was left as an open problem whether the same holds true in the more realistic setting where a constant fraction of photons are lost. Another open problem in Aaronson and Brod (2016) was to generalise the loss model to include loss within and after the linear optical circuit. In Berry and Lvovsky (2010) it is shown that uniform (i.e. mode-independent) loss channels can commute with linear optics. Here we prove the slightly stronger result, that as long as the overall transfer matrix is proportional to a unitary, loss can always be considered at the input even if the physical loss channels, wherever they occur, are not uniform (i.e are mode-dependent).

Consider a boson sampling device consisting of an ideal unitary linear optical transformation $U$ on a set of $m$ optical modes which is preceded or succeeded by path-independent loss. This loss can be modelled by considering a set of $m$ additional virtual ancilla modes such that the optical transfer matrix on all $2m$ modes remains unitary. A uniform transmission probability of $\eta$ can then be described by beamsplitters coupling each mode to its corresponding ancilla, resulting in the transfer matrix

[TABLE]

and so including the interferometer, the full transfer matrices for input and output losses are

[TABLE]

Any m-mode optical state can be expressed in a coherent state basis Gerry and Knight (2005):

[TABLE]

where $|\bm{\alpha}\rangle=|(\alpha_{1},\alpha_{2},...,\alpha_{m})^{T}\rangle\equiv\bigotimes_{i=1}^{m}|\alpha_{i}\rangle$ is an $m$-mode coherent state and $\langle\bm{\beta}|=\langle(\beta_{1},\beta_{2},...\beta_{m})^{T}|=\bigotimes_{i=1}^{m}\langle\beta_{i}|$. A coherent state evolves under a transfer matrix $T$ as

[TABLE]

It can then be shown that when the initial state contains vacuum in all ancilla modes, $\rho=\rho_{m}\otimes|\bm{0}\rangle\!\langle\bm{0}|_{\text{an}}$, the same state is produced in the $m$ system modes under the transformations $M_{O}$ and $M_{I}$

[TABLE]

More generally, wherever loss occurs in the experiment, the overall transfer matrix $K$ on the system modes can be efficiently characterised Laing and O’Brien (2012); Rahimi-Keshari et al. (2013). Since path-dependent loss is usually small in experiments, and can be mitigated by interferometer design Clements et al. (2016), the matrix $K/||K||_{2}\approx U$. The matrix $K$ can then be embedded into a larger unitary matrix acting on additional modes as before. We note all unitary dilations of $K$, $M_{K}\in U(m+p)$ where $p\geq m$, can be parameterised using the Cosine-Sine decomposition as

[TABLE]

where $K=A\cos{(\Theta)}B^{{\dagger}}$, with $A,B\in U(m)$ and $\cos{(\Theta)}=\text{diag}(\cos{\theta_{1}},...,\cos{\theta_{m}})$ with $\theta_{1}\leq\theta_{2}\leq..\leq\theta_{m}$, is a singular value decomposition of $K$ and $X,Y\in U(p)$. In fact, all unitary dilations are related by the choice of $X$ and $Y$ Horn and Johnson (1991); Allen and Arceo (2006). Since $\mathcal{U}(Y)\left|\mathbf{0}\right\rangle=\left|\mathbf{0}\right\rangle$ and the choice of $X$ does not affect $\rho_{K}$ using the cyclic property of the trace as above, setting $\eta=||K||_{2}^{2}$, we see that $\rho_{K}=\rho_{I}$. Moreover, we have shown that all unitary dilations of a transfer matrix produce the same output state and therefore any boson sampling experiment with overall path-independent losses is equivalent to introducing uniform loss channels with transmission probability $\eta$ at the input, followed by the ideal unitary evolution.

Our MIS method can readily be adapted to deal with loss at the input, by inserting an initial step for each tuple to be output, which generates a uniformly random input subset $T$. This would be followed by the usual MIS method with permanents of $k\times k$ submatrices computed. We described a similar adaptation in Sec. IV.3 to deal with scattershot boson sampling. The core part of the classical sampling procedure for both the lossy and scattershot variants therefore follows precisely that of standard boson sampling. From our analysis of the required burn in period for MIS, we can see that the performance of our sampler will be similar to the standard boson sampling case. That is, it is likely that lossy and scattershot boson sampling is no more difficult classically than standard boson sampling.

The modified MIS sampler can then be compared with experimental performance. Using the fixed loss regime of Aaronson and Brod (2016), we adjusted the quantum and classical runtimes and optimised QA for the case where up to 2 photons were allowed to be lost. This resulted in the dashed lines in fig.2d) and e) of the main text, and in fig. 9. We see that allowing this loss relaxes the requirements for quantum supremacy, but does not significantly change the experimental regimes required.

In current and future experiments, where $\eta\ll 1$, more photons than this are usually lost. Again using expected loss parameter values from Wang et al. (2016a), we find that the optimal number of photons to lose (in terms of performance enhancement relative to the classical sampler) increases with $n$. The estimated experimental time to produce a tuple when $n=20$, $k=12$ and $m=80$ is $11.8\mu$s, where this time encompasses the expected number of repetitions required to detect exactly $k$ photons. With the same parameter values and problem size, the estimated MIS runtime on a personal computer is $3.9$ms. This amounts to a factor of $333$ performance enhancement for the experiment over the personal computer.

However, the MIS sampler is not the only way in which one could attempt to sample from the lossy boson sampling distribution: indeed, it could even be the case that enough loss could render the boson sampling distribution easy to sample from classically. For example, it was shown in Bremner et al. (2016b) that a class of quantum circuits whose output probability distributions are likely to be hard to sample from classically becomes easy in the presence of noise. We also remark that Rahimi-Keshari, Ralph and Caves Rahimi-Keshari et al. (2016) have proven classical simulability of boson sampling under various physically motivated models of errors (e.g. loss, mode mismatch, and dark counts), and that Aaronson and Brod have speculated Aaronson and Brod (2016) that each probability in the small-$k$ lossy boson sampling distribution tends to the product of the squared 2-norms of the rows of $A_{S}$ – an efficiently computable quantity.

Lossy boson sampling also suffers from an additional difficulty when one considers performing likelihood based verification techniques. In particular, computing the likelihood of a sample according to the indistinguishable boson hypothesis requires computing eqn. (19) for each tuple $S$ in a sample. The sum over all $T\in\Lambda_{n,k}$ necessitates a factor of $\binom{n}{k}$ slowdown in computing the likelihood relative to a standard (i.e. lossless) $k$-photon boson sampling experiment. The worst case for this is when $k=n/2$ photons survive, amounting to $\sim 2^{n}/\sqrt{\pi n/2}$ matrix permanent computations per likelihood calculation. To avoid this exponential slowdown, one would be forced to use a nonstandard verification technique Aaronson and Arkhipov (2013); Walschaers et al. (2016), or devise a more efficient means of computing (19) than naïve evaluation of the sum.

VII Discussion

We finish by discussing the limitations of our techniques, and prospects for future experimental and classical improvements.

VII.1 Limitations

Our sampling methods have some limitations that we now set out. Our rejection sampling method which is applicable for $n$ up to approximately $12$, is guaranteed to give independent samples. Our approach to improving the efficiency of the method requires us first to estimate the maximum probability $\mu$ of the distribution. If this estimate is very poor, that is if there is a large proportion of the probability mass greater than $\mu$, then there is a risk we will not sample from the boson sampling distribution correctly. From our computational experiments, we can see that this appears to be very unlikely to occur at least in the range we are able to test.

For our MIS method, we have to be concerned both with how long to wait until convergence and the related question of potential dependence between samples. It is unfortunately not possible to prove that the samples we take are entirely independent of each other or indeed of the initial state. Although this question of dependence in boson sampling is not one that has arisen previously, it is desirable for samples to be taken independently.

We have mitigated this problem in two ways. The first is by our choice of MIS as a sampling method, which means that samples from our proposal distribution are taken independently. The second is by our use of a thinning procedure. In fig. 7 we see that this thinning procedure successfully reduces the auto-correlation of the computed probabilities of the sampled submatrices to near zero. In addition to this, in fig. 6 we see that extending the thinning procedure beyond a certain point causes no improvement in the results of a likelihood ratio test as outlined in Bentivegna et al. (2014). However, it is still possible that there remains dependence on some other detailed property of the sampled submatrices which is not exposed by looking at the computed probabilities or likelihood ratios.

VII.2 Prospects for experimental boson sampling

The fact that producing boson sampling machines requires low loss rates has been known to be a hurdle since the conception of the problem Aaronson and Arkhipov (2011). However, the extent of this hurdle has not been well understood, not least due to the lack of an optimised classical boson sampling algorithm.

Our results suggest that with loss parameters and photon numbers expected in the near future, no straightforward or significant quantum supremacy is achievable via boson sampling. This is not to say that it is impossible to build a boson sampling machine which demonstrates quantum supremacy. However, it is absolutely necessary to further decrease loss rates, or develop efficient fault-tolerance techniques targeted at linear optics in order to achieve this. A machine performing boson sampling with tens of photons in thousands of modes is likely required if we are to witness an experiment significantly outperforming a classical simulator. We note the engineering feat required to build an interferometer of this size alone may prohibit its realisation in the near future.

One approach to mitigating loss in an experiment would be to reduce circuit depth. In fig. 1e) of the main text and fig. 9 we have assumed a circuit depth of $d=m$ to infer $\eta$ for the specific experimental points plotted (solid and open white points). Using more exotic structures such as the 3D platform used in Crespi et al. (2016), it may be possible to reduce this scaling, and hence increase $\eta$.

In addition to the problems caused by photon loss, it is also necessary to build reconfigurable linear optical circuits capable of performing high precision transformations on modes Leverrier and García-Patrón (2015); Kalai and Kindler (2014); Arkhipov (2015b), to generate highly indistinguishable photons Tichy (2015); Shchesnovich (2014) and to limit detector noise and higher order photon number terms in the input state Rahimi-Keshari et al. (2016) in order to perform a classically intractable task. We note that although there exist methods to aid in achieving these things – such as the spectral filtering of single photons or adding additional components to perform higher fidelity beamsplitter transformations Miller (2015); Wilkes et al. (2016) – they typically result in an increased amount of photon loss.

There have also been proposals for performing boson sampling using other hardware platforms, such as trapped ions Shen et al. (2014) and superconducting qubits Peropadre et al. (2016). Although currently experimentally untested, it is possible these approaches could provide better scalability. However, coherence and gate operation times are likely to limit the rates and sizes of implementations possible in these architectures, e.g. in ref. Peropadre et al. (2016) current estimates suggest up to 20 photons.

VII.3 Prospects for classical boson sampling

Our current approach to simulating a boson sampling experiment assumes a Fock state input of the form given in eqn. (9). Another proposal for boson sampling type experiments is based on Gaussian input states Lund et al. (2014); Hamilton et al. (2016). In this case, Hamilton et al. Hamilton et al. (2016) showed that the probabilities in the output distribution are related to the hafnian Minc (1978) of a submatrix, rather than the permanent. It has been shown by Björklund that hafnians can be computed using an algorithm with runtime matching Ryser’s algorithm for the permanent, up to polynomial factors Björklund (2012). We are hopeful that the techniques we have presented in this paper can be carried over to Gaussian boson sampling, potentially allowing for the classical simulation of experiments designed to generate molecular vibronic spectra Huh et al. (2015).

In terms of performance, it is possible that improvements can be made in a number of ways. Firstly, it is likely that MIS with the distinguishable particle proposal distribution is not optimal for classical boson sampling. Although computing matrix permanents will remain the bottleneck, requiring fewer of these computations per output tuple could result in up to a $\sim 100$ times speed-up. Secondly, the implementation could be scaled up to run on high performance computing hardware. This would require careful parallelisation of both MIS and Ryser’s (or Balasubramanian–Bax/Franklin–Glynn) algorithm. We remark that parallelisation of the most computationally expensive parts of MIS is relatively trivial. As the proposal distribution is independent of the current state, we can distribute the permanent computations across multiple processors. In addition to this, it may be possible to build a more realistic model of the experiment (including, for example, partial distinguishability of photons) to allow for some approximations of matrix permanents to be made, or for permanents of smaller submatrices to be computed with some probability.

VIII Conclusion

We have shown that boson sampling can be simulated classically for a range of parameters that was previously considered likely to be computationally intractable Aaronson and Arkhipov (2011); Preskill (2012); Bentivegna et al. (2015); Latmiral et al. (2016). Our results pose a challenge for both experimentalists and theorists. The challenge for experimentalists is to build low-loss, large-scale linear optical circuits going substantially beyond current technology, while the challenge for theorists is either to develop efficient linear-optical loss-tolerance schemes, or to find some way to distinguish our classical boson sampling simulators from the real boson sampling distribution.

Note that one way of implementing a large-scale loss-free boson sampling experiment is simply to simulate boson sampling on a fault-tolerant universal quantum computer, perhaps based on some other technology than optics. However, in the near term this is unlikely to be a more efficient means of demonstrating quantum supremacy than an approach targeted directly at the underlying hardware Boixo et al. (2016); Bremner et al. (2016a); Lund et al. (2017); Bermejo-Vega et al. (2017).

We close by highlighting another interpretation of our results: as providing an additional tool for classical verification of future boson sampling experiments. Due to the computational difficulty of simulating the boson sampling distribution, existing verification procedures have not been designed to take advantage of statistical information about the true boson sampling distribution. Indeed for $n>7$ this information has simply not been available within a realistic time frame. Our new efficient computational sampling procedure allows us for the first time to test directly whether an experiment is sampling from the correct distribution. Setting the null hypothesis to be that we are sampling from the boson sampling distribution, we can sample directly from the true distribution using our MIS procedure and then perform standard statistical hypothesis tests on the results of any boson sampling experiment. MIS can also be used to certify lossy versions of boson sampling, where previous methods based on computing likelihoods are inefficient, as discussed in Section VI. This will potentially give a new level of confidence that was not previously available in the accuracy of any experimental design and implementation.