Connection between BosonSampling with quantum and classical input states

BosonSampling is a problem of sampling events according to the transition probabilities of indistinguishable photons in a linear optical network. Computational hardness of BosonSampling depends on photon-number statistics of the input light. BosonSampling with multi-photon Fock states at the input is believed to be classically intractable but there exists an efficient classical algorithm for classical input states. In this paper, we present a mathematical connection between BosonSampling with quantum and classical light inputs. Specifically, we show that the generating function of a transition probability for Fock-state BosonSampling (FBS) can be expressed as a transition probability of thermal-light inputs. The closed-form expression of a thermal-light transition probability allows all possible transition probabilities of FBS to be obtained by calculating a single matrix permanent. Moreover, the transition probability of FBS is shown to be expressed as an integral involving a Gaussian function multiplied by a Laguerre polynomial, resulting in a fast oscillating integrand. Our work sheds new light on computational hardness of FBS by identifying the mathematical connection between BosonSampling with quantum and classical light. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
When indistinguishable photons enter a linear optical network, multimode interference draws various output-photon distributions. If the input photons are prepared in a multi-photon Fock state, each transition probability is calculated with the absolute square of a matrix permanent [1]. Even though the computational complexity of computing the absolute square of permanents is still not rigorously proven, it is strongly believed to be in the same complexity class as computing permanents, which is a #P-hard problem [2,3]. Thus, generating samples according to the transition probabilities, namely Fock-state BosonSampling (FBS), is considered to be a classically intractable problem [4][5][6][7]. To demonstrate that a quantum device can outperform a classical computer experimentally, large-scale FBS devices are under development [8][9][10][11].
Besides Fock-state inputs, BosonSampling with various input sources has been studied in terms of the computational complexity [12][13][14][15][16][17][18][19][20][21]. For instance, the transition probability from thermal light to a specific output photon distribution can be estimated from the permanent of a Hermitian positive semidefinite matrix (HPSM) [22][23][24][25]. Because there is an efficient classical algorithm for thermal-light BosonSampling (TBS), by using Stockmeyer's approximate counting algorithm, the computational complexity of approximating the transition probability, i.e. the permanent of an HPSM, is classified as BPP NP [24,26]. In addition, inspired by TBS, a classical sampling algorithm for approximating the permanent of an HPSM within an additive error has been proposed [25].
In this paper, we present a mathematical connection between BosonSampling with quantum and classical light inputs. Specifically, we show, in Sec. 2, that the generating function of transition probabilities for FBS can be interpreted as a transition probability of TBS. Because of their relations with matrix permanents, the result also draws a connection between the permanents of submatrices formed from a unitary matrix and an HPSM. The closed-form expression of a transition probability of TBS allows all possible transition probabilities of FBS to be obtained by calculating a single matrix permanent. However, in Sec. 3, we show that the method turns out to be exponentially slower than the brute-force permanent calculations for all possible transitions. Furthermore, by exploiting the connections, we express the transition probabilities of FBS in a definite integral involving a Gaussian function multiplied by a Laguerre polynomial in Sec. 4. Considering the fast oscillating integrand, we discuss the computational hardness of FBS and the absolute square of a matrix permanent. Finally, we discuss the potential application of our connections for estimating output mode correlations of FBS, in Sec. 5, which is expected to be used to certify whether a given BosonSampling device is properly working [27,28].

