Speeding up the classical simulation of Gaussian boson sampling with limited connectivity

Gaussian boson sampling (GBS) plays a crucially important role in demonstrating quantum advantage. As a major imperfection, the limited connectivity of the linear optical network weakens the quantum advantage result in recent experiments. In this work, we introduce an enhanced classical algorithm for simulating GBS processes with limited connectivity. It computes the loop Hafnian of an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n \times n$$\end{document}n×n symmetric matrix with bandwidth w in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(nw2^w)$$\end{document}O(nw2w) time. It is better than the previous fastest algorithm which runs in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(nw^2 2^w)$$\end{document}O(nw22w) time. This classical algorithm is helpful on clarifying how limited connectivity affects the computational complexity of GBS and tightening the boundary for achieving quantum advantage in the GBS problem.

The central issue in GBS experiments is to verify the quantum advantage of the result.The time cost of the best known classical algorithm for simulating an ideal GBS process grows exponentially with the system size [12].Therefore a quantum advantage result might be achieved when the system size is large enough [13][14][15].
However, there are always imperfections in real quantum setups, and hence the time cost of corresponding classical simulation will be reduced.When the quantum imperfection is too large, the corresponding classical simulation methods can work efficiently [16][17][18].In this situation, a quantum advantage result won't exist even if the system size of a GBS experiments is very large.Therefore, finding better methods to classically simulate the imperfect GBS process is rather useful in exploring the tight criteria for quantum advantage of a GBS experiment.
A major imperfection of recent GBS experiments is the shallow optical circuit [5,7,17,18].A shallow optical circuit has limited connectivity and its transform matrix will deviate from the global Haar-random unitary [18,19].In the original GBS protocol [2,3], the unitary transform matrix U of the passive linear optical network should be randomly chosen from Haar measure.However, due to the photon loss which increases exponentially with the depth of the optical network, the optical network might not be deep enough to meet the requirements of the full connectivity and the global Haarrandom unitary [19].Classical algorithms can take advantage of the bandwidth structure and the deviation from the global Haar-random unitary to realize a speedup in simulating the whole sampling process [17,18].Actually, with limited connectivity in the quantum device, the speed-up of corresponding classical simulation can be exponentially [18].
The most time-consuming part of simulating the GBS process with limited connectivity is to calculate the loop Hafnian of banded matrices.A classical algorithm to calculate the loop Hafnian of a banded n × n symmetric matrix with bandwidth w in time O (nw4 w ) is given in Ref. [17].Later, an algorithm that takes time O nw 2 2 w is given in Ref. [18].Here we present a classical algorithm to calculate the loop Hafnian of a banded n × n matrix with bandwidth w in time O (nw2 w ).We also show that this algorithm can be used to calculate the loop Hafnian of sparse matrices.
Our algorithm reduces the time needed for classically simulating the GBS process with limited connectivity.This is helpful in clarifying how limited connectivity affects the computational complexity of GBS and tightening the boundary of quantum advantage in the GBS problem.
This article is organized as follows.In Sec.II, we give an overview of the background knowledge which will be used later.In Sec.III, we give our improved classical algorithm for calculating the loop Hafnian of banded matrices.Finally, we make a summary in Sec.IV.

