Approximating outcome probabilities of linear optical circuits

Quasiprobability representation is an important tool for analyzing a quantum system, such as a quantum state or a quantum circuit. In this work, we propose classical algorithms specialized for approximating outcome probabilities of a linear optical circuit using $s$-parameterized quasiprobability distributions. Notably, we can reduce the negativity bound of a circuit from exponential to at most polynomial for specific cases by modulating the shapes of quasiprobability distributions thanks to the norm-preserving property of a linear optical transformation. Consequently, our scheme renders an efficient estimation of outcome probabilities with precision depending on the classicality of the circuit. Surprisingly, when the classicality is high enough, we reach a polynomial-time estimation algorithm within a multiplicative error. Our results provide quantum-inspired algorithms for approximating various matrix functions beating best-known results. Moreover, we give sufficient conditions for the classical simulability of Gaussian boson sampling using the approximating algorithm for any (marginal) outcome probability under the poly-sparse condition. Our study sheds light on the power of linear optics, providing plenty of quantum-inspired algorithms for problems in computational complexity.


I. INTRODUCTION
Quantum computers are believed to provide significant advantages in solving computational problems beyond classical power [1,2].Despite the potential advantage, determining which types of quantum circuits are classically simulable is still a longstanding open problem.A viable approach for the problem is to examine quasiprobability distributions in the phase space, a standard method in quantum optics [3].Specifically, on the one hand, if the output distributions of a quantum circuit can be described by nonnegative quasiprobability distributions, there is a classically efficient simulation of the circuit [4].On the other hand, one can estimate the outcome probabilities of a quantum circuit with the convergence rate depending on the negativity bound of the circuit, a measure of negativity in the phase space [5].In both cases, the negativity of quasiprobability distributions plays a crucial role, which has been extensively studied [6][7][8][9].
Meanwhile, boson sampling has recently attracted lots of attention due to its feasible quantum advantage using a linear optical circuit [10][11][12][13][14]. Also, there have been numerous studies of boson sampling using quasiprobability distributions [15][16][17][18][19].In particular, while Fock-state and Gaussian boson sampling are believed to be hard to classically simulate [10,20], boson sampling with a classical state input, represented by a nonnegative Pfunction, is efficiently simulated because its output distribution can be expressed by nonnegative quasiprobability distributions [16].Interestingly, this quantum optical observation leads to the fact that the permanent of Hermitian positive-semidefinite (HPSD) matrices can be approximated within a multiplicative-error more easily (in BPP NP ) than that of arbitrary matrices (#Phard) [15,21].In addition, a quantum-inspired algorithm using quasiprobability distributions has been proposed for estimating the permanent of an HPSD matrix within an additive-error, which outperforms Gurvits' algorithm under certain eigenvalue conditions [22,23].Thus, studying the quasiprobability representation of quantum circuits often provides new insight into computational problems.
However, previous studies are limited to the cases of classical input states with no negativity [15,16,18,23] or negativity bound that is polynomial in the system size [5], which cannot cover typical quantum circuits with exponential negativity bound.A central question is how to generalize the quasiprobability methods to handle cases of exponential negativity bound.Such a generalization can also lead to improved quantum-inspired algorithms for approximating matrix functions, e.g., permanent and hafnian.Moreover, it is still open to finding an efficient algorithm for the multiplicative-error approximation of a matrix function beyond the method using sampling from nonnegative quasiprobability distributions [15].Since a multiplicative-error approximation is significantly more powerful than an additive-error one and only a few examples have been known [24][25][26], such findings have intriguing applications in computational complexity.
In this work, we provide algorithms specialized for approximating the outcome probabilities of a linear optical circuit.First, we choose the s-parameterized quasiprobability distributions (s-PQDs), a generalization of P -, Q-, and Wigner distributions [29].An advantage of adopting s-PQDs over previous approaches is that we can always obtain a nonnegative representation for a Gaussian input state by choosing an appropriate parameter s.Consequently, for a Gaussian boson sampling (GBS) circuit, we can significantly reduce the negativity bound of the circuit, which depends only on the maximum peak of additive-error multiplicative-error (Eq.( 9)) i W (1/e) 2 ( * ) (Th. 1) #P-hard [27] Per(B) λmin = 0 : ϵ M i 4λ 2 max e(2λmax−λ i ) ( * ) (Th. 2) λmin = 0 : NP-hard [28] λmin > 0 : ϵ M i H B i (λi) ( * ) Eq. (A59) λmin > 0 : λmax λ min ≤ 2 [26] ( * ) the s-PQDs of measurement operators.Furthermore, the negativity bound is determined by the maximum possible s, called classicality of the input Gaussian state [30], which has a clear physical meaning: the more classicality in the input state, the lower negativity bound.
Our main technical contribution is to provide a way of manipulating the shape of s-PQDs using the symmetry of the circuit transformation in the phase space.This method can considerably lower the negativity bound, from exponential to at most polynomial in several cases, which renders efficient estimations of outcome probabilities within 1/poly additive-error.Strikingly, when the classicality of the input state is high enough, we introduce a fully polynomial-time randomized approximation scheme (FPRAS) by an efficient sampling from a log-concave function [26,31], which can efficiently (in BPP) approximate the corresponding outcome probability within a multiplicative-error.
Our results have several intriguing applications to problems in computational complexity.First, we give estimating algorithms with additive-errors for matrix functions represented by the outcome probabilities of a linear optical circuit, beating the best-known classical algorithms, e.g., the hafnian of a complex symmetric matrix and the permanent of an HPSD matrix.Second, we provide efficient multiplicative-error algorithms for those matrix functions with certain structured matrices, which are not considered in the literature to the best of our knowledge.Last but not least, we present sufficient conditions on the classical simulability of GBS, by applying our estimation scheme to any (marginal) outcome probability of a GBS circuit with 1/poly additive-error under a poly-sparsity condition.

II. SUMMARY OF RESULTS
We place a background material in Section III to provide preliminary information of quasiprobability and estimation schemes.The main results are following: • Our key technique for manipulating quasiprobability distribution in the phase space (Section IV A) • Scheme for additive/multiplicative-errors approximations of outcome probabilities of a linear optical circuit depending on the parameter regime (Section IV B) • Providing additive-error approximation algorithms for various matrix functions beating best-known classical algorithms (Section IV C 1), such as the hafnian of a complex symmetric matrix (Theorem 1) and the permanent of HPSD matrix (Theorem 2) • Providing multiplicative-error approximation algorithms for various matrix functions (Section IV C 2), such as the hafnian of a structured matrix (Theorem 3) • Sufficient conditions for the efficient simulability of lossy Gaussian boson sampling under a polysparsity assumption (Section IV D and Theorem 4) We summarize additive-and multiplicative-errors quantum-inspired algorithms for various matrix functions including Torontonian in Table I.

A. Quasiprobability representation of outcome probability
Let us consider the Born rule probabilities of a quantum optical circuit using quasiprobability distributions in the phase space.For an M -mode input state ρ in , a quantum channel E, and a measurement Π ν , the probability for a measurement outcome ν := (ν 1 , . . ., ν M ) can be written as [16] where Πν (β), and T (t,s) E (β|α) are s-PQDs of the input state, measurement, and the transition function of the circuit channel, respectively.Here, α is a quadrature variable in the 2M -dimensionnal phase space and β is a transformed quadrature variable by the transition function T (t,s) E (β|α).Specifically, W (s) ρ (α) is the s-PQD for a Hermitian operator ρ defined by where D(α ′ ) = e α ′ â † −âα ′ † is the M -mode displacement operator.Note that s = −1, 0, 1 of s-PQDs correspond to the Q-, Wigner, and P-distribution, respectively.Also, (−s)-PQDs of the measurement operators satisfy the normalization condition such that π M ν W (−s) Πν (β) = 1 for any β.The transition function of a quantum channel E is defined by [16] where for a coherent state |γ⟩.In this work, we are concerned with a linear optical circuit represented by a unitary matrix U with a product input state ρ in = ⊗ M i=1 ρ i and a product measurement Π ν = ⊗ M j=1 |ν j ⟩ ⟨ν j |.In that case, if we choose t = s, then the transition function becomes a simple form such as T E (β|α) = δ(β−U α).Consequently, Eq. ( 1) can be written in the simpler form where U ji α j by the transition function.Suppose we have a product Gaussian input state with the covariance matrix V i and a product photon number measurement Π mj .The s-PQD of a single-mode Gaussian state with zero-displacement having covariance matrix V is given by [32] Note that for any physical covariance matrix V , there exists a critical value s max such that W (s) V (α) is a proper Gaussian distribution for s < s max and has a delta function singularity for s = s max [18].We call s max 'classicality' of the Gaussian input state [30].Meanwhile, the (−s)-PQD for the photon number measurement operator |m⟩ ⟨m| can be represented [9,33] as where L m (x) is the m th Laguerre polynomial and s > −1.

Additive-error approximation
The first method is to estimate the Born rule probability within an additive-error.From the result in Ref. [5], one can estimate the outcome probability p(ν) within error ϵ with success probability 1 − δ for a given number of samples N = Here, M → is the (forward) negativity bound of the circuit defined as where ρi (α i ) is the negativity of the state ρ i in the phase space.When the negativity bound M → increases at most polynomially in the number of modes M , we can efficiently estimate the corresponding outcome probability with additive-error ϵ within running time T = poly(M, 1/ϵ, log δ −1 ).
Although this method is generally applicable for circuits having negativity, we usually encounter exponential negativity bound, i.e., M → = c M with c > 1, due to the input state negativity or a high maximum peak of the quasiprobability distribution of the measurement part.To circumvent this problem, finding a good quasiprobability representation for a given circuit with a small negativity is critical [34].For instance, let us consider a GBS circuit with the input of a pure squeezed vacuum state and photon number measurement.By choosing s-PQDs, the negativity of the input state is removed when s ≤ s max .However, the negativity bound is still exponential in the number of modes because of the high peaks of s-PQDs for the measurement operators except for a small squeezing parameter.

Multiplicative-error approximation
Remarkably, for special cases, we can reach a much stronger approximation, namely FPRAS.An FPRAS can estimate the target function p(ν) within a multiplicativeerror ϵ, which means for any 0 < ϵ < 1, 0 < δ < 1, the algorithm outputs µ such that with the running time poly(M, 1/ϵ, log δ −1 ).A sufficient condition for FPRAS is that the target function can be written as an integral of a log-concave function f (t), i.e., log [26,31].Let us consider a linear optical circuit as in Eq. ( 5).Since a multivariate Gaussian distribution is log-concave, s-PQD of a Gaussian input state fulfills the condition of log-concavity.If we choose a Gaussian measurement, the integrand satisfies the log-concavity but it is a trivial case.In this work, we develop a technique for making the quasiprobability of a non-Gaussian measurement logconcave, by manipulating quasiprobability in the phase space.

A. Manipulating quasiprobability in the phase space
In the previous section, we introduced two approximation schemes for the Born rule probability of an optical circuit.Those methods themselves do not seem to give powerful results for physically relevant situations.To obtain more interesting results, we propose a method manipulating the shapes of quasiprobability distributions using the symmetry of circuit transformation in the phase space.
In a general situation, one can rewrite the Born rule probability Eq. (1) as where W Πν (β)g(β) with appropriate functions h(α) and g(β).For the second equality, the transition function should satisfy a condition resulting by a symmetry in the phase space, which is given by Here, the auxiliary functions h(α) and g(β) have to be chosen so that the modified functions W ′ (t) ρin (α), T ′ (t,s) U (β|α), and W ′ (−s) Πν (β) are well-behaved in the phase space.An important point is that the modified functions do not need to be proper quasiprobability distributions of physical operators because we only need to exploit their shapes in the phase space.Consequently, these additional degrees of freedom allow us to further optimize the shapes of s-PQDs for our purposes.

B. Improved approximation for outcome probability of linear optical circuit
Let us first focus on improving the approximation scheme with additive-error using our method.Since P i (α i , V i , s, γ)'s are nonnegative distributions for a Gaussian input state, the modified negativity bound is given by M ′ → = M j=1 max βj f j (β j , Π mj , s, γ) .The advantage of our method is manifest especially when M → grows exponentially in the number of modes whereas M ′ → grows at most polynomially in the number of modes (see Fig. 1 (b)).Let us present a simple example for which our method works well.Consider an M -mode identical pure squeezed vacuum states input with squeezing parameter r > 0 and all single-photon outcomes m = (1, . . ., 1).Since the negativity of the photon number measurement operator is monotonically decreasing with growing s, we choose s = s max = e −2r [30].Then the negativity bound M → is exponential in M when the squeezing is high, i.e., r > (β j ) > 1.However, we can shift the Gaussian factor by choosing γ as (Appendix A 1) where W (x) is the Lambert W function.Note that γ * < 0, which means an inverse Gaussian function acts on the s-PQD of the measurement operator.From the fact that max βj |f j (β j , Π 1 , s max , γ * )| < 1, the modified negativity bound M ′ → is exponentially small in M for any squeezing r, which renders an efficient approximation with additiveerror ϵ within running time T = poly(M, 1/ϵ, log δ −1 ).Furthermore, our method provides an intriguing result on the approximation scheme with multiplicativeerror (see Fig. 1 (c)).For instance, we consider an M -mode identical squeezed thermal state (r, n th ) having high enough classicality s max > 1 and all single-photon outcomes m = (1, . . ., 1).In this case, an appropriate γ > 0 can make s-PQD of the measurement operator a log-concave function by adding Gaussian smoothing.

C. Application I: quantum-inspired algorithms for matrix functions
Permanent and hafnian are important matrix functions in computational complexity.Although computing these matrix functions is generally hard [21,[35][36][37], there are still efficient methods for matrices that have specific structures or restrictions [24,25,38,39].Developing algorithms for estimating matrix functions of particular classes of matrices is a highly nontrivial problem and might enable us to understand the hardness of the problem better.
One approach is using quasiprobability representations of matrix functions.In general, there can be several ways to match the matrix functions with the outcome probability of quantum circuits (Fig. 1 (a)).For example, an outcome probability of a linear optical circuit with photon number measurements and a Gaussian input state with zero-displacement having covariance matrix V can be written using a hafnian as [20] where and A S is a submatrix of A with repeated rows and columns depending on the detected photons.Meanwhile, if the measurements are threshold detectors, i.e., Π off = |0⟩ ⟨0| and Π on = 1 − |0⟩ ⟨0|, the corresponding probability is given in terms of Torontonian as [40] p where and m ′ is an M -element binary vector representing on/off measurement outcomes.Therefore, our approximating methods for an outcome probability p(m) are closely related to estimating the above matrix functions with additive or multiplicativeerrors.

Estimating algorithms with additive-errors
Suppose we have a Gaussian input state as a squeezed thermal state.When the thermal part is absent, corresponding to the standard GBS circuit, we can obtain an algorithm for estimating the absolute square of the hafnian of a complex symmetric matrix: with a success probability 1 − δ using the number of samples O(log δ −1 /ϵ 2 ) within the additive-error where λ max is the largest singular value of R.
Proof.First, we embed the hafnian of a complex symmetric matrix to an outcome probability of a GBS circuit after rescaling the matrix so that its singular values are on the interval [0, 1).More specifically, when the input is an M -mode product of pure squeezed states with the squeezing parameters {r i } M i=1 , the probability of a GBS circuit with all single-photon outcome m = (1, . . ., 1) is [20].Since any complex symmetric matrix can be decomposed as U DU T [41] with a unitary matrix U , our algorithm can be applied to general complex symmetric matrices.Meanwhile, this probability can be also written by using s-PQDs in the form of Eq. ( 14) such as where V sq,i is the covariance matrix of the squeezed state on the ith mode and N sq i 's are the normalization factors for P sq,i (α i , r i , s, γ)'s.Note that γ ∈ (− 2 s+1 , 2 e 2rmax −s ) for a given s and r max := max i r i .An appropriate choice of γ (see Appendix A 1) gives an upper bound on |f sq,j (β j , r j , γ, s)|: where λ j = tanh r j , and λ max := max j λ j .Then by Hoeffding inequality [42], where Given the success probability of the estimation 1 − δ and the number of samples O(log δ −1 /ϵ 2 ), we arrive at the result by substituting λ j = λ max to obtain an upper bound.
Note that the above algorithm gives the finest precision, to the best of our knowledge, i.e., ϵ(eλ max ) M in Ref. [43].However, the error is still larger than what we need for the hardness conjecture in GBS [44], so it does not lead to a contradiction.Next, if the input of the GBS circuit is a thermal state, we have an algorithm for the permanent of an HPSD matrix: Theorem 2. (Estimating permanent of HPSD matrices) For an M × M HPSD matrix B, one can approximate Per(B) with a success probability 1 − δ using the number of samples O(log δ −1 /ϵ 2 ) within the error where λ i are singular values of the matrix B and λ max is the largest one.
Proof.This corresponds to the case where the input is an M -mode product of thermal states with the average photon numbers {n i } M i=1 .The probability of all singlephoton outcomes matches the permanent of an HPSD matrix B such as Then we obtain the result by a similar procedure in the hafnian case.A detailed derivation is given in Appendix A 2.
The precision of our result outperforms the best-known one [23] because the latter is a special case of ours when s = 1 and γ = 0 (see Fig. 1 (b)).Specifically, when λ max ∈ (0, 1/2), our method's precision is better than the previous one (see Appendix A 2).Moreover, when the input is a squeezed thermal state, we obtain an algorithm for the hafnian of a structured matrix (Appendix A 1).Similarly, we provide algorithms for the Torontonian of some structured matrices within additive-error by substituting the photon number measurement by a threshold detector.The detailed results are in Appendix A 3.

Estimating algorithms with multiplicative-errors
Recall that we have an FPRAS when the estimate function is log-concave.Thus our goal is to make the estimate function log-concave by controlling the parameters s and γ in our scheme.For the permanent of an HPSD matrix, this is possible when λ max /λ min ≤ 2 (a proof in Appendix B 1), which reproduces the existing result of Ref. [26] in the case 1 ≤ λ i ≤ 2 after a normalization.Especially for an HPSD matrix with λ min = 0, estimating the permanent within a multiplicative-error is NPhard [28]; thus λ i > 0 is the crucial condition for the efficient approximation.We emphasize that our formulation comes from a physical setup, which is essentially different from the method in Ref. [26], where the technique is restricted to the properties of the permanent of a positive definite matrix.Thus our result can be readily extended to a more general situation other than the permanent of a positive definite matrix.When the input state is a product of squeezed thermal states {r i , n i } M i=1 , we have an FPRAS for the hafnian of a matrix having a specific form, such that 2 × 2 block matrix whose diagonal elements are symmetric matrices and off-diagonal elements are HPSD matrices.
where n = n i for all i and n, r i ≥ 0. Then Haf(A) can be approximated by FPRAS when the parameters satisfy a condition as We also consider an on/off measurement instead of the photon number measurement, which corresponds to the Torontonian (detailed analysis is given in Appendix B 3).
Moreover, we give nontrivial lower and upper bounds on the values of various matrix functions including the permanent and hafnian by adjusting the parameter s, which have independent interests (Appendix C) [45].

D. Application II: Sparse Gaussian boson sampling
Our algorithm has an interesting application to the simulation of GBS.Recall that if the input is a classical state, the corresponding GBS can be efficiently simulated [16].In our language, this is the case when s max ≥ 1.For a non-classical input state (s max < 1), however, a classically efficient simulation may not be possible from the hardness of GBS.Nevertheless, we can approximate any (marginal) outcome probability of the circuit using our algorithm under certain conditions.Then one can simulate the GBS when its output distribution is polysparse, in which a probability distribution with polynomially many of the most likely outcomes well approximates the true distribution [46,47].

Theorem 4. (Estimating outcome probabilities of GBS)
For a lossy GBS circuit with squeezing {r i } M i=1 and a transmissivity η, one can efficiently approximate any (marginal) outcome probability within 1/poly(M ) additive-error when ηe Proof.Let us first consider a GBS circuit with a product of lossy squeezed input states with a transmissivity η having the covariance matrix on ith mode as V η,i = We take s = s max = ηe −2rmax + 1 − η.If we examine the probability of all single-photon outcomes m = (1, . . ., 1), a condition for an efficient estimation is max βj |πW for n ≥ 2 and s ≥ 0. Lastly, we consider πW Π0 (β j ) for zero-photon detection and mj πW (−s) Πm j (β j ) for the marginalized probability on mode j, where the latter sum equals to one by the normalization condition.In both cases, the integrals for β j 's can be easily computed because β j components in W (s) Vη (α) constitute a multivariate normal distribution.Therefore, we can first perform the integrals on β j 's corresponding to zero-photon or marginalized ones and estimate the remaining terms.
We also consider GBS with threshold detectors [40], whose s-PQD for "click" is given by Since the range of πW 1] in the Wigner representation (s = 0), it allows any squeezing of the input state even for the lossless case (η = 1).Then we can estimate any (marginal) probability by the same argument of the photon-number measurement case.We give a detailed analysis in Appendix D.
In both cases, those algorithms of estimating (marginal) probability with inverse-polynomial additiveerror precision with the poly-sparsity condition imply classically efficient simulations of GBS.Since it is difficult to expect the sparsity condition for current experiments, our results do not lead to a direct simulation of them.However, it can be useful for application-targeted GBS experiments whose outcomes might satisfy the polysparsity condition.
Furthermore, we might consider an additional source of noise to allow s max > 1 by introducing thermal noise with the average photon number n th .In that case, the covariance matrix of ith mode input state is V η,n th ,i = .
By the log concavity condition, we can compute any (marginal) probability within a multiplicative-error if the following condition is satisfied: For instance, if η = 0.5 and r max = 1, then n th ≥ n * th ≃ 3.79, and the minimum value of n * th is 1 when η → 0. In Fig. 2, we depict our results of approximating probability and simulation for GBS circuits with threshold detectors via the classicality s max .Remarkably, there are two transition points of complexities via the classicality s max .i) s max = 1: This point indicates the transition from ϵ-simulation with the sparse condition to exact simulation.Whether ϵ-simulation without sparse condition can exist when s max < 1 is an interesting open question.ii) s max = s m : Note that s m is the point saturating the inequality Eq. (32), which indicates the transition of complexities from BPP NP to BPP for approximating the probability within a multiplicative error.

