  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 $k$th
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
$k$th order approximation would simply result in the uniform distribution, in
contrast to GBS.