A. Gaussian Boson sampling protocol
In the GBS protocol, K single-mode squeezed states (SMSSs) are injected into an M -mode passive linear optical network and detected in each output mode by a photon number resolving detector.The detected photon number of each photon number resolving detector forms an output sample which can be denoted as n = n 1 n 2 ...n M .The schematic setup of Gaussian boson sampling is shown in Fig. 1 Before detection, the output quantum state of the passive linear optical network is a Gaussian quantum state [20][21][22][23].A Gaussian state is fully determined by its covariance matrix and displacement vector.Denote operator vector ξ = (â † 1 , . . ., â † M , â1 , . . ., âM ) T , where â † i and âi are the creation and annihilation operators in the ith (i ∈ {1, 2, . . ., M }) optical mode, respectively.The covariance matrix σ and displacement vector r of the Gaussian state ρ is defined as where Notice that ( ξ − r)( ξ − r) † is the outer product of column vector ( ξ − r) and row vector ( ξ − r) † and ( ξ − r)( ξ − r) † ̸ = (( ξ − r)( ξ − r) † ) † as the operators in the vector ξ might not commute with each other.As an example, assume M = 1 so that ξ = (â † 1 , â1 ) T , we have But, Usually, a matrix (say matrix A) is used to calculate the output probability distribution of a GBS process [2,3,14] as identity matrix with rank 2M (or M ).The matrix A is fully determined by the output Gaussian sate as follows: where The probability of generating an output sample n is where lhaf(A) = M ∈SPM(n) (i,j)∈M A i,j is the loop Hafnian function of matrix A, SPM is the single-pair matchings [12] which is the set of perfect matchings with loops.The matrix A n is obtained from A as follows n [14,24]: for ∀i = 1, . . ., M , if n i = 0, the rows and columns are deleted from matrix A; if n i ̸ = 0, rows and columns i and i + M are repeated n i times.

B. Continuous variables quantum systems
If a Gaussian state is input into a linear optical network, the output quantum state is also a Gaussian state.Denote the unitary operator corresponding to the passive linear optical network as Û .A property of the passive linear optical network is: where The output Gaussian state and the input Gaussian state is related by where σ in and rin are the covariance matrix and displacement vector of the input Gaussian state.The SMSS in input mode i has squeezing strength s i and phase ϕ i .
For simplicity, assume that ϕ i = 0 for i = 1, . . ., M .The covariance matrix of the input Gaussian state is The matrix Ã is thus Ã = B B * , where If a part of the optical modes in the Gaussian state are measured with photon number resolving detectors, the remaining quantum state will be a non-Gaussian state [25,26].A Gaussian state can be represented in the following form: β i and α i for i = 1, . . ., M are complex variables.For our convenience, define the permutation matrix P , such that γ Suppose the first N modes of an M -mode Gaussian state are measured and a sample pattern n = n where Rdd is a 2N ×2N matrix corresponding to modes 1 to represents the correlation between modes 1 to N and modes N + 1 to M .The remaining non-Gaussian state is [26] ρ where If all optical modes of the Gaussian state are measured (N=M), then Eq. ( 17) gives the probability of obtaining the output sample n, which is equivalent to Eq. (7).