Generating function for Fock-state BosonSampling
We consider a lossless M-mode linear optical interferometer for FBS, described by a unitary operatorÛ. The transition probability from a multi-photon Fock input state |n = M k=1 |n k to the output state |m = M k=1 |m k is given by P f (n, m) = | m|Û|n | 2 . Here, n k and m k represent the number of photons in the k-th input and output modes, respectively, and the total number of photons is conserved as N = M k=1 n k = M k=1 m k because of the lossless condition. Through the use of the unitary relation for the bosonic creation operator ofÛâ † jÛ † = m k=1 U jkâ † k , the transition probability can be rewritten in the absolute square of a matrix permanent as follows: where [U] n,m is the N × N submatrix formed by taking the entry U jk 'm j ' times for row and 'n k ' times for column from the unitary matrix U [1].
Here, we construct a generating function of all possible transition probabilities for the input Fock state |n . Through the introduction of an indeterminate z = (z 1 , . . . , z M ), the generating function can be built as [29][30][31] where every possible M-mode output configuration is assigned to the summation index. Note that if the total photon number of an output state differs from N, the transition probability is zero. A specific transition probability P f (n, m) can be extracted by taking the m k -th partial derivative of the generating function with respect to each z k and making z = 0, i.e., When P f (n, m) is substituted with its bra-ket notation | m|Û|n | 2 , the generating function is recast into The term inside the parentheses is proportional to a multi-mode thermal-light stateρ th = M k=1ρ th k , withρ th k = (1 − z k ) ∞ m k =0 z m k k |m k m k | following Bose-Einstein statistics represented by a mean photon number of z k /(1 − z k ). Consequently, the generating function can also be expressed as Intriguingly, the above equation shows that the generating function is associated with TBS [22][23][24][25] because the numerator n|Û †ρthÛ |n corresponds to the transition probability P th r (ρ th , n) from thermal lightρ th to |n via the reversed linear optical interferometer ofÛ † . Here, the subscript r is introduced to emphasize that the linear optical interferometer is reversed.
From Eqs. (2) and (5), we can find a mathematical relationship between FBS and TBS: where the subscript r can be dropped by replacing the linear optical interferometerÛ † witĥ U. Herein, we use the fact that time reversal gives the same transition probability, i.e., P f (n, m) = P f r (m, n). Figure 1 shows a physical interpretation of how the transition probabilities of the two different input sources are related. The multi-mode thermal light state is a statistical mixture of photon-number states,ρ th = ∞ m=0 P BE (m)|m m|, following the Bose-Einstein The transition of thermal light can therefore be considered as the incoherent sum of the transitions of multi-photon Fock states |m satisfying the photon number conservation, with a weight of P BE (m), as described in Eq. (6).
That is, the transition of thermal light can be considered as the incoherent sum of infinite transitions of multi-photon Fock states |m , and the thermal-light transition probability P th (ρ th , n) can be expressed as the sum of the Fock state transition probability P f (m, n) with a weight of P BE (m). The connection can be summarized as The relationship in Eq. (6) can be easily extended for the permanents of two types of matrices because each of the transition probabilities is expressed with the permanent of different matrix types. Similar to Eq. (1), the transition probability of TBS can be written as a matrix permanent as follows: where Z = diag(z) and UZU † is a Hermitian positive semidefinite matrix (HPSM), given z k ≥ 0 [24,25]. Because of Eq. (6), the permanents of submatrices formed from the unitary matrix and the HPSM are connected as The set N for the summation index consists of m satisfying the photon number conservation M k=1 m k = M k=1 n k .

One-shot calculation of all transition probabilities for Fock-state BosonSampling
The generating function in Eq.
(2) encodes all P f (n, m) as the coefficients of a power series M k=1 z m k k . Interestingly, through the combination of Eqs. (5) and (7), the generating function can be expressed in a single matrix permanent, Therefore, every transition probability of FBS can be extracted from the coefficients of M k=1 z m k k in the polynomial expansion of Perm([U † ZU] n,n ). At a glance, it seems that there might be an exponential improvement in the computational cost because it only requires a single matrix permanent, Perm([U † ZU] n,n ), for evaluating all P f (n, m). To compare the one-shot calculation with a brute-force method, we investigate the computational costs of both methods. Firstly, we analyze the cost of Perm(A) based on Ryser's formula, known for being the most efficient classical algorithm to date for calculating a permanent of a matrix [32,33], Here, a ij is an element of an N × N matrix A and the index s of the outer summation is assigned by all possible subsets of {1, . . . , N}. |s| refers to the cardinality or number of elements of the subset s. The calculation of Ryser's formula can be optimized by using a Gray code, which picks the subsets so that only a single element is changed at a time [34]. If a ij is a number, the cost of the inner summation of a ij over j ∈ s can be lowered  [35]. We note here that the symmetries of the submatrix structures due to the multiplicities in n and m are not considered in the evaluation of the permanents [36][37][38], so the computational cost can be further reduced.
In case of the one-shot method, the elements of [U † ZU] n,n are multivariate irreducible polynomials of z k where k runs from 1 to M. Therefore, we need to analyze the cost of computing the permanent newly. The cost of the inner summation over j ∈ s in Eq. (10) cannot be reduced to less than O(M) because each coefficient of z k has to be separately calculated. When the cost of the multiplication over i and the outer summation over s is considered without polynomial expansions, the total cost for calculating Perm([U † ZU] n,n ) is given as O(MN2 N ). The cost is exponentially more efficient than that of the brute-force method O(e N (M/N) N N 1/2 2 N ) because M should be larger than N for FBS [3]. However, the polynomial expansion of Perm([U † ZU] n,n ) is necessary to extract the transition probability P f (n, m) from the coefficient of M k=1 z m k k , as seen in Eq. (9). Unfortunately, the main computational cost arises from the polynomial expansion.

