Efficient approximation of experimental Gaussian boson sampling

Two recent landmark experiments have performed Gaussian boson sampling (GBS) with a non-programmable linear interferometer and threshold detectors on up to 144 output modes see Refs. onlinecite{zhong_quantum2020,zhong2021phase}. Here we give classical sampling algorithms with better total variation distance and Kullback-Leibler divergence than these experiments and a computational cost quadratic in the number of modes. Our method samples from a distribution that approximates the single-mode and two-mode ideal marginals of the given Gaussian boson sampler, which are calculated efficiently. One implementation sets the parameters of a Boltzmann machine from the calculated marginals using a mean field solution. This is a 2nd order approximation, with the uniform and thermal approximations corresponding to the 0th and 1st order, respectively. The kth order approximation reproduces Ursell functions (also known as connected correlations) up to order k with a cost exponential in k and high precision, while the experiment exhibits higher order Ursell functions with lower precision. This methodology, like other polynomial approximations introduced previously, does not apply to random circuit sampling because the kth order approximation would simply result in the uniform distribution, in contrast to GBS.