V. DISCUSSION
We propose a method for calculating the outcome probability of a linear optical circuit, by introducing modified functions with a lower negativity bound than FIG. 2. Regime of efficient classical algorithms for approximating outcome probabilities and the simulation of a GBS circuit with threshold detectors via the classicality smax.When the output distribution is poly-sparse, an approximate simulation is possible by estimating the probabilities within 1/poly additive-error (orange arrow) [46].On the other hand, a multiplicative-error approximation of probability in BPP NP is possible by using exact simulation with classical input state smax ≥ 1 (blue arrow) [15].
the quasi-probability distribution.This leads to various improved approximating algorithms for the outcome probabilities of the circuit.Furthermore, we suggest an FPRAS using our method modulating s-PQDs and the efficient sampling of log-concave functions with multiplicative-error. Our results provide a helpful tool for controlling the negativity of the circuit in the phase space and interesting quantum-inspired algorithms in computational complexity.
Although we focus on Gaussian input states and photon number or threshold measurements in a linear optical circuit in this work, our scheme can be also applied to other systems, for example, Clifford circuits [48] with a dimension of odd prime d.If one exploits the symmetry of the transition function in the phase space, a nontrivial approximation algorithm can be rendered for the corresponding matrix function.Since there can be several equivalent choices of circuits for the same matrix function, e.g., permanent of a unitary matrix [49,50], finding the optimal circuit is still an interesting open problem.After a proper circuit choice, there are other optimization problems, such as choosing quasiprobability representation and optimizing parameters for manipulating them in the phase space.Therefore, there might be more improvements for quantum-inspired algorithms approximating matrix functions.

