Simulating lossy Gaussian boson sampling with matrix-product operators

Gaussian boson sampling, a computational model that is widely believed to admit quantum supremacy, has already been experimentally demonstrated and is claimed to surpass the classical simulation capabilities of even the most powerful supercomputers today. However, whether the current approach limited by photon loss and noise in such experiments prescribes a scalable path to quantum advantage is an open question. To understand the effect of photon loss on the scalability of Gaussian boson sampling, we analytically derive the asymptotic operator entanglement entropy scaling, which relates to the simulation complexity. As a result, we observe that efficient tensor network simulations are likely possible under the $N_{\text{out}}\propto \sqrt{N}$ scaling of the number of surviving photons $N_{\text{out}}$ in the number of input photons $N$. We numerically verify this result using a tensor network algorithm with $\mathrm{U}\left(1\right)$ symmetry, and we overcome previous challenges due to the large local Hilbert-space dimensions in Gaussian boson sampling with hardware acceleration. Additionally, we observe that increasing the photon number through larger squeezing does not increase the entanglement entropy significantly. Finally, we numerically find the bond dimension necessary for fixed accuracy simulations, providing more direct evidence for the complexity of tensor networks.

Figure 1

Operator entanglement entropy vs the number of input squeezed modes for different photon survival scaling $N_{\text{out}}=\beta N^{\gamma }$ at $r=0.88$. (a) $\gamma =\frac{1}{4}$. (b) $\gamma =\frac{1}{2}$. (c) $\gamma =1$. (d) Convergence of MPO EE with increasing $n_{\text{max}}$ for $N=50,\beta =1,\gamma =\frac{1}{2},r=0.88$.

Figure 2

Operator entanglement entropy vs the number of input squeezed modes for different photon survival scaling $N_{\text{out}}=\beta N^{\gamma }$ at $r=0.88$. Details of experiment configurations can be found in Methods. (a) $\gamma =\frac{1}{4}$. (b) $\gamma =\frac{1}{2}$. (c) $\gamma =1$. Dots are results obtained from full simulations using $\mathrm{U}\left(1\right)$ symmetry. Dashed lines are estimates using asymptotic assumptions.

Figure 3

Operator entanglement entropy vs the number of input squeezed modes for different squeezing parameters $r$. Dashed lines are guides to the eye. (a) $r=0.88$, averaging approximately one photon per mode. (b) $r=1.146$, averaging approximately two photons per mode. (b) $r=1.44$, averaging approximately four photons per mode.

Figure 4

Analysis of bond dimension, system size, and error. Details can be found in Methods. (a) Bond dimension needed to reach accuracy $1-\text{Tr}\left(\widehat{\rho }\right)=0.02$ vs the number of input squeezed modes' photon survival scaling $N_{\text{out}}=0.4N^{\gamma }$ at $r=0.88$. Dots are individual estimates of the bond dimension obtained from full simulations using $\mathrm{U}\left(1\right)$ symmetry. Dashed lines are the means. Traces have higher estimated bond dimensions in ascending order in $\gamma$. (b) Reduction in $1-\text{Tr}\left(\widehat{\rho }\right)$ error as bond dimension increases for three different experimental configurations.

Figure 5

Algorithms for computing $\mathrm{\Theta }$ matrices. (a) CPU-based implementation. A subset of bonds is selected from ${\mathrm{\Gamma }}^{\left[k\right]},{\mathrm{\Gamma }}^{[k+1]}$ that have the correct selected charge values $c^{[k-1]},c^{[k+1]}$. A subset of $\mathrm{\Theta }$ is computed. (b) GPU-based implementation. All bonds are used and the entire $\mathrm{\Theta }$ matrix is computed at once.

Figure 6

Matrix multiplication with (a) inner products and (b) outer products.

Figure 7

Memory access pattern (a) without sorting and (b) with sorting.

Figure 8

Illustration of insertion of empty bonds. (a) Worst case scenario of charge value changes within a single fragment without empty bond insertion. The thread has to store 16 values of the unitary. (b) Generic case of a fragment at a charge change boundary. Less than 16 values need to be stored, but this is not known a priori and 16 values of the unitary still need to be stored. (c) With empty bond insertion, each thread only needs one unitary value.

Figure 9

High-level parallelization. Independent unitary gate updates are distributed to different nodes. Within each unitary update, $\mathrm{\Theta }\left(c^{\left[k\right]}\right)$'s are computed and decomposed with SVD independently on different GPUs.

Figure 10

(a) CPU and GPU simulation time. Traces have higher simulation time in ascending order in $d$. (b) Contribution to the CPU simulation time from subroutines.

Figure 11

Contribution to simulation time of the GPU algorithm from subroutines. (a) Bond dimension $\chi =4096$. (b) Local Hilbert-space size $d=14$.