1. 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–7]. To demonstrate that a quantum device can outperform a classical computer experimentally, large-scale FBS devices are under development [8–11].

Besides Fock-state inputs, BosonSampling with various input sources has been studied in terms of the computational complexity [12–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–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$^{\textrm {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].

2. Generating function for Fock-state BosonSampling

We consider a lossless $M$-mode linear optical interferometer for FBS, described by a unitary operator $\hat {U}$. The transition probability from a multi-photon Fock input state $\vert \mathbf {n}\rangle =\bigotimes _{k=1}^{M}\vert n_{k}\rangle$ to the output state $\vert \mathbf {m}\rangle =\bigotimes _{k=1}^{M}\vert m_{k}\rangle$ is given by $P^{f}(\mathbf {n},\mathbf {m})=|\langle \mathbf {m}|\hat {U}|\mathbf {n}\rangle |^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=\sum _{k=1}^{M}n_{k}=\sum _{k=1}^{M}m_{k}$ because of the lossless condition. Through the use of the unitary relation for the bosonic creation operator of $\hat {U}\hat a_j^{\dagger }\hat {U}^{\dagger }=\sum ^{m}_{k=1}U_{jk}\hat {a}_k^{\dagger }$, the transition probability can be rewritten in the absolute square of a matrix permanent as follows:

Here, we construct a generating function of all possible transition probabilities for the input Fock state $\vert \mathbf {n}\rangle$. Through the introduction of an indeterminate $\mathbf {z}=(z_{1},\ldots ,z_{M})$, the generating function can be built as [29–31]

From Eqs. (2) and (5), we can find a mathematical relationship between FBS and TBS:

 

Fig. 1. A pictorial description of the connection between thermal-light BosonSampling and Fock-state BosonSampling. The multi-mode thermal light state $\hat {\rho }_{th}$ is a statistical mixture of photon-number states $|\mathbf {m}\rangle$ following the Bose-Einstein statistics $P_{BE}(\mathbf {m})=\prod _{k=1}^{M}(1-z_k)z_k^{m_k}$. That is, the transition of thermal light can be considered as the incoherent sum of infinite transitions of multi-photon Fock states $|\mathbf {m}\rangle$, and the thermal-light transition probability $P^{th}(\hat {\rho }^{th},\mathbf {n})$ can be expressed as the sum of the Fock state transition probability $P^{f}(\mathbf {m},\mathbf {n})$ with a weight of $P_{BE}(\mathbf {m})$. The connection can be summarized as $P^{th}(\hat {\rho }^{th},\mathbf {n})=\sum _{\mathbf {m}=\mathbf {0}}^{\boldsymbol {\infty }}P_{BE}(\mathbf {m})P^{f}(\mathbf {m},\mathbf {n})$.

Download Full Size | PDF

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:

3. One-shot calculation of all transition probabilities for Fock-state BosonSampling

The generating function in Eq. (2) encodes all $P^{f}(\mathbf {n},\mathbf {m})$ as the coefficients of a power series $\prod ^M_{k=1}z_k^{m_k}$. Interestingly, through the combination of Eqs. (5) and (7), the generating function can be expressed in a single matrix permanent,

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],

According to Eq. (1), each transition probability can be calculated from the permanent of an $N\times N$ matrix $[U]_{\mathbf {n},\mathbf {m}}$ and the number of possible output states is given as a binomial coefficient, $\binom {M+N-1}{N}$, for a lossless $M$-mode linear optical interferometer and $N$ indistinguishable photons. As a result, the brute-force method, which calculates all the transition probabilities directly, has a computational cost of $\mathcal {O}(\binom {M+N-1}{N}N2^N)\approx \mathcal {O}(e^N(M/N)^N N^{1/2}2^N)$ [35]. We note here that the symmetries of the submatrix structures due to the multiplicities in $\mathbf {n}$ and $\mathbf {m}$ are not considered in the evaluation of the permanents [36–38], so the computational cost can be further reduced.

In case of the one-shot method, the elements of $[U^{\dagger }ZU]_{\mathbf {n},\mathbf {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\in s$ in Eq. (10) cannot be reduced to less than $\mathcal {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 $\textrm {Perm}([U^{\dagger }ZU]_{\mathbf {n},\mathbf {n}})$ is given as $\mathcal {O}(MN2^N)$. The cost is exponentially more efficient than that of the brute-force method $\mathcal {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 $\textrm {Perm}([U^{\dagger }ZU]_{\mathbf {n},\mathbf {n}})$ is necessary to extract the transition probability $P^{f}(\mathbf {n},\mathbf {m})$ from the coefficient of $\prod ^M_{k=1}z_k^{m_k}$, as seen in Eq. (9). Unfortunately, the main computational cost arises from the polynomial expansion. For the polynomial expansion of the product in Eq. (10), the multiplication of $N$ polynomials having $M$ irreducible terms costs $\mathcal {O}(M^N)$ for each $s$. Bacause there are $2^N$ of subsets $s$, the total cost of the multiplication becomes $\mathcal {O}(M^N 2^N)$. In addition, the outer summation for $\binom {M+N-1}{N}$ of irreducible terms resulting from the multiplication requires a computational cost of $\mathcal {O}(\binom {M+N-1}{N}2^N)$. Consequently, the total cost of computing $\textrm {Perm}([U^{\dagger }ZU]_{\mathbf {n},\mathbf {n}})$ to extract all configurations of the transition probabilities of the FBS is $\mathcal {O}(M^N2^N)$ when the most dominant cost is considered. According to a comparison of the cost of one-shot calculation of $\textrm {Perm}([U^{\dagger }ZU]_{\mathbf {n},\mathbf {n}})$ with that of brute-force method, the one-shot calculation requires an order of $\mathcal {O}((N/e)^N N^{-1/2})$ times the computational cost of the brute-force method.

4. Definite integral form of transition probabilities for Fock-state BosonSampling

The computational hardness of TBS was investigated by expressing $P^{th}(\hat {\rho }^{th},\mathbf {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}(\mathbf {n},\mathbf {m})$ by exploiting the relation in Eq. (6) and a integral form of $P^{th}(\hat {\rho }^{th},\mathbf {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.

Similar to Fig. 1, the transition of thermal light can be decomposed into an infinite number of transitions of coherent states $|\boldsymbol {\alpha }\rangle$. Multi-mode interference of $|\boldsymbol {\alpha }\rangle$ via a lossless linear optical interferometer $\hat {U}^{\dagger }$ results in the output coherent state of $|\boldsymbol {\beta }\rangle =\otimes ^M_{j=1}|\beta _j\rangle$, where $\beta _j=\sum ^M_{k=1}U_{kj}^*\alpha _k$ and $\sum ^M_{k=1}|\alpha _k|^2=\sum ^M_{j=1}|\beta _j|^2$. The transition probability from $|\boldsymbol {\alpha }\rangle$ to a specific multi-photon Fock state $|\mathbf {n}\rangle$ is therefore given as $|\langle \mathbf {n}|\boldsymbol {\beta }\rangle |^2$, i.e.,

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 $\mathbf {z}=\mathbf {0}$ on each side.

Unlike the transition probability of thermal light in Eq. (13), the integrand of $P^{f}(\mathbf {n},\mathbf {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 $\gamma _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].

5. 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^{\dagger }ZU]_{\mathbf {n,n}}$)/$\prod _{k=1}^{M}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^{\dagger }ZU]_{\mathbf {n,n}}$)/$\prod _{k=1}^{M}n_k!$ with respect to $z_i$ and $z_j$ and making $\mathbf {z=1}$, one can estimate the mode correlation between $i$ and $j$ modes, i.e., $\langle \hat {m}_i\hat {m}_j\rangle =\sum _{\mathbf {m=0}}^{\boldsymbol {\infty }}m_im_jP^f(\mathbf {n},\mathbf {m})$ with the bosonic number operator $\hat {m}_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^{\dagger }ZU]_{\mathbf {n},\mathbf {n}}$).