DATA AVAILABILITY
The data supporting the results of this manuscript are given in the article and the appendix.Extra data are available upon reasonable request. as where V sq,i is the covariance matrix of a squeezed vacuum state in mode i and γ ∈ [0, 1) is the (normalized) parameter shifting the Gaussian factor such that γ → 1 (γ = 0) means maximum (no) shifting.To obtain a bound on |f sq,j (β j , r j , γ, s)|, let us set s = s max = e −2rmax and r j = tanh −1 λ j .The extreme points of f sq,j are 0, ±β * with Note that for s max < 1, f sq,j (0, λ j , γ, s max ) < 0 and f sq,j (β * , λ j , γ, s max ) > 0. Since the extreme values are changing monotonically as γ, we choose γ satisfying the condition −f sq,j (0, λ j , γ, s max ) = f sq,j (β * , λ j , γ, s max ), which is given by where W (x) is the Lambert W function, and W (1/e) 1−W (1/e) ≃ 0.386.Then an upper bound is obtained by substituting γ * into f sq,j (β * , λ j , γ, s max ), such that min Although this bound is valid only for a certain range of λ max , we can also find the same bound out of the range by shifting the Gaussian factor in the reverse direction.Specifically, where γ ′ ∈ [0, 1) is the parameter shifting the exponential term in the reverse direction such that γ ′ → 1 (γ ′ = 0) means maximum (no) shifting.Similarly, the extreme points are 0, ±β ′ * with (A17) Then the γ ′ satisfying the condition −f ′ sq,j (0, λ j , γ ′ , s max ) = f sq,j (β ′ * , λ j , γ ′ , s max ) is given by Finally we obtain an upper bound such that which covers the remaining range of λ max .Note that Consequently, by Eq. (A3), for a given complex M × M matrix R, and for the number of samples N = O(log δ −1 /ϵ 2 ), we can estimate |Haf(R)| 2 with a success probability 1 − δ within an additive-error as where λ ′ i = λi aλmax .We find upper and lower bound by setting λ i = λ max and λ i = 0, respectively, such as (A22) Next, we consider a squeezed thermal input state {r i , n} M i=1 to allow s max ≥ 1.For simplicity, the average thermal photon number n is fixed.Then the probability of all single-photon outcome is written by hafnian as where V st Q = V st + I 2M /2 with V st the covariance matrix of a squeezed thermal state and A = R B B T R * with a symmetric matrix R and an HPSD matrix B decomposed as Eqs.( 26) and (27).Let us define a ± (r i , n) = (2n+1)e ±2ri , a min = (2n + 1)e −2rmax , and a max = (2n + 1)e 2rmax .Then the probability p st is written by s-PQDs as where we use the reverse shifting of the Gaussian factor.To obtain an upper bound on |f ′ st,j (β j , r j , n, γ ′ , s)|, let us set s = s max = a min .The extreme points are 0, ±β ′ * with Finally, for a matrix A satisfying Eqs. ( 26) and ( 27), we can estimate Haf(A) with a success probability 1 − δ using number of samples N = O(log δ −1 /ϵ 2 ) within the additive-error given by 2. Permanent of Hermitian positive semidefinite matrices (Proof of Theorem 2) When a thermal state input with average photon numbers {n i } M i=1 goes through a linear optical circuit instead of a squeezed vacuum state, the probability of all single-photon outcomes corresponds to the permanent of HPSD matrices [23].In Ref. [23], an algorithm for estimating the permanent of an HPSD matrix is proposed.Here, we improve the precision of the estimation using s-PQDs and shifting Gaussian factors.The probability of all singlephoton outcomes is connected to the permanent of an HPSD matrix such that [23] where Thus if we have an M × M HPSD matrix B, firstly we rescale the matrix with the largest eigenvalue λ max as B ′ = B/(aλ max ) with a > 1, and find its unitary diagonalization such as B ′ = U DU † so that we can find a GBS circuit U with thermal input state {n i } M i=1 , whose probability matches Per(B ′ ).Then, If an estimator of Per(B) lies in the interval [−C M , C M ], by Hoeffding's inequality [42], where µ is the sample mean of p th .Thus for the number of samples N = O(log δ −1 /ϵ 2 ), we can estimate Per(B) with a success probability 1 − δ within an additive-error ϵ(aλ max CZ ′1/M ) M .Now we use the same method as in the hafnian case by introducing the shifting parameter γ ∈ [0, 1), such as f th,j (β j , n j , γ, s).
Next, we compare our algorithm's precision with Gurvits' randomized algorithm for the permanent of a general complex matrix A giving the additive-error as ϵ||A|| M = ϵλ M max with samples N = O(M 2 /ϵ 2 ).Thus from Eq. (A55), the necessary and sufficient condition for beating Gurvits' precision is written as To get lower and upper bounds, we put λ i = 0 and λ i = λ max for all i, respectively, such as Also for 0 < λ min < 1/2, the additive-error for Per(B) is given by where Similarly, the necessary and sufficient condition for beating Gurvits' is given by

