Polynomial-Time Classical Simulation of BosonSampling under Any Constant-Fraction Photon Loss

The paper in detail

In BosonSampling, \(N\) photons enter the first \(N\) of \(M\) modes of a Haar-random interferometer \(U\) and are counted at the output. In the fixed-loss model exactly \(N-n\) photons are lost, the surviving set \(S\) being uniformly random, so the probability of output pattern \(\mathbf t\) is \[P_U^{N,n}(\mathbf t)=\frac1{\binom Nn}\sum_{|S|=n}\frac{|\operatorname{Per}(U_{\mathbf t,S})|^2}{\mathbf t!}.\] Permanents are hard to compute, which underlies the claims of quantum advantage. Loss cannot be corrected in this linear-optical setting. Previous algorithms were provably efficient only when the surviving fraction vanished, for example \(n=O(\sqrt{\varepsilon N})\). Many experiments instead keep a fixed fraction of photons.

The key observation is that averaging over which photons survive damps high-order interference. Written as a sum over sets \(J\) of input columns, the order-\(j\) part of the output probability carries the factor \(\binom nj/\binom Nj\le(n/N)^j\). Truncating at degree \(d\) leaves a mean-square error at most \((n/N)^{d+1}/(1-n/N)\) (Proposition 2).

Theorem 2 turns this into a sampler. It runs in \(MnN^3\varepsilon^{-2}\delta^{-1}2^{O(d)}\) arithmetic operations with \(d=O\bigl(1+\log(\varepsilon^{-1}\delta^{-1}(1-n/N)^{-1})/\log(N/n)\bigr)\). For a fixed fraction \(n/N\), \(d=O(\log(1/\varepsilon\delta))\), so \(2^{O(d)}\) is polynomial. The pipeline has four steps. A truncated, signed approximation of the distribution is built. Its prefix marginals are evaluated by random color coding. A signed-marginal sampler then produces a proper distribution over ordered output strings, and sorting gives Fock outcomes. Markov’s inequality converts an average-case error bound into a statement for a \(1-\delta\) fraction of interferometers.

Under the extra condition \(M=\Omega(n^2/(\varepsilon\delta))\), collisions are rare enough that the result covers the full Fock distribution (Corollary 2). For uniform transmission loss, one first samples the number of survivors from a binomial distribution and then runs the fixed-loss algorithm (Corollary 3). The authors list the saturated-mode regime, partial distinguishability and non-ideal detectors as open.

Transferable methods

1. Decompose an average over random survivors

   The output probability is an average over surviving subsets. A Hoeffding (ANOVA) decomposition over input columns shows how strongly each interference order is weighted.

   For an \(n\times N\) matrix \(B\), let \(G_L(B)=(n-|L|)!\sum_{|I|=|L|}|\operatorname{Per}(B_{I,L})|^2\) and define \(H_J=\sum_{L\subseteq J}(-1)^{|J|-|L|}G_L\). Möbius inversion gives \(G_S=\sum_{J\subseteq S}H_J\), because \(\sum_{L\subseteq J\subseteq S}(-1)^{|J|-|L|}=(1-1)^{|S|-|L|}\) vanishes unless \(L=S\). Inserting this into the average over surviving sets \(S\), and counting the \(\binom{N-j}{n-j}\) sets that contain a given \(J\), gives \[\Phi(B)=\frac1{\binom Nn}\sum_{|S|=n}|\operatorname{Per}(B_S)|^2=n!+\sum_{j=1}^{n}\frac{\binom nj}{\binom Nj}\sum_{|J|=j}H_J(B).\] The coefficient \(\binom nj/\binom Nj\) is the probability that \(j\) given photons all survive.

   For a Gaussian matrix \(X\), each \(H_J\) has zero mean in every column it involves. So distinct components are orthogonal, and \(\mathbb E|H_J|^2=(n!)^2/\binom n{|J|}\) by Wick’s theorem. Hence \[\frac{\mathbb E|\Phi-\Phi_{\le d}|^2}{(n!)^2}=\sum_{j>d}\frac{\binom nj}{\binom Nj}\le\sum_{j>d}\Bigl(\frac nN\Bigr)^j\le\frac{(n/N)^{d+1}}{1-n/N}.\] For example, with half the photons surviving, each extra order of truncation halves the bound. Lemma 4 transfers this to Haar submatrices, which are close to Gaussian when \(M\ge C_0n^2N\).

   How to reuse it

   Whenever an observable is an average over random subsets (lost particles, erased qubits, random subsampling), decompose it by subset order. The weight of order \(j\) is then the probability that a given \(j\)-set survives, and truncation is controlled once the components are orthogonal under some reference measure. The cost is that the error bound is an average: here it holds for most interferometers, not all.

2. Sum over subsets with color coding

   Sampling needs prefix marginals of the truncated function. Enumerating all \(j\)-subsets of inputs and outputs would cost \(N^{O(d)}\), which is only quasipolynomial when \(d\) grows like \(\log(1/\varepsilon\delta)\).

   For each order \(j\), color every input column and every output mode uniformly at random with one of \(j\) colors. A fixed \(j\)-set receives distinct colors with probability \(j!/j^j\). Keeping only terms whose input set and output set are both colorful, and multiplying them by \((j^j/j!)^2\), leaves the expectation unchanged: \[\mathbb E_{\mathrm{col}}\bigl[\widehat q^{(1)}_{d,U}(\mathbf r)\bigr]=q_{d,U}(\mathbf r).\] The overhead factor obeys \(j^j/j!\le e^j\), so it is \(2^{O(j)}\).

   With colors fixed, the colorful sums are coefficients of an inverse determinant in formal variables with \(x_c^2=y_c^2=0\): \[\mathcal P_j=\det\Bigl[\delta_{uv}-x_{g(u)}z^{\mathbf 1\{u\in\text{prefix}\}}\sum_l y_{h(l)}B^{(l)}_{uv}\Bigr]^{-1},\] read off at \(x_1\cdots x_jy_1\cdots y_j z^p\). Since each \(x_c\) can appear at most once, the inverse is a finite polynomial, and the coefficient selects exactly one mode and one column of each color. Averaging many independent colorings controls the variance. Prefix queries then cost \(N^3\varepsilon^{-2}\delta^{-1}2^{O(d)}\) after preprocessing, with the exponent of \(N\) independent of \(d\).

   How to reuse it

   Color coding replaces a sum over all \(j\)-subsets by \(2^{O(j)}\) passes over colorful configurations, and reweighting keeps it unbiased. Use it whenever a low-degree expansion is accurate but the degree must grow slowly with the target accuracy. Pair it with marginal sampling to turn a signed approximation into a sampler.

Source: Sections III–VIII: Propositions 1–2, Lemmas 1–3, Theorem 2; Appendix B (color coding). Based on the linked paper version. Derivations marked as editorial, and suggested uses in other problems, are editorial guidance rather than results of the paper.

Context & limitations

Why it matters. Earlier provable simulations needed only \(O(\sqrt{\varepsilon N})\) photons to survive. Here, adding photons at a fixed transmission rate does not by itself escape efficient classical simulation.

Limitations. The guarantee is average-case over Haar-random interferometers. It needs many more modes than photons, \(M=\Omega(n^2N)\), and \(M=\Omega(n^2/(\varepsilon\delta))\) for full Fock outputs. It assumes uniform loss and ideal photons. Saturated-mode experiments, partial distinguishability and non-resolving detectors are left open.

Description based on the full text of the linked arXiv version, read 2026-10-21; equations follow the paper’s notation.