Simulating Gaussian boson sampling on graphs in polynomial time

Abstract: We show that a distribution related to Gaussian Boson Sampling (GBS) on graphs can be sampled classically in polynomial time. Graphical applications of GBS typically sample from this distribution, and thus quantum algorithms do not provide exponential speedup for these applications. We also show that another distribution related to Boson sampling can be sampled classically in polynomial time.

• The paper demonstrates that any GBS-based sampling distribution on graphs can be simulated classically in polynomial time.
• It introduces a graph transformation using the Cartesian product and specific edge weighting to reduce the problem to sampling perfect matchings.
• The approach shows that quantum algorithms for graph-related GBS tasks yield at most a polynomial advantage over classical methods.

Polynomial-Time Classical Simulation of Gaussian Boson Sampling on Graphs

Introduction

This work resolves a central question in the interface between quantum algorithms and classical complexity for sampling problems linked to bosonic quantum systems. Specifically, it investigates the classical simulability of Gaussian boson sampling (GBS) on graphs, a task that had been widely conjectured to exhibit exponential quantum speedup for various graph-theoretic applications. The authors prove any GBS-based sampling distribution associated with graph parameters can be realized by a polynomial-time classical algorithm. As a consequence, quantum algorithms based on GBS do not yield exponential advantages for these tasks. This result also extends to a boson sampling variant under non-negative matrix inputs.

Background and Problem Formulation

Boson sampling, originally proposed by Aaronson and Arkhipov, involves simulating the output distribution of indistinguishable non-interacting bosons passing through a linear optical network, where output probabilities are proportional to the squared permanent of submatrices of a unitary transition matrix. GBS generalizes this by using Gaussian input states and hafnians rather than permanents, and has become prominent due to experimental accessibility and its connections to graph invariants.

For graph applications, the relevant GBS distribution for a graph $G = (V, E)$ and subset $S \subseteq V$ is

$\mu_{GBS,G}(S) \propto c^{2|S|} PM(G[S])^2,$

where $PM(G[S])$ is the number of perfect matchings of the induced subgraph $G[S]$, and $c$ is a normalization derived from the spectral properties of $G$.

Previous algorithms for classical sampling from this distribution were only efficient for extremely dense graphs or relied on coarser variants of the distribution (e.g., sampling proportional to $PM(S)$ instead of $PM(S)^2$). The general case for arbitrary graphs was open.

Main Algorithms and Theoretical Contributions

Classical Polynomial-Time GBS Simulation

The main construction leverages the connection between GBS distributions and weighted perfect matchings in carefully constructed graphs. The algorithm proceeds as follows:

1. Graph Transformation: Given a graph $G = (V, E)$, construct the Cartesian product $G \square K_2$, incorporating two copies of $G$ and additional matching edges between them.
2. Edge Weighting: Assign weights to edges from $G$ and its copy by $c^2$, and to cross edges by $1$.
3. Reduction to Matching Problem: Demonstrate that sampling perfect matchings from the weighted $G \square K_2$ and projecting to $V$ yields samples from the desired GBS graph distribution.
4. Application of Jerrum–Sinclair Markov Chain: Use the Markov chain algorithm of Jerrum and Sinclair for sampling (weighted) perfect matchings. The analysis confirms that the required ratio of near-perfect to perfect matchings is polynomially bounded, ensuring polynomial mixing and sampling time.

Formally, given any graph and precision parameter $\epsilon$, the algorithm samples from a distribution $\epsilon$-close in total variation distance to $\mu_{GBS,G}$ in time $$O(\overline{c}mn^4 \log^2 (n\overline{c}/\epsilon))$$.

Extension to Non-negative Boson Sampling

A similar approach addresses the simulation of boson sampling when the input matrix $A$ is non-negative, i.e., not the general complex case typical for Haar-random unitaries. The reduction constructs a bipartite graph whose perfect matchings correspond to the output values of the required boson sampling distribution. The algorithm uses the Jerrum–Sinclair–Vigoda method to sample in polynomial time. The runtime is $$O\left(\frac{m^7n^{14}}{\epsilon^7}\log^4\left(\frac{mn}{\epsilon}\right)\right)$$ for input matrix size $m \times n$ and accuracy $\epsilon$.

Complexity Separation and Limitation

The result only applies to non-negative (and, for GBS, graph-related) input instances. For generic boson sampling with general complex weights, the hardness stems from computing or approximating the permanent of complex matrices—known to be $\#P$-hard even for positive semi-definite instances. Thus, the classical tractability achieved here leverages structural restrictions on the matrices involved, not genericity.

Implications and Theoretical Significance

• No Exponential Quantum Speedup for GBS-Based Graph Applications: For any graph-based application reducible to sampling from the GBS-related distribution (including key combinatorial optimization and isomorphism tasks), quantum algorithms via GBS offer at most polynomial advantage.
• Bridging Quantum and Classical Complexity: The results provide an explicit classical algorithm that matches, in polynomial time, the sampling capabilities of photonic quantum computers for a widely studied family of quantum distributions—contradicting the presumption of quantum supremacy for GBS in these domains.
• Algorithmic Approach: The main technical novelty is the combinatorial reduction from GBS and (restricted) boson sampling to perfect matching sampling in augmented graphs. This reduction exploits the structure of hafnians and permanents arising from graph adjacency matrices and non-negative matrices, respectively.
• Separation from General Boson Sampling: The work clarifies the boundary—general boson sampling with arbitrary complex weights remains out of reach for classical simulators due to the hardness of approximating complex permanents, a key circuit complexity bottleneck. However, graph-theoretic and non-negative cases fall to efficient classical techniques.

Future Directions

• Extension to Structured Non-Negative Instances: Whether the classical algorithms can be generalized to larger families of inputs, such as real positive semi-definite or structured complex-valued matrices, is open.
• Conditional Sampling and Permanent Approximation Oracles: The reductions show that, given an oracle for (approximate) permanent computation, broader sampling is accessible. Advances in permanent approximation could hence lift the horizon for classical simulation.
• Quantum-Inspired Algorithms for Other Physical Models: The paradigm here suggests possible quantum-inspired reductions for other sampling tasks in quantum optics or many-body systems, leveraging Markov chain Monte Carlo techniques.

Conclusion

This paper establishes that sampling from GBS distributions on graphs, as well as boson sampling with non-negative matrices, is tractable in polynomial time by classical algorithms. This fundamentally limits the scope of quantum speedup achievable by current photonic architectures for a large class of graph-driven tasks and clarifies the complexity landscape for quantum vs. classical simulation of bosonic systems (2511.16558).