Torontonian
The Torontonian of a 2M × 2M complex matrix A ′ is defined as [40] Tor where P ([M ]) is the power set of [M ] := {1, 2, ..., M } and the matrix A ′ has block structure such that where B is HPSD and R is symmetric.Let us first consider a special case B = 0, which corresponds to a GBS with pure squeezed input state {r i } M i=1 .The probability of all threshold detectors "click" in a GBS circuit is related to the Torontonian as where , and D = ⊕ M i tanh r i .Meanwhile, the probability p on|sq can be written in terms of s-PQDs as where P sq,i (α i , r i , γ, s) is given in Eq. (A9) and f on|sq,j (β j , r j , γ, s) is defined as We set s = s max = e −2rmax and r j = tanh −1 λ j .The extreme points are at 0, ±β * with After we choose γ * = 1 2 (1 − λ max ), the corresponding upper bound is given by min Therefore, with success probability 1 − δ, we can estimate the Torontonian of a matrix 0 R * R 0 using the number of samples O(log δ −1 /ϵ 2 ) within an additive-error Next, we consider a special case where R = 0 corresponds to a thermal input state {n i } M i=1 .In that case, where Meanwhile, the probability p on|th can also be written as where P th,i (α i , n i , γ, s) is given in Eq. (A38) and f on|th,j (β j , n j , γ, s) is defined as We set s = s max = 2n min + 1 and n j = λj 1−λj .The extreme points are at 0, ±β * with Choosing γ * = 1 2 (1 − λ max ), an upper bound on |f on|th,j (β j , λ j , γ, s)| is given by min Therefore, we can estimate Torontonian of a matrix B T 0 0 B with success probability 1 − δ using number of samples O(log δ −1 /ϵ 2 ) within an additive-error Finally, we consider a matrix A ′ defined by Eq. (A63), satisfying Eqs. ( 26) and (27).Then the Torontonian of matrix A ′ is related with a squeezed thermal input state {r i , n} M i=1 such that Meanwhile, p on|st is also can be written as where a ± (r i , n) = (2n + 1)e ±2ri , a max = (2n + 1)e 2rmax , and a min = (2n + 1)e −2rmax .We set s = s max = a min and note the extreme points are 0, ±β * with Then by straightforward calculation, From Eq. (A36), Therefore, the condition for log-concavity of f th,j (β j , n i , γ, s) is written as After putting γ = 1 and s = s max = 2n min + 1, Substituting as desired.

