Quadratic Speed-Up 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 the computation of the pure-state probability amplitudes, for which the most computationally expensive step is the calculation of 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 events involving up to 50 photons in a single workstation, i.e., without the need of a supercomputer.