Universal bound on sampling bosons in linear optics and its computational implications

ABSTRACT In linear optics, photons are scattered in a network through passive optical elements including beam splitters and phase shifters, leading to many intriguing applications in physics, such as Mach–Zehnder interferometry, the Hong–Ou–Mandel effect, and tests of fundamental quantum mechanics. Here we present the fundamental limit in the transition amplitudes of bosons, applicable to all physical linear optical networks. Apart from boson sampling, this transition bound results in many other interesting applications, including behaviors of Bose–Einstein condensates (BEC) in optical networks, counterparts of Hong–Ou–Mandel effects for multiple photons, and approximating permanents of matrices. In addition, this general bound implies the existence of a polynomial-time randomized algorithm for estimating the transition amplitudes of bosons, which represents a solution to an open problem raised by Aaronson and Hance (Quantum Inf Comput 2012; 14: 541–59). Consequently, this bound implies that computational decision problems encoded in linear optics, prepared and detected in the Fock basis, can be solved efficiently by classical computers within additive errors. Furthermore, our result also leads to a classical sampling algorithm that can be applied to calculate the many-body wave functions and the S-matrix of bosonic particles.

Decision problems contain problems that has a definite answer, either yes or no, from an (classical or quantum) algorithm. The time complexity of an algorithm is measured by the number of elementary steps it takes to reach the answer, in the worst-case scenario. Specifically, the class of problems that can be solved by a classical computer with a polynomial number O(n k ) of steps is labeled as "polynomial time" (P). Here n is the size of the input of the problem. These problems are often regarded as "easy" problems for classical computers. It is because, for most of the practical problems in P, the degree of the polynomial k is usually not very large. Traditionally, a more precise definition of P is stated in terms of the language-theoretic framework: a language L (or decision problem) is in the complexity class P, if and only if there exists an algorithm (or classical circuit C n (x)) that runs in a polynomial time such that, (i) for all x ∈ L, the circuit outputs "yes", i.e., C n (x) = 1, and (ii) for all x / ∈ L, the circuit outputs "no", i.e., C n (x) = 0.
So far, we have discussed the idea of deterministic computation only; many more problems can be solved efficiently, if we are allowed to use randomness as well, which allows mistakes to occur occasionally. In this way, we might be able to deduce the correct answer by repeating the nondeterministic algorithm multiple times. More precisely, a language L is in the complexity class "bounded-error probabilistic polynomial time" (BPP), if and only if there exists a probabilistic classical circuit C n (x) that runs in a polynomial time, such that (i) for all x ∈ L, the circuit outputs "yes" with a probability larger than or equal to 2/3, i.e., Pr (C n (x) = 1) 2/3, and (ii) for all x / ∈ L, the circuit outputs "no" with a probability less than or equal to 1/3, i.e., Finally, the class of problems that can be solved efficiently by a quantum computer, namely "bounded-error probabilistic polynomial time" (BQP), can be defined in a similar way. A language L is in BQP if and only if there exists a polynomial-time quantum circuits Q n (x), which takes n qubits and output 1 bit, such that (i) for all x ∈ L, Pr(Q n (x) = 1) ≥ 2/3 , and

BACKGROUND ON COMPLEXITY CLASS BQP AND THE LINEAR-OPTICS VERSION
Let us consider a concrete example. Suppose we apply a quantum computer to simulate the time evolution of a many-spin system, initialized in the state |x 1 x 2 x 3 ...〉 with x i ∈ {0, 1}. Then, a quantum circuit is applied to simulate the time operator, i.e., U ≈ e −iHt . At the end, one would be interested in certain observables, e.g., the two-point correlation function between the first two Note that we can always express it in terms of the probabilities, P (a 1 a 2 ) ≡ a 3 ··· |〈a 1 a 2 a 3 · · ·| U |000 · · · 0〉| 2 , where a i ∈ {0, 1}, and the summation is over all indices other than a 1 and a 2 . In other words, we can express the correlation function in the following way: 〈σ 1 z σ 2 z 〉 = P (00) + P (11) − P (01) − P (10), and can focus on each of them separately.
For example, let us consider P (00) and express it through P (00) = 0.p 1 p 2 p 3 · · · for p i ∈ {0, 1}. Here x = x 1 x 2 x 3 · · · represents one of the possible inputs to the quantum circuit. The language L contains all the the strings such that p 1 = 1. The quantum algorithm can then be repeated until error is reduced to a small value. Then, the whole procedure is repeated for different p i 's, and for other P 's. In addition, one may also be interested in the probability of a particular outcome |y〉. In this case, we have Pr (Q n (x) = 1) = |〈y| U x |0 ⊗n 〉| 2 .
In the case of quantum optics, the input string x can be represented by the different initial states |t 1 t 2 t 3 · · · t m 〉, where |t k 〉 ≡ (t k !) −1/2 (a † k ) t k |vac〉. The quantum circuit is now replaced with a linear optical circuit, denoted by Q LO m , which can be decomposed into a polynomial number of the elementary linear-optical components. The unitary matrix associated with the optical circuit is denoted by U LO x . Furthermore, the language L now contains all of the strings such that the optical circuit outputs specific state |s 1 s 2 s 3 · · · s m 〉, which means that, we are interested in the following probability, Pr Q LO m (x) = 1 = 〈s 1 s 2 · · · s m | U LO x |t 1 t 2 · · · t m 〉 2 . One of the main goals of this work is to show that such a probability can in fact be calculated through a classical non-deterministic algorithm, making the decision problems for linear optics a problem in BPP.

BACKGROUND ON ERRORS IN DECISION VERSUS SAMPLING PROBLEMS
Sampling problems are very different from decision problems. Each sampling problem is associated with a probability distribution D x for an input x. A sampling problem is solved by an algorithm, if it can produce an approximating distribution D est In particular, D = 1. As the number of D x is exponentially many, 2 n for n-bit strings, typically, the values of each member in D x is exponentially small. If one is satisfied only with an additive error with an algorithm, it may just output zero, which tells us nothing. In this case, we will need to increase the accuracy for each probability; instead of additive error, one requires an accuracy within a multiplicative error. An estimation A est of a quantity A to within an multiplicative error means that A (1 − ε) A est A (1 + ε), or |A est − A| Aε. Consequently, we have D est x − D x D x ε = ε and the sampling problem can be solved with such accuracy.
In this work, we are not interested in the sampling problem for linear optics, as boson sampling is now widely regarded as a computationally-hard problem. Instead, we are interested in knowing if the decision problems, to within additive errors, of linear optics is still hard. Our results suggest that there exist a polynomial-time classical algorithm for solving the decisions problems encoded in linear optics.

BACKGROUND ON ESTIMATORS
It is known that any m × m matrix W = (w i,j ) permanent can be calculated exactly with a scaling O(m 2 2 m ) using Ryser's formula. Glynn suggested a different algorithm requiring a similar computational cost that the Glynn's formula is given as a normalized form: where is Glynn's estimator, and x = (x 1 , x 2 , ..., x m ) ∈ {±1} m . Based on the Glynn's formula, Gurvits proposed a polynomial-time randomized algorithm to produce an approximation Perm (W ) to the value of the permanent, with an additive error ± W m , i.e., where W ≡ sup v∕ =0 W v / v .
The main idea of Gurvits is that one can convert Glynn's formula into an expectation value: An approximation of the permanent is obtained by randomly and uniformly picking T strings x k ∈ {±1} m , for k = 1, 2, ..., T , and evaluate the average value: Gly (x k ) .