# A general quantum algorithm for numerical integration

### 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**\({{\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.

$${\int }_{a}^{b} f(x)dx\approx {\sum }_{k=0}^{n} {A}_{k}f\left({x}_{k}\right)$$

(4)

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**\({{\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*_{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*.

$$x,y) = \left\{ {\left( {x_{i} ,y_{ij} } \right)|i = {1},{2}, \ldots ,{2}^{{k_{1} }} ; j = {1},{2}, \ldots ,{2}^{{k_{2} }} } \right\}$$

(5)

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)

$$f\left( x \right) = a_{0} + a_{{1}} * x + a_{{2}} * x^{{2}} + \cdots + a_{n} * x^{n} + R_{n} \left( x \right)$$

(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 \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*_{1} + *r* qubits in that way, where *k*_{1} stands for the number of qubits for variable *x*, and \(\left|{\phi }_{u}\right.\rangle\) is the output of the auxiliary qubits.

$$H^{{ \otimes k_{1} }} |0\rangle ^{{ \otimes k_{1} }} |0\rangle ^{{ \otimes r}} \mathop \to \limits^{{U_{f} }} \left( {\frac{1}{{\sqrt 2 }}} \right)^{{k_{1} }} \sum\limits_{{u = 0}}^{{2k_{1} – 1}} | u\rangle \left| {f_{u} } \right.\rangle _{q} \left| {\phi _{u} } \right.\rangle _{{r – q}}$$

(7)

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 2*n* + 2.

$$\begin{gathered} U_{{cmp}} U_{f} H^{{ \otimes k_{1} + k_{2} }} |0\rangle ^{{ \otimes k_{1} + k_{2} }} |0\rangle ^{{ \otimes r + 2 + \left| {q – k_{2} } \right|}} \hfill \\ \quad = U_{{cmp}} U_{f} |x\rangle |y\rangle |0\rangle ^{{ \otimes r + 2 + \left| {q – k_{2} } \right|}} \hfill \\ \quad = U_{{cmp}} |x\rangle \left| {f_{u} } \right.\rangle _{q} \left| {\phi _{u} } \right.\rangle _{{r – q}} |y\rangle |0\rangle ^{{ \otimes 2 + \left| {q – k_{2} } \right|}} \hfill \\ \quad = |x\rangle \left| {f_{u} } \right.\rangle _{q} \left| {\phi _{u} } \right.\rangle _{{r – q}} |y\rangle |0\rangle ^{{ \otimes 1 + \left| {q – k_{2} } \right|}} \sqrt {1 – \lambda } |0\rangle \hfill \\ \quad \quad + |x\rangle \left| {f_{u} } \right.\rangle _{q} \left| {\phi _{u} } \right.\rangle _{{r – q}} |y\rangle |0\rangle ^{{ \otimes 1 + \left| {q – k_{2} } \right|}} \sqrt \lambda |1\rangle \hfill \\ \end{gathered}$$

(8)

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^{33}.

#### Extraction of results

For a quantum state in amplitude amplification (AA) algorithm^{29},

$$|\Psi \rangle =\mathcal{A}|0{\rangle }_{n}|0\rangle =\sqrt{1-\lambda }{\left|{\psi }_{0}\right.\rangle }_{n}|0\rangle +\sqrt{\lambda }{\left|{\psi }_{1}\right.\rangle }_{n}|1\rangle$$

(9)

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}}\),

$${\text{Q}}=-\mathcal{A}{{\text{S}}}_{0}{\mathcal{A}}^{-1}{ \, {\text{S}}}_{\chi }$$

(10)

