A general quantum algorithm for numerical integration

Quantum algorithms have shown their superiority in many application fields. However, a general quantum algorithm for numerical integration, an indispensable tool for processing sophisticated science and engineering issues, is still missing. Here, we first proposed a quantum integration algorithm suitable for any continuous functions that can be approximated by polynomials. More impressively, the algorithm achieves quantum encoding of any integrable functions through polynomial approximation, then constructs a quantum oracle to mark the number of points in the integration area and finally converts the statistical results into the phase angle in the amplitude of the superposition state. The quantum algorithm introduced in this work exhibits quadratic acceleration over the classical integration algorithms by reducing computational complexity from O(N) to O(√N). Our work addresses the crucial impediments for improving the generality of quantum integration algorithm, which provides a meaningful guidance for expanding the superiority of quantum computing.

functions of Holder class 22 .Heinrich proposed a quantum integration algorithm on the Lebesgue space, which proves the algorithm is optimal 23 .Rebentrost et al. presented an optical quantum multi-dimensional integration algorithm, which demonstrates quadratic acceleration over MC method 24 .These algorithms mentioned above have limitations on the type of integration function.In addition, some quantum and classical hybrid integration algorithms have also been proposed, such as Suzuki et al. raised a hybird integration method to simplify the circuit depth [25][26][27] .Moreover, some are controversial over acceleration capability, such as the MCI method based on quantum amplitude estimation (QAE) can bring quadratic acceleration 28,29 , but Kaneko et al. confirmed that the probability distribution of the initial state coding prepared by the Grover-Rudolph method does not reflect the quantum advantage 30 .
Existing classical integration algorithms face high computational complexity and quantum algorithms face low generality problems, as shown in Table 1.Here, we propose a general quantum integration algorithm (GQIA) that eliminates these shortages showing strong universality and operability.Firstly, our algorithm quantizes the classical Monte Carlo integration process and utilizes quantum superposition to possess exponential representation capabilities beyond classical methods.Secondly, to construct a quantum oracle with polynomial approximation of the integration function, use the parallelism of quantum to count the points on the integration region, and store them in the phase of the quantum state.Finally, amplitude amplification and phase estimation are used to obtain phase information with high probability and accuracy, whereupon the integration value is calculated.Compared with classical integration algorithms, GQIA demonstrates the quadratic acceleration and higher computational accuracy.Due to the use of Monte Carlo and polynomial approximations, GQIA has no limitations on the formal properties of integration functions, making it more versatile than other quantum algorithms.

Methods
For more comprehensive discussions of GQIA, we begin here by briefly reviewing some relevant results from quantum algorithm and quantum computation theory.

Amplitude amplification
The key idea of the amplitude amplification (AA) algorithm 24 originally came from an unordered database-search algorithm, known as Grover's quantum algorithm 5 .By amplifying the amplitude of a given pure state, the algorithm achieves the purpose of adjusting the measurement probability of the pure state.The specific implementation method of the algorithm is given as follows, first given an oracle Â and initial state |0� ⊗n (1) ψ� = Â|0� ⊗n = Â 00 . . .0� = cos θ/2|ψ 0 � + sin θ/2|ψ 1 � Table 1.Performance analysis of different integration algorithms. 1 N stands for the number of points of area S, N = 2 k 1+ k 1 , and has the same meaning as other n in the table.M represents the number of point in the area S 1 . 2α , d is real number. 3r , d is real number. 4δ is a constant that can be arbitrarily small. 31older function

