Quadratic speedup for simulating Gaussian boson sampling

We introduce an algorithm for the classical simulation of Gaussian boson sampling that is quadratically faster than previously known methods. The complexity of the algorithm is exponential in the number of photon pairs detected, not the number of photons, and is directly proportional to the time required to calculate a probability amplitude for a pure Gaussian state. The main innovation is to use auxiliary conditioning variables to reduce the problem of sampling to computing pure-state probability amplitudes, for which the most computationally-expensive step is calculating a loop hafnian. We implement and benchmark an improved loop hafnian algorithm and show that it can be used to compute pure-state probabilities, the dominant step in the sampling algorithm, of up to 50-photon events in a single workstation, i.e., without the need of a supercomputer.