Exact simulation of Gaussian boson sampling in polynomial space and exponential time

We introduce an exact classical algorithm for simulating Gaussian boson sampling (GBS). The complexity of the algorithm is exponential in the number of photons detected, which is itself a random variable. For a fixed number of modes, the complexity is in fact equivalent to that of calculating output probabilities, up to constant prefactors. The simulation algorithm can be extended to other models such as GBS with threshold detectors, GBS with displacements, and sampling linear combinations of Gaussian states. In the specific case of encoding non-negative matrices into a GBS device, our method leads to an approximate sampling algorithm with polynomial runtime. We implement the algorithm, making the code publicly available as part of Xanadu's The Walrus library and benchmark its performance on GBS with random Haar interferometers and with encoded Erdős-Renyi graphs.

Figure 1

Running times of the exact simulation algorithm for a GBS system with 100 modes. The interferometer is Haar-random and the squeezing levels are set uniformly in the interval [0,1] for each mode. The algorithm is implemented in a cluster of 56 Intel(R) Xeon(R) CPUs operating at 2.6 GHz. The top points (circles) indicate the upper bound on the runtime, which occurs when $N$ photons are each observed in the first $N$ modes. The bottom points (stars) show a lower bound on the runtime, when $N$ photons are each observed in the last $N$ modes. Since we fix the detection pattern and also the total number of photons detected, the running time of these extreme cases is independent of the squeezing parameters used for the simulation.

Figure 2

Running times of the approximate simulation algorithm for non-negative kernel matrices. The kernel is the adjacency matrix of a random Erdős-Renyi graph with 100 vertices and edge probability $1/2$. The algorithm is implemented in a cluster of 56 Intel(R) Xeon(R) CPUs operating at 2.6 GHz. The top points (circles) show the upper bound on the runtime, which occurs when $N$ photons are each detected in the first $N$ modes. The bottom points (stars) indicate a lower bound on the runtime, when $N$ photons are each observed in the last $N$ modes.

Figure 3

(a) Bound on the total variation distance for a GBS setup as a function of the maximum number of photons $k$ that can be resolved in a detector. We assume 40 squeezed states with input mean photon number $\overline{n}=1$ are sent into a $100\times 100$ interferometer. The black dotted lines consider the lossless case and the gray lines consider applying 3 dB of loss in each of the modes. Note the exponential decay in the bound on the total variation distance. The results are averaged over 1000 unitaries drawn from the Haar measure and the bars denote one standard deviation. (b) We show the distribution of the total photon number $P\left(n\right)$ for pure and lossy GBS with the same parameters used in the top figure.