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 just 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 employ the Williamson decomposition to express a mixed Gaussian state as a probability distribution over pure displaced Gaussian states. By first sampling this random displacement, it is possible to reduce 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.