C. Limited connectivity
In the Gaussian Boson Sampling protocol [2,3], the transform matrix U of the passive linear optical network should be randomly chosen from the Haar measure.The circuit depth needed to realize an arbitrary unitary transform in a passive linear optical network is O(M ), where M is the number of optical modes [19,[27][28][29].However, due to the photon loss, the circuit depth of the optical network might not be deep enough to meet the requirements of the full connectivity and the global Haarrandom unitary [19].This is because photon loss rate ϵ will increase exponentially with the depth of the optical network, i.e., ϵ = ϵ D 0 , where ϵ 0 is the photon loss of each layer of the optical network.If the photon loss rate is too high, the quantum advantage result of GBS experiments will be destroyed [16,30,31].
The shallow circuit depth lead to a limited connected interferometer [17][18][19].Assume that the beam splitters in the optical network is local, which means they can only connect the adjacent optical modes as shown in Fig. 1.If D < M , the transform matrix U of the passive linear optical network will have a bandwidth structure [17], i.e., and where w U = D is the bandwidth of the transform matrix U .An example of a matrix with bandwidth structure is shown in Fig. 2. According to Eq. ( 10) and (11), the covariance matrix σ of the output Gaussian state will have a bandwidth structure.Thus, if beam splitters are local, the circuit depth D must be no less than M/2 to reach the full connectivity.Note that, according to Eq. ( 12), matrix B (and A in Eq. ( 7)) has bandwidth Recently, a scheme known as "high-dimensional GBS" has been proposed [4].This scheme suggests that by interfering non-adjacent optical modes, the connectivity can be improved while maintaining a relatively shallow circuit depth.However, due to the limited circuit depth, the transformation matrix in this scheme cannot represent an arbitrary unitary matrix.Consequently, the transformation matrix deviates from the global Haarrandom unitary.To address this, the scheme introduces a local Haar-random unitary assumption, which says that the transformation matrix corresponding to each individual beam splitter is randomly selected from the Haar measure.Under this local Haar-random unitary assumption, Ref. [18] demonstrates that when the circuit depth is too shallow, the high dimensional GBS process can be approximate by a limited connected GBS process with a small error.
As a result of the shallow circuit depth and the deviation from the Haar measure, a speed-up can be realized in simulating the corresponding GBS process [17,18].The speed-up is attributed to the faster computation of the loop Hafnian for matrices with bandwidth as compared to the computation of general matrices.

D. Classical simulation of GBS
To date, the most efficient classical simulation method for simulating a general GBS process has been presented in Ref. [14].This classical algorithm, which samples from an M -modes Gaussian state ρ, operates as follows.
1.If ρ is a mixed state, it can be decomposed as a classical mixture of pure displaced Gaussian states.We randomly select a pure displaced Gaussian state based on this classical mixture.This can be done in polynomial time [14,21].Its covariance matrix and displacement vector are denoted as σ and r, respectively.
2. If ρ is a pure state, denote its covariance matrix and displacement vector as σ and r, respectively.

Compute the conditional covariance matrix σ
The most time-consuming part for the above classical simulation method is to calculate the probability p (n 1 , . . ., n k ) of an output sample pattern n 1 , . . ., n M .According to Eq. ( 7), we have where the r(B) A and A (B) A are the corresponding conditional matrices for subsystem contains modes 1 to k (denoted as A) when modes k + 1 to M are measured by heterodyne measurements with outcome denoted as A can be computed by the following equation: Note that if the Gaussian state ρ is a pure state, the conditional Gaussian state is still a pure state.So, we have Ã(B) As pointed in Ref. [18], we have BA = [U ( i tanh s i )U T ] A .This shows that BA has the same bandwidth structure as B. A classical algorithm to calculate the loop Hafnian of a banded matrix is thus needed for the classical simulation of GBS with limited connectivity.

III. SIMULATE THE SAMPLING PROCESS WITH LIMITED CONNECTIVITY A. Loop Hafnian algorithm for banded matrices
An algorithm to calculate the loop Hafnian of an n × n symmetric matrix with bandwidth w in time O (nw4 w ) is given in Ref. [17].Then, a faster algorithm with time complexity O nw 2 t 2 wt to calculate the loop Hafnian of an n×n symmetric matrix is proposed in Ref. [18], where w t is the treewidth of the graph corresponding to the matrix [32].For a banded matrix, the smallest treewidth w t is equal to the bandwidth w, i.e., w = w t .So the time complexity for this algorithm to calculate the loop Hafnian of an n × n symmetric matrix with bandwidth w is O nw 2 2 w .Here we show that the loop Hafnian of an n × n symmetric matrix with bandwidth w can be calculated in O (nw2 w ).
Our algorithm to calculate the loop Hafnian for banded matrices is outlined as follows.
Algorithm.To calculate the loop Hafnian of an n × n symmetric matrix B with bandwidth w: 2. Let t w = min (t + w, n), and P ({t + 1, . . ., t w }) be the set of all subsets of {t + 1, . . ., t w }.
3. For every Z t ∈ P ({t + 1, . . ., t w }) satisfying Z t ̸ = ∅ and |Z t | ≤ min (t, w), let and if t w ∈ Z t , then During the above iterations, if C t−1 {... } is not given in the previous steps, it is treated as 0.

Let
The loop Hafnian of matrix B is obtained in the final The time complexity of this algorithm is O(nw2 w ) as shown in theorem 1.
Theorem 1.Let B be an n × n symmetric matrix with bandwidth w.Then its loop Hafnian can be calculated in O(nw2 w ).
Proof.As shown in our algorithm, the number of coefficients (C t Zt , C t ∅ and C t {t} ) needed to be calculated for each t ∈ {1, . . ., n} is at most 2 w .As shown in Eq. ( 27) (28) (29), in each iteration, we need O(w) steps to calculate each coefficient C t Z t .So, for each t, the algorithm takes O(w2 w ) steps.The overall cost is thus O (nw2 w ).
An example of calculating a 4 × 4 matrix with bandwidth w = 1 using our algorithm can be found in Appendix A.
Note that, our algorithm for calculating loop Hafnian function of matrices with bandwidth structure can be easily extended to cases where the matrices is sparse (but not banded).A description of this can be found in Appendix B By combining our algorithm with the classical sampling techniques described in Ref. [14], as summarized in Sec.II D, a limited connected GBS process with a bandwidth w can be simulated in O(M nw2 w ) time, where M represents the number of optical modes and n denotes the maximum total photon number of the output samples.A proof of this statement is presented in Appendix C.

B. Validity of the algorithm
By sequentially computing the states (usually non-Gaussian) of the remaining M − k modes with the measurement outcome of modes 1 to k (k = 1, . . ., M ) in photon number basis, we can demonstrate the validity of the algorithm introduced in Sec.III A. Recalling matrix R defined in Eq. ( 15), we denote the bandwidth of matrix R by w.According to the definition of R and l, we know that lhaf(R) = lhaf(A), (31) where Assuming that the measurement results are n i = 1 for every i = 1, . . ., M , we obtain p(11 . . . 1) = P 0 lhaf(R). (33) Although we make the assumption that R corresponds to a Gaussian state, the subsequent proof remains valid in the more general case as discussed in Appendix D.
According to Eq. ( 17), after measuring the mode 1 in photon number basis, the remaining non-Gaussian quantum state is: We then have where D 2 x is the coefficient for γ2 x .
(38) Thus we prove This demonstrates the validity of the algorithm given in Sec.III A.
Intuitively, as shown in the previous derivation, after measuring the first t modes of the output Gaussian state, the state of the remaining modes is: where Poly t (R, γht ) represents the sum of polynomial terms formed by complex variables in γht and ) When all modes of the output state are measured, we have t = M for Eq.(41), and hence Poly M (R, γh M ) will not contain any complex variables, ensuring that the constant term C 2M ∅ equals to lhaf(R).Terms { x∈Z 2t γx } of different Z 2t influence the constant term C 2M ∅ in the subsequent computing.High-order terms containing γ k x with k ≥ 2 do not affect the final result.Consequently, the iteration of C t Z t appearing in step 3 of our algorithm is required from t = 1 to t = 2M to eventually determine the value of lhaf(R) = C 2M ∅ .

IV. SUMMARY
We present an algorithm to calculate the loop Hafnian of banded matrices.This algorithm runs in O (nw2 w ) for n × n symmetric matrices with bandwidth w.Our result is better than the prior art result of O nw 2 2 w [18].Our classical algorithm is helpful on clarifying how limited connectivity reduces the computational resources required for classically simulating GBS processes and tightening the boundary of quantum advantage in GBS problem.

(A1)
The algorithm given in Sec.III A works as follows: Step 1.

FIG. 1 .
FIG. 1.A schematic setup of Gaussian boson sampling.In this example, K = 8 single-mode squeezed states are injected into a passive linear optical network with M = 10 optical modes.Then photon number resolving detectors detect the photon number in each output mode.An output pattern n = n1n2...nM is generated according to the detected results.
. Denote that X 2M =   0 I M I M 0   , and I 2M (or I M )

FIG. 2 .
FIG.2.A symmetric matrix with a bandwidth structure.X represents an arbitrary non-zero entry in the matrix.The bandwidth in this example is w = 3.