A general quantum algorithm for numerical integration

A General quantum integration algorithm (GQIA)

We define S as the area consisting of integration area S₁ and non-integration area S₂. N and M denote the number of points in S and S₁ 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*\({{\text{sin}}}^{2}(\theta /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 \({\int }_{a}^{b}f(x)=F\left(b\right)-F\left(a\right)=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₁. 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₁ 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*\({{\text{sin}}}^{2}(\theta /2)\).

The components of GQIA

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₁ and k₂ qubits, respectively. This means using bits of length k₁ to represent numbers between 0 and 1, similarly, to using bits of length k₂ 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⁸, etc. If the function f(x) has n-order derivatives on interval (a, b) including x₀, function f(x) can be approximated with a Taylor expansion at point x₀, which the coefficients are the n-order derivatives of the function at x₀. 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ⁿ⁺¹). Polynomial approximation is a linear combination of the power of variables and when the value of k₁ 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₁ and S₂ by the function curve, and only the points in S₁ 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 \left(x\right)=\sum_{i=0}^{n}{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 \({\left|{f}_{u}\right.\rangle }\). The oracle needs k₁ + r qubits in that way, where k₁ stands for the number of qubits for variable x, and \(\left|{\phi }_{u}\right.\rangle\) is the output of the auxiliary qubits.

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³², 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\rangle ={\left(\frac{1}{\sqrt{2}}\right)}^{{k}_{2}}{\sum }_{\nu =0}^{{{2}{k}_{2}}-1} |\nu \rangle\), similarly, \(|x\rangle ={\left(\frac{1}{\sqrt{2}}\right)}^{{k}_{1}}{\sum }_{u=0}^{{{2}{k}_{1}}-1} |u\rangle\). If \(q\) is not equal to \({k}_{2},2+\left|q-{k}_{2}\right|\) auxiliary qubits \(|0{\rangle }^{\otimes 2+\left|q-{k}_{2}\right|}\) are needed for this comparator. \(\lambda\) 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 \(\lambda\), and the necessary qubit number of marking oracle is \({k}_{1}+{k}_{2}+r+2+\left|q-{k}_{2}\right|\).

To sum up, the marking oracle \(\mathbf{U}\) of GQIA is made up with computation and comparation, which can be recorded as \(\mathbf{U}={U}_{cmp}{U}_{f}\). Different from direct measurement, the GQIA obtains the results by phase estimation, which brings quadratic acceleration³³.

Extraction of results

For a quantum state in amplitude amplification (AA) algorithm²⁹,

where \({\left|{\psi }_{1}\right.\rangle }_{n}\) is the required target state, \({\left|{\psi }_{0}\right.\rangle }_{n}\) is non-target state. Defining \(\theta \in (0,\pi )\), so that \({{\text{sin}}}^{2}\theta /2=\lambda\) and a unitary operator \({\text{Q}}\),

\({\mathbf{S}}_{\chi }\) adds a negative phase before \({\left|{\psi }_{1}\right.\rangle }_{n}|1\rangle ,{\left|{\psi }_{0}\right.\rangle }_{n}|0\rangle\) keeps unchanged.

Using \({\text{Q}}\) repeated j times for the quantum state \(|\Psi \rangle\) gives

Similarly, from Eq. (8) we can get a quantum state \(|\widetilde{\Psi }\rangle\) and a unitary operator \(\widetilde{\mathbf{Q}}\) as

where

In the \({\text{AA}}\) algorithm, the unitary operator \(\widetilde{\mathbf{Q}}\) is equivalent to rotate the superposition states by angle \(\theta\), which can be expressed as a 2*2 dimensional matrix in the single bit case

It is easy to get the eigenvalues of \(\widetilde{\mathbf{Q}}\) from following formula

The eigenvalues are \({\gamma }_{1}={e}^{i\theta }\) and \({\gamma }_{2}={e}^{i(2\pi -\theta )}\), either one is feasible because \(\theta\) has a small value and \(2\pi -\theta\) is large, making it easy to observe. In this article, it may be assumed that the value is \(\theta\) and we take the case of \({\gamma }_{1}\). Phase estimation is to get \(s\) in \(\widetilde{\mathbf{Q}}|\psi \rangle ={e}^{2\pi is}|\psi \rangle\), where \(s=\frac{\theta }{2\pi }\). Supposing that we denote \(s\) with \(t\) qubits, that is, \(s=0.{s}_{1}{s}_{2}\cdots {s}_{t}\). We can obtain the result of with inverse quantum Fourier transform(Eq. 19).

Therefore, \(s\) and \(\lambda\) can be obtained by measuring the first \(t\) qubits in phase estimation circuit, and the integration can be approximated with \({\lambda }^{2}\cdot N\).