Hafnian (Proof of Theorem 3)
When the input is a squeezed thermal state {r i , n} M i=1 , we can find out the condition for multiplicative estimation of the hafnian of a particular matrix.Note that the estimated function in the case of all single-photon outcomes is given by f st,j (β j , r j , n, γ, s) ∝ where a max = (2n + 1)e 2rmax .Then from Lemma. 1, the condition for log-concavity of f st,j (β j , r j , n, γ, s) is where we put γ = 1 and s = s max = (2n+1)e −2rmax .Thus for a matrix A = R B B T R * satisfying conditions Eqs. ( 26) and (27), Haf(A) can be estimated within a multiplicative-error efficiently when Eq. (B8) holds.

Torontonian
We can apply the same method to the Torontonian.Here, we use the following lemma: Lemma 2. Let a, b, c ≥ 0 and q : R M → R + is a positive semidefinite quadratic form.Then the function (a − be −bq(x) )e −cq(x) is log-concave when a ≥ b 2 +2bc c .
Proof.By the same argument in Lemma. 1, what we need to check is g ′′ (τ ) ≤ 0 for all τ ∈ R, where By a straightforward calculation, First we consider thermal input state {n i } M i=1 .From Eq. (A73), By Lemma.2, the function (a − be −b|β| 2 )e −c|β| 2 is log-concave when a ≥ b 2 +2bc c .Consequently, the condition for log-concavity of f on|th,j (β j , n j , γ, s) is written as When γ = 1 and s = s max = 2n min + 1, this condition yields where λ j = nj 1+nj .Thus for an HPSD matrix B, Tor B T 0 0 B can be estimated within a multiplicative-error when eigenvalues {λ i } of C satisfy the condition Eq. (B16).Note that this condition is more stringent than the permanent case in which there is no restriction for λ max when λ min ≥ 1 2 .Next, when the input is a squeezed thermal state {r i , n} M i=1 , f on|st,j (β j , r j , n, γ, s) is given by where a max = (2n + 1)e 2rmax .If we set γ = 1 and s = s max = (2n + 1)e −2rmax , the condition for log-concavity of f on|st,j (β j , r j , n, γ, s) is given by using Lemma.
Tor(B ′ ) TABLE II.Upper and lower bounds on various matrix functions.

