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\H{o}s-Renyi graphs.