Suitable Area Complexity Precision
Gauss forum [7][8][9] High precision with fewer nodes Taylor expansion 13,14 Each derivative of the function is known Romberg forum [7][8][9] less computation, high precision requirements Monte carlo 10 Complex function or the form is unknown Composite rule [7][8][9] Degree < 8 and integration interval is large Newton-cotes forum [7][8][9] Degree < 8 where θ ∈ [0, π ], |ψ 1 � is the target pure state and |ψ 0 � is the non-target pure state.The state vector |ψ� can be expressed as a parameterized vector (cosθ/2, sinθ/2) on the space stretched by two basis vectors |ψ 0 � and |ψ 1 � .The purpose is to promote the amplitude of the |ψ 1 � first, while reducing the amplitude of the non-target state, because the sum of the squares of the two probability amplitudes is 1.At this time, the measurement probability of the target state will also increase.This is done by flipping the quantum state |ψ� towards the direction of the target state |ψ 1 � .i.e., making θ as much larger as possible.
The overall operation Q is defined as where (2|00 . . .0��00 . . .0| − I) operation means to flipping the amplitude of all quantum states but keeping the amplitude of |00 . . .0� unchanged.By repeating this process many times, the quantum state |ψ� can be flipped to the target state |ψ 1 � with an angle of θ each time.What we need to pay attention to is that the real value of θ does not need to be known in advance.

Phase estimation
Quantum phase estimation is an algorithm for estimating the phase information of quantum states, and it is the core of many quantum algorithms.Two quantum registers are necessary in the phase estimation circuit, which the first register requires t qubits whose initial state is |0� and the value of t depends on two aspects: the number of digits for precision and successful probability of phase estimation result.
The initial state of the second register is the prepared state |u� and the number of qubits is as many as possible to store |u� .The process of phase estimation is roughly divided into three steps.First, the circuit begins with t Hadamard gates that are applied to all qubits in the first register, and the controlled-U operations are applied simultaneously on the second register, where the U-gate appears in consecutive integer powers of 2. The result- ing state is The second step of phase estimation is the application of inverse quantum Fourier transform on the first register, and this step can be accomplished in O t 2 stages.The third and final step of phase estimation is to measure the state in the first register and get an estimation of ϕ .The fianl state of the first register can be written as where φ denotes the estimator for ϕ when measured.

A General quantum integration algorithm (GQIA)
We define S as the area consisting of integration area S 1 and non-integration area S 2 .N and M denote the number of points in S and S 1 respectively (Fig. 1a).The core idea of GQIA is to convert a definite integration problem on a continuous interval [a, b] into the problem of getting the value of N*sin 2 (θ/2) with quantum computing.The phase θ that contained in the amplitude of a superposition state is obtained by phase estimation, as shown in Fig. 1b.To achieve integration calculation with quantum method, we need to address the following challenges.
The first problem is to encode and calculate integrable functions with quantum gate circuits.Supposing the function f (x) is integrable in the interval [a, b], it can be written as b a f (x) = F(b) − F(a) = I and I is a known constant.However, when functions are complex or cannot be formulated, they are usually approximated by the linear combination of function values under some discrete points in the interval.
The integration function can be quantized by quantum coding under discreate variables after selecting appropriate approximation methods.
The second problem is to get the distribution of points in the area S 1 .The point x that meets f (x) ≤ y is required, and the number is recorded as M. By constructing quantum gate circuits of computation and comparison, the distribution of these points in S 1 is acquired and stored in the amplitude of a superposition state.
The third problem is to get the ratio λ of distribution.There is no acceleration advantage to measure result from quantum circuit directly.Thus, by phase estimation circuit, we can acquire the phase θ with quadratic quantum advantage and the approximate value of the integration is S*sin 2 (θ/2).

Quantization of integration functions
Quantization of integration functions refers to the quantum representation of classical integration functions and the construction of the quantum circuits.The integrations over any interval [a, b] can be transformed into interval [0, 1], and the real number points (x, y) in the area S can be discretized by using k 1 and k 2 qubits, respectively.This means using bits of length k 1 to represent numbers between 0 and 1, similarly, to using bits of length k 2 to represent the value of the vertical axis y and these points evenly divide the integration interval S.
Thus, the number of discreate points is N = 2 k 1+ k 1 (Fig. 1a).The researches on quantum integration algorithms are rare, and one of the obstacles is to represent integrable functions with quantum coding.In numerical analysis, continuous and bounded functions can be polynomials approximated, such as Chebyshev approximation, best square approximation 8 , etc.If the function f(x) has n-order derivatives on interval [a, b] including x 0 , function f(x) can be approximated with a Taylor expansion at point x 0 , which the coefficients are the n-order derivatives of the function at x 0 .Otherwise, the polynomial coefficients can be determined by taking the function values from multiple points and using the method of undetermined coefficients, and the function is approximated by a simple polynomial, as shown in Eq. ( 6) If the first n + 1 terms are used to approximate the function, the precision of approximation is O(1/x n+1 ).Polynomial approximation is a linear combination of the power of variables and when the value of k 1 of variable discretization is determined, the power circuit is easy to construct, that is the reason for approximating integration functions with Eq. (6).For example, the quadratic power quantum circuit (Fig. 2) can be constructed with Table 2.

