Source-linked AI summary
The Boundary for Quantum Advantage in Gaussian Boson Sampling
Jacob F. F. Bulmer, Bryn A. Bell, Rachel S. Chadwick, Alex E. Jones, Diana Moise, Alessandro Rigazzi, Jan Thorbecke, Utz-Uwe Haus, Thomas Van Vaerenbergh, Raj B. Patel, Ian A. Walmsley, Anthony Laing
TL;DR
The paper addresses limited understanding of classical simulation complexity for GBS with high collision rates, a regime relevant to threshold-detector experiments. It develops faster probability and sampling methods and a classically tractable distribution. The methods reduce estimated simulation time by nine orders of magnitude for idealised Jiuzhang-type experiments, while the authors identify photon-number-resolving detectors as a route to greater complexity.
Problem
The classical complexity of GBS in regimes with high collision rates is poorly understood, including for Jiuzhang, which uses threshold detectors in a high-collision regime.
Source: §I · p. 1Method
The paper improves collision-aware loop-hafnian calculations, develops faster threshold-detector sampling, and introduces the Independent Pairs and Singles distribution.
Source: §I · p. 1, §II · p. 2, §VI · p. 5Results
Nine orders of magnitude reduction in estimated simulation time for idealised Jiuzhang-type GBS experiments with threshold detectors.
Source: §I · p. 1Takeaways & Limitations
Independent Pairs and Singles passes a variety of GBS verification tests, suggesting it should be used as a classically simulable adversary in stronger verification.
Source: §VI · p. 5A 100-mode PNRD simulation was impractically slow on the benchmarking system; generating one 92-photon event took 82 minutes.
Source: §VI · p. 6
Abstract
from arXiv · showhide
Identifying the boundary beyond which quantum machines provide a computational advantage over their classical counterparts is a crucial step in charting their usefulness. Gaussian Boson Sampling (GBS), in which photons are measured from a highly entangled Gaussian state, is a leading approach in pursuing quantum advantage. State-of-the-art quantum photonics experiments that, once programmed, run in minutes, would require 600 million years to simulate using the best pre-existing classical algorithms. Here, we present substantially faster classical GBS simulation methods, including speed and accuracy improvements to the calculation of loop hafnians, the matrix function at the heart of GBS. We test these on a $\sim \! 100,000$ core supercomputer to emulate a range of different GBS experiments with up to 100 modes and up to 92 photons. This reduces the run-time of classically simulating state-of-the-art GBS experiments to several months -- a nine orders of magnitude improvement over previous estimates. Finally, we introduce a distribution that is efficient to sample from classically and that passes a variety of GBS validation methods, providing an important adversary for future experiments to test against.
I. INTRODUCTION
The classical complexity of GBS with high collision rates and threshold detectors remains poorly understood, complicating efforts to locate the boundary for quantum advantage. This paper develops faster classical simulation methods and identifies a classically tractable distribution that passes several GBS verification tests.
I. INTRODUCTION
- Jiuzhang ran in 200 s, while the best available classical algorithms were estimated to need 600 million years to simulate it.The experiment comprised 100 optical modes and observed up to 76 photon detection events.Source: §I · p. 1
- The classical complexity of GBS in regimes with high collision rates is poorly understood, even though Jiuzhang uses threshold detectors and operates where collisions are highly probable.The original proposal treated collisions as improbable and neglected them.Source: §I · p. 1
- The paper reports a nine-orders-of-magnitude reduction in estimated simulation time for idealised Jiuzhang-type threshold-detector experiments.The methods enabled simulations on a roughly 100,000-core supercomputer, including 100 modes and up to 60 click events.Source: §I · p. 1
- The authors introduce a classically tractable distribution that passes a variety of canonical GBS verification tests.They argue this sharpens the boundary for quantum advantage and highlights the need for stringent adversarial tests.Source: §I · p. 2
II. LOOP HAFNIAN ALGORITHMS
The paper accelerates loop-hafnian calculations by exploiting repeated photon pairings in collision-heavy outcomes, and presents a finite-difference sieve to improve numerical accuracy. These methods support faster probability calculations for GBS with photon-number-resolving detectors.
II. LOOP HAFNIAN ALGORITHMS
- The loop hafnian of a matrix derived from the photon pattern and Gaussian state determines the probability of a GBS detection event.For pure Gaussian states, the matrix is block diagonal, allowing calculation using a smaller matrix; sampling mixed states can use randomly chosen displacement vectors.Source: §II · p. 2
- Grouping photons into repeated pairs reduces the inclusion/exclusion term count to ∏_i(η_i + 1), between ∏_i√(n_i + 1) and 2^(N/2).A greedy algorithm selects pairings to approach the lower bound; the figure illustrates repeated pairs and the resulting speed-up.Source: §II · p. 2, Figure 1 · p. 2
- For n = (2, 1, 0, 3), the repeated-pairing calculation uses 6 terms, versus 8 for eigenvalue-trace and 24 for repeated-moment.One of the two photon pairings is repeated.Source: §II · p. 2
- The finite-difference sieve has relative error < 10^-8 through N = 60, addressing eigenvalue-trace accuracy issues above N = 50 with only a minor time penalty.The authors use this formula for all benchmarking results, with conventional 128-bit complex floating point.Source: §II · p. 2
III. THRESHOLD DETECTORS
For threshold-detector sampling, the paper avoids calculating the full Torontonian probability sum by modeling detection through a fan-out of sub-detectors. The resulting sampling calculation is based on an N_c × N_c loop hafnian.
III. THRESHOLD DETECTORS
- Threshold-detector probabilities can be calculated using a Torontonian sum with 2^N_c terms, but that need not be the fastest route to sampling.N_c is the number of clicks, or outputs with one or more photons.Source: §III · pp. 2–3
- The proposed sampling approach reduces the calculation to an N_c × N_c loop hafnian with 2^(N_c/2) terms, for any photon density.At sufficiently low photon density, simulating PNRDs and converting nonzero photon counts to clicks may be faster than direct threshold-probability calculation.Source: §III · p. 3
- The fan-out model samples sub-detectors sequentially and stops once a photon is detected, because later sub-detectors in that mode can be ignored.The detected photon's normalized position x corresponds to a fraction x of sub-detectors ignored, equivalent to loss before single-photon detection.Source: §III · p. 3
A. Chain-rule sampling
The chain-rule method samples mode photon counts sequentially, conditioning each mode on the preceding counts. Batching reuses calculations across candidate counts and sub-detectors, while the figure compares MIS estimates with exact outcome probabilities.
A. Chain-rule sampling
- For each mode, chain-rule sampling calculates joint probabilities over photon counts up to a cutoff, conditioned on counts in earlier modes.The cutoff is chosen so that the probability of more photons can be neglected, and is generally several times the expected photon number.Source: A. Chain-rule sampling · p. 3
- Batching calculates loop hafnians for candidate counts up to n_cut with approximately the runtime of the largest calculation.For threshold detectors, n_cut is reduced to 1 per sub-detector, and batching across sub-detectors largely offsets their added sampling overhead.Source: A. Chain-rule sampling · p. 3
- Figure 2 compares MIS-estimated probabilities with exact probabilities for all 6-photon outcomes in 8-mode PNRD and 8-mode, 3-click threshold-detector simulations.Blue bars denote MIS estimates and orange bars denote exact probabilities.Source: Figure 2 · p. 3
B. Metropolis Independence Sampling
The section describes Metropolis independence sampling (MIS) for GBS and an IPS-based proposal distribution, then extends the sampling representation for threshold detectors.
B. Metropolis Independence Sampling
- MIS proposes samples from a distribution Q and accepts them with a probability based on the target P and proposal Q distributions.Rejected proposals repeat the previous sample; burn-in and thinning affect efficiency.Source: B. Metropolis Independence Sampling · p. 3
- The proposal pairs a displacement vector drawn for the mixed state with a photon pattern sampled from the Independent Pairs and Singles distribution.IPS is efficiently sampleable and assigns probabilities using N × N loop hafnians of positive matrices.Source: B. Metropolis Independence Sampling · pp. 3–4
- IPS alone is sufficient to pass many GBS verification methods, and its samples supply proposals for MIS.Source: B. Metropolis Independence Sampling · p. 4
- For threshold-detector MIS, the method adds a continuous variable for the first detected photon and calculates probabilities with Nc × Nc loop hafnians.Tracing out unused sub-detectors reintroduces mixture, so the target-probability calculation uses an adjusted displacement vector.Source: B. Metropolis Independence Sampling · p. 4
V. BENCHMARKING
The benchmarks test the simulation methods on small distributions and larger GBS experiments. They report substantial collision-related speedups and runtime estimates for chain-rule sampling and MIS, while identifying chain-rule sampling as preferable when independent samples are required.
V. BENCHMARKING
- The benchmarks use Haar-random interferometers, M/4 two-mode squeezed-vacuum sources, squeezing r = 1.55, and transmission η = 0.3, with parameters similar to Zhong et al..Source: §V · p. 4
- ∼10 s per sample is the estimated average above 80 photons for the M = 60 PNRD chain-rule simulation; Fugaku could reduce this to 130 ms.The estimate extrapolates from the measured runtime fit.Source: Figure 4 · p. 4, §V · p. 5
- For M = 8, MIS TVDs were 0.0153 for PNRD and 2.9 × 10^-3 for the Nc = 3 threshold distribution; chain-rule TVDs were 0.0554 and 0.0138, respectively.The chain-rule post-selected distributions had fewer samples, which explains their larger TVDs.Source: §V · p. 4
- 8.4 s per sample is the predicted mean chain-rule runtime for an M = 100 threshold-detector experiment; MIS averages 0.45 s per step, but independent samples favor chain-rule sampling.On Fugaku, the estimates are around 127 ms per chain-rule sample and 7 ms per MIS step; MIS thinning is likely to exceed 600.Source: Figure 5 · p. 5, §V · p. 5
- Collision-aware loop-hafnian calculations generally run one to two orders of magnitude faster, making an 80-photon probability comparable in time to a 60-photon probability without collision speed-up.Runtime varies substantially between samples with the same photon number, depending on their collision pattern.Source: §V · p. 5
VI. CONCLUSION
The conclusion sharpens estimates of classical GBS simulation runtimes and identifies photon and detector-click scales relevant to the classical threshold. It also presents a polynomial-time classical adversary for verification and states the scope of the near-exact simulations.
VI. CONCLUSION
- Independent Pairs and Singles (IPS) samples in polynomial time and largely passes quantitative GBS tests, making it a better approximation than standard adversarial models.The authors suggest using IPS as a classically simulable adversary and call for stronger verification methods.Source: §VI · p. 5
- Quadratic complexity reduction brings the 100-mode chain-rule simulation of Zhong et al.’s 50 million samples in 200 s to ~73 days; threshold GBS needs around 100 correlated clicks to cross the classical threshold.For PNRD GBS, the photon threshold depends on collisions and must be at least 100.Source: §VI · pp. 5–6
- These near-exact simulations exploit no experimental imperfections beyond collisions; faster classical methods may exploit photon loss, distinguishability, or limits on Haar-random transformations.Source: §VI · p. 6
Appendix A: Gaussian Boson Sampling
- An M-mode Gaussian state’s Wigner representation uses a 2M-component mean vector and a 2M × 2M covariance matrix; PNRD pattern probabilities are given using loop hafnians.Source: Appendix A · p. 6
- For pure states, probabilities can use loop hafnians of half-size matrices; B_n is symmetric and N × N, where N is the total photon number.Source: Appendix A · p. 6
1. Sampling pure Gaussian states from mixed Gaussian states
- Williamson decomposition expresses a mixed-state covariance as pure-state covariance plus Gaussian noise; sampling draws a pure state’s mean from that noise distribution.Source: §1 · pp. 6–7
2. Chain-rule GBS sampler
The chain-rule sampler builds GBS outcomes sequentially, while specialized loop-hafnian algorithms target repeated pairings, numerical accuracy, and batched probability calculations. Together these methods improve the efficiency and accuracy of classical GBS simulation.
2. Chain-rule GBS sampler
- The chain-rule sampler replaces coherent-state amplitudes with photon numbers sequentially, approximating each conditional distribution by normalized joint probabilities up to a cutoff.The cutoff is chosen so the probability of exceeding it can be neglected; these joint-probability calculations dominate sampling effort.Source: §2 · p. 7
1. Eigenvalue-trace
- Eigenvalue-trace reduces trace-of-matrix-power calculations to eigenvalues of CX and performs inclusion/exclusion over a fixed pairing.Its cost includes 2^(N/2) terms, each dominated by finding eigenvalues in O(N^3) time.Source: §1 · p. 8
2. Repeated pairs
- For mixed-state matrices with repeated pairings, lhafmix sums over multiplicities zi from 0 to ni rather than all pair-inclusion choices.Pairs with zero multiplicity can be deleted to speed up the eigenvalue calculation.Source: §2 · p. 8
3. Matching algorithm
- A greedy matching algorithm chooses fixed pairings to minimize inclusion/exclusion steps, improving on eigenvalue-trace when repeated pairs occur.The stated improvement conditions are at least two photon-number entries greater than 1, or one entry at least 4.Source: §3 · p. 9
4. Finite difference sieve
- The finite-difference sieve offers significantly better numerical accuracy than inclusion/exclusion by dividing accumulated error by the number of sum terms.The comparison uses random block-diagonal constructions and is shown in Fig. 6.Source: §4 · p. 9, Figure 6 · p. 9
5. Batching probability calculations
- Batching computes probabilities for every outcome up to ncut with the eigenvalue work needed for the largest outcome in the batch.The same approach batches calculations across sub-detectors within one threshold-detector mode.Source: §5 · pp. 9–10
6. Implementation details
The implementation uses compiled, parallel loop-hafnian code and develops MIS proposals for photon-number-resolving and threshold-detector GBS. The proposals use efficient sampling constructions and express threshold-click probabilities through loop hafnians.
6. Implementation details
- The loop-hafnian code uses Python with Numba and MPI, and was tested across major operating systems and x86-64 and arm64 architectures.The Fugaku and Isambard systems used for runtime predictions have arm64 chips.Source: §6 · p. 10
1. Independent Pairs and Singles GBS distribution
- Independent Pairs and Singles (IPS) proposes samples from independent single-photon and photon-pair processes without interference between them.It is a better approximation to GBS than thermal states or distinguishable squeezed states; evaluating both loop hafnians adds at most a factor of 2 to runtime.Source: §1 · p. 10
2. PNRD GBS
- For PNRD MIS, the sampler augments photon-number patterns with a displacement variable and chooses its proposal distribution to match the target displacement distribution.This lets the acceptance probability avoid dependence on the displacement probability density.Source: §2 · p. 10
- Post-selecting MIS proposals on a fixed total photon number does not appreciably add runtime, because proposal generation is efficient and acceptance calculations dominate.The acceptance probability remains identical when fixing the photon number.Source: §2 · p. 11
3. Threshold detector GBS
- For threshold-detector MIS, click patterns with loss variables can be expressed as single-photon patterns whose probabilities use an Nc × Nc loop hafnian.The resulting calculation has a quadratic speedup compared with a Torontonian.Source: §3 · p. 11
4. Thinning interval and burn-in time scaling
The section estimates how thinning intervals and burn-in times scale with mode count, while emphasizing that the estimates depend on the desired closeness to the target distribution. It also describes CHOG as a likelihood-based comparison of trial and adversarial samples against an ideal distribution.
4. Thinning interval and burn-in time scaling
- τthin = 6.81M −65.3 estimates the interval needed to reduce repeat probability to 0.1; the M = 36 point is anomalous, likely from an extreme target/proposal ratio.Source: Figure 8 · p. 13, §4 · p. 13
- The likelihood test's fitted burn-in relation is τburn = 1.43M + 12.86, extrapolating to τburn = 155 for 100 modes.Source: Figure 11 · p. 14, §4 · p. 14
- The burn-in tests yield significantly different estimates; finite-sample noise limits both, and the chosen convergence threshold is arbitrary.Source: §4 · pp. 14–15
- The acceptance-rate test gives a different fitted relation, τburn = 8.03M −17.06, using the point where acceptance is within 0.001 of its estimated final value.Source: §4 · p. 14, Figure 13 · p. 15, Figure 12 · p. 15
- Burn-in time scales approximately linearly with mode count, and its fitted gradient depends on how closely the sampled distribution must match the target.Source: §4 · p. 15
1. CHOG ratio
- CHOG compares trial and adversarial samples by their likelihood under an ideal distribution; convergence to 1 favors the trial samples.Source: §1 · pp. 15–16
a. IPS vs thermal
An efficiently classically sampled IPS distribution passes the CHOG test against thermal samples, challenging that test's usefulness for validating GBS complexity against a thermal adversary. The authors propose IPS and finite-purity distributions as adversaries for future tests.
- IPS samples converge to a CHOG ratio of 1 against thermal samples, indicating they are more likely than thermal samples to come from the ideal distribution.Source: a. IPS vs thermal · p. 16
- Because IPS can be sampled efficiently classically yet passes this test, the authors argue it challenges CHOG's usefulness against thermal adversaries and recommend IPS as a future adversary model.Source: a. IPS vs thermal · p. 16
- The authors also propose finite-purity adversaries, with squeezing across at least two Schmidt modes, as a more challenging model to test.Source: a. IPS vs thermal · p. 16
b. MIS vs IPS
The MIS chain is tested against its IPS proposal distribution using CHOG, while two-point correlators provide a separate comparison among ideal, thermal, and IPS distributions.
b. MIS vs IPS
- The comparison uses a 24-mode interferometer with six squeezed-vacuum sources, 10^5 IPS proposals, burn-in 50, and thinning interval 10.Source: b. MIS vs IPS · p. 16
- MIS samples converge to a CHOG ratio of 1 against IPS samples, indicating they are more likely to have come from the ideal distribution than the starting IPS samples.Source: b. MIS vs IPS · p. 16
2. Two-point correlators
- At 32 modes, ideal and IPS two-point correlator values overlap, while both differ significantly from the thermal distribution.Source: Figure 16 · p. 17