Permanent of Hermitian positive semidefinite matrices
From Eq. (A33), the permanent of an HPSD matrix is connected to a GBS circuit with a thermal state input.The probability of all single-photon measurements is written as where the inequality is valid when s ≥ 1.Note that we take s = 2n min + 1 and n i = λi 1−λi in the last equality.
Thus the permanent of an HPSD matrix B are bounded from above as For the lower bound, our method gives the same result in Ref. [23].

Hafnian
Let us set a i,± (r i , n) = (2n + 1)e ±2ri , a min = (2n + 1)e −2rmax , and a max = (2n + 1)e 2rmax .When the input state is a squeezed thermal state {r i , n} M i=1 with a fixed n, the probability of all single-photon detection is where the inequality holds for s ≥ 1, and we put s = 1 for the last equality.Here, assume a min ≥ 1.Then for a matrix A = R B B T R * satisfying Eqs. ( 26) and ( 27), a lower bound of the Hafnian is written as Similarly, an upper bound can be obtained as where we put s = 1 for the last inequality.Consequently, an upper bound is obtained as

Torontonian
First, we consider thermal state input {n i } M i=1 and the probability of all "click" outcomes is written as where the inequality is valid when s ≥ 1, and we take s = 1 for the last equality.Thus for an M × M HPSD matrix B, the Torontonian of Similarly, we also obtain an upper bound as where we take s = 1 for the last equality.Consequently, Next, when the input state is a squeezed thermal state {r i , n} M i=1 .Then the probability of all "click" detection is where a i,± (r i , n) = (2n + 1)e ±2ri , a min = (2n + 1)e −2rmax , a max = (2n + 1)e 2rmax , and we put s = 1 for the last inequality.Assume a min ≥ 1.Then for a matrix A ′ = B T R * R B satisfying Eqs. ( 26) and ( 27), a lower bound of the Torontonian is given by Tor An upper bound can be obtained by a similar method, such as where we put s = 1 for the last inequality.Consequently, an upper bound of the Tor(A ′ ) is obtained as Appendix D: Simulability of Gaussian boson sampling (Proof of Theorem 4) Our approximation algorithm for outcome probabilities of a linear optical circuit have applications not only for the matrix functions, but also for Gaussian boson sampling, which is crucial for the demonstration of quantum supremacy [20].From the results in Refs.[47], we have three level of a hierarchy of notions of classical simulation as following: 1. Poly-box : Inverse-polynomial additive-error approximation of any outcome probabilities including any marginals.
2. ϵ-simulation : Approximate sampling simulation of probability distributions with ϵ-error in total variation distance.
3. Multiplicative precision estimation : multiplicative-error approximation of any outcome probabilities including any marginals.
A poly-box can be promoted to ϵ-simulation when the outcomes are poly-sparse, and multiplicative precision estimator implies ϵ-simulation [47].In our work, we can investigate this hierarchy in the GBS via a degree of classicality of the input state, s max .To do that, we consider a lossy GBS, in which the input state is product of lossy squeezed state having the covariance matrix on ith mode with V i = 1 2 ηe 2ri + 1 − η 0 0 ηe −2ri + 1 − η := 1 2 a i+ (η, r i ) 0 0 a i− (η, r i ) .
From Eq. ( 6), −1 < s ≤ s max ≤ 1 for a lossy squeezed state, and s max = a − (η, r max ).Note that the s max goes to 0 as the maximum squeezing parameter r max → ∞, which is consistent with the fact that a general Gaussian state can be for given η, r max , and γ = 0, s = s max for the inequality.Then from the condition max βj |f j (β j , η, r max , 1, 0, s max )| ≤ 1, s max satisfies s max ≥ √ 5 − 2 ≃ 0.236.This corresponds to r max ≤ 1 2 log(2 + √ 5) ≃ 0.722 for an ideal GBS (η = 1).However, if we allow the photon loss, any squeezed input state is possible when η ≤ 3 − √ 5 ≃ 0.764, which is much higher transimissivity than those used in current experiments [12,13].Next, we need to check whether this condition is valid for any other outcomes.From the behavior of f j (β, η, r, m, 0, s), we can find out that for m ≥ 2 and s ≥ 0. Finally, we consider n = 0 for zero-photon detection and f j = 1 for the marginalized probability owing to the normalization of measurement operators.In both cases, the integrals for β j 's can be easily computed because f j (β j ) and β j components in P (α) are just Gaussian distributions.Therefore, we can always perform the integrals including β j 's corresponding to zero-photon or marginalized one, and estimate remaining terms.Furthermore, we examine the case of threshold detectors instead of number resolving measurements [40].The corresponding f on,j for a 'click' event is written as f on,j (β j , η, r j , γ, s) = a + (η, r max ) − s a + (η, r max ) − s − γ(a j+ (η, r j ) − s) a + (η, r max ) − s a + (η, r max ) − s − γ(a j− (η, r j ) − s) Then for γ = 0 and s = 0, the range of f on,j (β j , η, r j , 0, 0) = 1 − 2e −2|βj | 2 is on [−1, 1], thus the poly-box condition is satisfied for all input squeezing and loss parameter.Now we investigate whether an efficient estimation of GBS probability within multiplicative-error is possible.To do that, we consider Gaussian states which can have s max > 1, where the covariance matrix of ith mode state is given by a i+ (η, n th , r i ) 0 0 a i− (η, n th , r i ) . (D13) These are squeezed thermal states, in which pure squeezed states undergo a thermal noise with average photon number n th instead of the vacuum loss.In this case −1 < s ≤ a − (η, n th , r max ) for given η, n th , and r max .Then we need to check the log-concavity of f on,j such that f on,j (β j , η, n th , r i , γ, s) For instance, if η = 0.5 and r max = 1, then n th ≥ n * th ≃ 3.79 for the multiplicative-error estimation of the probability, and the minimum value of s max is 3 when η → 0.

FIG. 1 .
FIG. 1. Schematic diagram of quantum-inspired classical algorithms for approximating matrix functions.For a given matrix function, (a) find an embedding of the matrix function onto an outcome probability of a quantum circuit (ρ, U, Π) and choose a quasiprobability representation of the probability.(b) We depict an example of a linear optical circuit, for the approximation scheme with additive-error.Using s-PQDs for the linear optical circuit, one can significantly reduce the negativity bound by appropriately choosing γ < 0. (c) Approximation scheme with multiplicative-error.When the classicality of the input state is large, one can make the s-PQDs of the measurement operator a log-concave function by choosing a suitable γ ′ > 0.

Theorem 3 .
(FPRAS for hafnian) Suppose we have a block matrix A = R B B T R * with an M × M complex symmetric matrix R and an M × M HPSD matrix B, which have decompositions by a unitary matrix U as U DU T and U D ′ U † , respectively, with