Construction of marking oracle
The area S is divided into S 1 and S 2 by the function curve, and only the points in S 1 are required.Here, a method is proposed to mark these points by constructing two quantum oracles.The first one is U f computing f(x), which f (x) = n i=0 a i x i is the first n terms of the approximate polynomial.Supposing r denotes the number of auxiliary qubits for this oracle, including q qubits representing the superposition output of f (x) under the initial input f u .The oracle needs k 1 + r qubits in that way, where k 1 stands for the number of qubits for variable x, and |φ u � is the output of the auxiliary qubits.
(5) www.nature.com/scientificreports/ The other oracle is to compare f(x i ) with y ij for each x i .The point (x i , y ij ) that meets f (x i ) ≤ y ij needs to be marked, and a comparison circuit including n CGC gates and n ICGC gates can make it 32 , which CGC and ICGC gates are quantum circuits composed of CNOT and Toffoli gates, and the CGC gate and its inverse ICGC are constructed with two CNOT gates and one Toffoli gate in different orders respectively (Fig. 3).Thus, the number of qubits is 2n + 2.
where k 2 is the number of qubits required to represent y, |y� = 1 |u� .If q is not equal to k 2 , 2 + q − k 2 auxiliary qubits |0� ⊗2+|q−k 2| are needed for this (7)   www.nature.com/scientificreports/comparator. is a real number between 0 and 1 evaluated by measurement, which representing the proportion of points that meet the comparison criteria and corresponding to ampulitude of the auxiliary bits, when the output is 1.Thus, the comparison results are obtained by evaluating the value of , and the necessary qubit number of marking oracle is To sum up, the marking oracle U of GQIA is made up with computation and comparation, which can be recorded as U = U cmp U f .Different from direct measurement, the GQIA obtains the results by phase estimation, which brings quadratic acceleration 33 .
Using Q repeated j times for the quantum state | � gives Similarly, from Eq. ( 8) we can get a quantum state | � and a unitary operator Q as where In the AA algorithm, the unitary operator Q is equivalent to rotate the superposition states by angle θ , which can be expressed as a 2*2 dimensional matrix in the single bit case It is easy to get the eigenvalues of Q from following formula The eigenvalues are γ 1 = e iθ and γ 2 = e i(2π−θ) , either one is feasible because θ has a small value and 2π − θ is large, making it easy to observe.In this article, it may be assumed that the value is θ and we take the case of γ 1 .Phase estimation is to get s in Q|ψ� = e 2πis |ψ� , where s = θ 2π .Supposing that we denote s with t qubits, that is, s = 0.s 1 s 2 • • • s t .We can obtain the result of with inverse quantum Fourier transform(Eq.19).
Therefore, s and can be obtained by measuring the first t qubits in phase estimation circuit, and the integra- tion can be approximated with 2 • N.

Figure 1 .
Figure 1.Diagram of quantum integration algorithm GQIA.(a) The classical monte carlo integration(MCI), where S 1 represents integration area, S 2 stands for non-integration area.M stands for the number of points in the area S 1 , and N is the total points in S 1 + S 2 .(b) Quantum circuit diagram of GQIA, Step1-Step5 corresponds to quantization of the MCI.

Figure 3 .
Figure 3. Structure and detailed circuits of comparator U cmp .