\({\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.

$$|\Psi \rangle =\mathcal{A}|0{\rangle }_{n}|0\rangle ={\text{cos}}\theta /2{\left|{\psi }_{0}\right.\rangle }_{n}|0\rangle +{\text{sin}}\theta /2{\left|{\psi }_{1}\right.\rangle }_{n}|1\rangle$$

(11)

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

$${{\text{Q}}}^{j}|\Psi \rangle ={\text{cos}}((2j+1)\theta /2){\left|{\psi }_{0}\right.\rangle }_{n}|0\rangle +{\text{sin}}((2j+1)\theta /2){\left|{\psi }_{1}\right.\rangle }_{n}|1\rangle$$

(12)

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

$$\begin{gathered} |\tilde{\Psi }\rangle \;\; = {\mathbf{U}}H^{{ \otimes k_{1} + k_{2} }} |0\rangle ^{{ \otimes k_{1} + k_{2} }} |0\rangle ^{{ \otimes r + 2 + \left| {q – k_{2} } \right|}} \hfill \\ \quad \quad = |x\rangle \left| {f_{u} } \right.\rangle _{q} \left| {\phi _{u} } \right.\rangle _{{r – q}} |y\rangle |0\rangle ^{{ \otimes 1 + \left| {q – k_{2} } \right|\sqrt {1 – \lambda } |0\rangle }} \hfill \\ \quad \quad \quad \; + |x\rangle \left| {f_{u} } \right.\rangle _{q} \left| {\phi _{u} } \right.\rangle _{{r – q}} |y\rangle |0\rangle ^{{ \otimes 1 + \left| {q – k_{2} } \right|\sqrt \lambda |1\rangle }} \hfill \\ \quad \quad = |a\rangle |0\rangle + |b\rangle |1\rangle \hfill \\ \end{gathered}$$

(13)

$$\widetilde{\mathbf{Q}}=-\mathbf{U}H{\mathbf{S}}_{0}{H}^{-1}{\mathbf{U}}^{-1}{\mathbf{S}}_{\chi }$$

(14)

where

$$|a\rangle =|x\rangle {\left|{f}_{u}\right.\rangle }_{q}{\left|{\phi }_{u}\right.\rangle }_{r-q}|y\rangle |0{\rangle }^{\otimes 1+\left|q-{k}_{2}\right|}\sqrt{1-\lambda }=|x\rangle {\left|{f}_{u}\right.\rangle }_{q}{\left|{\phi }_{u}\right.\rangle }_{r-q}|y\rangle |0{\rangle }^{\otimes 1+\left|q-{k}_{2}\right|}{\text{cos}}(\theta /2)$$

(15)

$$|b\rangle =|x\rangle {\left|{f}_{u}\right.\rangle }_{q}{\left|{\phi }_{u}\right.\rangle }_{r-q}|y\rangle |0{\rangle }^{\otimes 1+\left|q-{k}_{2}\right|}\sqrt{\lambda }=|x\rangle {\left|{f}_{u}\right.\rangle }_{q}{\left|{\phi }_{u}\right.\rangle }_{r-q}|y\rangle |0{\rangle }^{\otimes 1+\left|q-{k}_{2}\right|}{\text{sin}}(\theta /2)$$

(16)

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

$$\widetilde{\mathbf{Q}}=\left(\begin{array}{cc}{\text{cos}}\theta & -{\text{sin}}\theta \\ {\text{sin}}\theta & {\text{cos}}\theta \end{array}\right)$$

(17)

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

$$|\gamma I-\widetilde{\mathbf{Q}}|=\left|\begin{array}{cc}\gamma -{\text{cos}}\theta & {\text{sin}}\theta \\ -{\text{sin}}\theta & \gamma -{\text{cos}}\theta \end{array}\right|=(\gamma -{\text{cos}}\theta {)}^{2}+{{\text{sin}}}^{2}\theta =0$$

(18)

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).

$$\frac{1}{\sqrt{{2}^{t}}}\sum_{j=0}^{{2}^{t}-1} {e}^{2\pi isj}|j\rangle \to |s\rangle$$

(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\).

$${\int }_{a}^{b} f(x)dx\approx \left({S}_{1}+{S}_{2}\right)\cdot {\lambda }^{2}={2}^{{k}_{1}+{k}_{2}}\cdot {{\text{sin}}}^{2}\left(\frac{\pi s}{{2}^{t}}\right)$$

(20)