View article

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_quantum_2020,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 th order approximation reproduces Ursell functions (also known as connected correlations) up to order with a cost exponential in 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 th order approximation would simply result in the uniform distribution, in contrast to GBS.