Quadratic speedup for simulating Gaussian boson sampling

Abstract

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.

The Boundary for Quantum Advantage in Gaussian Boson Sampling

“…This does not diminish the experimental achievement of large-scale GBS from Zhong et al [6], which remains faster than classical methods on supercomputers, if the time required for circuit programming (or in the present case fabricating a new, fixed interferometer) is not included. However, it has previously been reported that in boson sampling with Fock state inputs, at least 50 photon events are required to extend beyond the reach of an exact classical simulation in a reasonable time-scale [9]; for collision-free GBS, this threshold has been reported as being around 100 photons [8]; we have now demonstrated that for GBS with threshold detectors, the number of correlated detector clicks should also be around 100. For GBS with PNRDs, the number of photons required to surpass this classical threshold will depend on the amount of collisions, but must be ≥ 100.…”

“…B n is an N × N matrix, so it is considerably faster to calculate its loop hafnian compared to A n . While a realistic GBS experiment will not produce a pure state, a Gaussian mixed state can be expressed as a statistical ensemble of pure states with differing displacement vectors [8,10]; so for the purposes of a sampling algorithm, it is generally possible to randomly choose a complex displacement vector α from the correct distribution, then sample from the corresponding pure state. Hence the computational complexity of generating a sample is set by the calculation of an N × N loop hafnian, lhaf(B n ).…”

“…This significantly improves the numerical accuracy with only a minor time penalty. Whereas accuracy is an issue for eigenvalue-trace when N > 50 [8], the finitedifference sieve method has relative error < 10 −8 when tested up to N = 60 (Appendix B 4). This allows us to maintain accuracy for large loop hafnians while using a conventional 128-bit complex floating point data type, which is desirable for speed and portability.…”

“…These methods can be applied directly to the chainrule for simulating GBS described in [8] and Appendix A 2. Here, the photon number in each mode is sampled sequentially, conditioned on the photon numbers in the previous modes.…”

“…Futhermore, we introduce a new classical method to generate samples for GBS simulations with threshold detectors, which runs orders of magnitude faster than classical methods to generate samples with PNRDs, when collisions dominate. We apply these results to two sampling algorithms: a probability chain-rule method [8] and Metropolis independence sampling (MIS) [9]. We report nine orders of magnitude reduction in the time taken to simulate idealised Jiȗzhāngtype GBS experiments with threshold detectors.…”