Definite integral form of transition probabilities for Fock-state BosonSampling
The computational hardness of TBS was investigated by expressing P th (ρ th , n) in a definite multi-dimensional integral with a non-negative integrand. Because the non-negative integrand can be considered as a probability density function, it allows an efficient sampling algorithm for TBS [24]. Similarly, we construct a definite integral form of P f (n, m) by exploiting the relation in Eq. (6) and a integral form of P th (ρ th , n), and discuss the computational hardness of FBS. Firstly, we recast the transition probability of thermal light in a definite integral form. For this, the multi-mode thermal light state is written as the Glauber-Sudarshan P-representation.
The representation shows that thermal-light state can be considered as a Gaussian mixture of multi-mode coherent states |α = ⊗ M k=1 |α k [39,40]. Similar to Fig. 1, the transition of thermal light can be decomposed into an infinite number of transitions of coherent states |α . Multi-mode interference of |α via a lossless linear optical interferometerÛ † results in the output coherent state of |β = ⊗ M j=1 | β j , where β j = M k=1 U * kj α k and M k=1 |α k | 2 = M j=1 | β j | 2 . The transition probability from |α to a specific multi-photon Fock state |n is therefore given as | n|β | 2 , i.e., By exploiting the relation in Eq. (11), the transition probability of thermal light to |n can be expressed as a Gaussian mixture of the transition probability P cs r (α, n) as follows: where the integration is over the whole complex plane C M for all α k . Note that the integrand in Eq. (13) is nonnegative and classically tractable, so that there is an efficient sampling algorithm for TBS [24].
To find a definite integral form for FBS, we substitute the thermal light transition probability of Eq. (6) with Eq. (13) and apply partial derivatives and set z = 0 on each side.
The m k -th partial derivatives with respect to z k in the second line are identified as where L m k (x) is a Laguerre polynomial, such that When Eqs. (12) and (15) are applied to Eq. (14), the probability is rewritten as To simplify the integrand, we rearrange the equation by supposing that every z k approaches to 0 with the same rate, i.e., z k = z for all k. The integral variable α k is then replaced with γ k = α k / √ z, and Eq. (17) becomes where η k = β k / √ z = M j=1 U * jk γ j . Because M k=1 z n k −m k = 1 for N = M k=1 m k = M k=1 n k , the transition probability of FBS becomes a definite multi-dimensional integral without z k , as follows: The integral can also be used for computing the absolute square of permanents for unitary submatrices by using the relation in Eq. (1). Note that this integral form of P f (n, m) is more generalized than the integral in Ref. [41], which was derived only for m k = 0 or 1. Unlike the transition probability of thermal light in Eq. (13), the integrand of P f (n, m) in Eq. (19) could be negative. Therefore, the integrand cannot be considered as a probability density function and the efficient sampling would be hard. Furthermore, the integrand is a highly oscillatory function of γ k because of the Laguerre polynomials, as seen in Eq. (16). The numerical integration of such oscillating function is known to be a difficult task even though some techniques have been developed for some specific classes [42] and the numerical integration in high dimension becomes intractable as the dimensionality of the integral space grows, which is the so-called "Curse of Dimensionality" [43]. The hardness of integrating the highly oscillatory function would be interpreted as the reason for the hardness of FBS. Note that similar highly oscillatory integrands occur in another quantum sampling problem [44].

Conclusion
We have determined that a transition probability of TBS can serve as a generating function of transition probabilities for FBS. The connection was explained via both mathematical description and physical interpretation. The closed-form expression of a transition probability of TBS allowed the calculation of all possible transition probabilities of FBS with a single permanent calculation, although the one-shot calculation consumes more time than the brute-force method does. Moreover, the connection was applied to derive a definite integral for the Fock-state transition probability. With the result, we discussed the computational complexity of FBS by associating it with the hardness of integration. Our work sheds new light on computational hardness of FBS by identifying the mathematical connection between BosonSampling with quantum and classical light.
Our connection could become very versatile by assigning certain numerical values to z k and differentiating Eq. (9) with respect to z k . For instance, the probability of excluding the photon-click events at an output mode i can be obtained from Perm([U † ZU] n,n )/ M k=1 n k ! with z i = 0 for the i-mode and z k = 1 for the other modes. In addition, by taking the derivatives of Perm([U † ZU] n,n )/ M k=1 n k ! with respect to z i and z j and making z = 1, one can estimate the mode correlation between i and j modes, i.e., m imj = ∞ m=0 m i m j P f (n, m) with the bosonic number operatorm i . Even with the identical linear interferometer, the mode correlation differs across types of input sources and the indistinguishability of particles. Based on this feature, the analysis of the mode correlation can be utilized to benchmark BosonSampling whether a device generates output-photon distributions according to the transition probabilities that is difficult to classically simulate [27,28]. Note that the more z k are assigned by numbers, the lower the computational cost of Perm([U † ZU] n,n ).