Quantum Algorithm for Solving a Quadratic Nonlinear System of Equations

Cheng Xue CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei, Anhui, 230026, P. R. China Xiao-Fan Xu CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Yu-Chun Wu wuyuchun@ustc.edu.cn CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Hefei National Laboratory, Hefei, Anhui, 230088, P. R. China Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei, Anhui, 230026, P. R. China Guo-Ping Guo CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Hefei National Laboratory, Hefei, Anhui, 230088, P. R. China Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei, Anhui, 230026, P. R. China Origin Quantum Computing Company Limited, Hefei, Anhui, 230026, P. R. China

Abstract

Solving a quadratic nonlinear system of equations (QNSE) is a fundamental, but important, task in nonlinear science. We propose an efficient quantum algorithm for solving $n$ -dimensional QNSE. Our algorithm embeds QNSE into a finite-dimensional system of linear equations using the homotopy perturbation method and a linearization technique; then we solve the linear equations with a quantum linear system solver and obtain a state which is $\epsilon$ -close to the normalized exact solution of the QNSE with success probability $\Omega(1)$ . The complexity of our algorithm is $O({\rm polylog}(n/\epsilon))$ , which provides an exponential improvement over the optimal classical algorithm in dimension $n$ , and the dependence on $\epsilon$ is almost optimal. Therefore, our algorithm exponentially accelerates the solution of QNSE and has wide applications in all kinds of nonlinear problems, contributing to the research progress of nonlinear science.

I Introduction

Nonlinear equations appear in many natural and social sciences, such as fluid dynamics Anderson and Wendt (1995), biology Hobbie and Roth (2007), atmospheric dynamics Ghil and Childress (2012), and nonlinear vibration mechanics Bishop and Johnson (2011). By solving nonlinear equations, we understand various nonlinear phenomena, such as turbulence Wilcox et al. (2006), chaos Schuster (1984), and fractal Vicsek (1992). Most nonlinear equations have no or hardly solvable analytical solutions, such that many numerical methods have been developed Dennis Jr and Schnabel (1996). When the dimension of the nonlinear equations is large, solving the nonlinear equations with classical computers requires too many computational resources and may exceed the ability of classical computers. There is great demand for developing more efficient algorithms for solving nonlinear equations.

Quantum computing is a new model of computation which provides a quantum advantage in some specific problems Shor (1999); Grover (1996); Harrow et al. (2009). A typical example is solving linear equations, where quantum computing provides exponential acceleration Harrow et al. (2009). There are already many quantum algorithms for solving various linear equations, such as systems of linear equations Harrow et al. (2009); Childs et al. (2017); Subaş ı et al. (2019); Xu et al. (2021) and linear differential equations Clader et al. (2013); Berry (2014); Montanaro and Pallister (2016); Berry et al. (2017); Xin et al. (2020); Cao et al. (2013); Costa et al. (2019); Fillion-Gourdeau et al. (2017); Engel et al. (2019); Arrazola et al. (2019); Linden et al. (2022); Childs and Liu (2020); Childs et al. (2021); Nielsen and Chuang (2002). A natural idea is to use quantum computing to accelerate the solution of nonlinear equations. In recent years, some quantum algorithms for solving nonlinear differential equations have been proposed, such as nonlinear ordinary differential equations Leyton and Osborne (2008); Lloyd et al. ; Liu et al. (2021); Kyriienko et al. (2021); Xue et al. (2021a); Krovi ; Lin et al. ; Jin and Liu , the nonlinear Schrödinger equation Lubasch et al. (2020), and Navier-stokes equations Chen et al. (2022); Ljubomir (2022).

However, there are still few quantum algorithms for solving a system of nonlinear equations. A related algorithm proposed by Qian $et\ al$ . Qian et al. is based on Grover’s algorithm Grover (1996) and provides polynomial acceleration. Another related work is the quantum Newton’s method proposed by Xue $et\ al$ Xue et al. (2021b). The quantum Newton’s method is a quantum-classical hybrid algorithm constructed using quantum random access memory Giovannetti et al. (2008a, b); Kerenidis and Prakash and $l_{\infty}$ tomography Kerenidis et al. (2020). Influenced by the sample complexity of $l_{\infty}$ tomography, the quantum advantage of the quantum Newton’s method is verified only by numerical simulation. Whether there are more effective quantum algorithms for solving a system of nonlinear equations requires further research.

In this paper, we focus on a special kind of system of nonlinear equations, the quadratic nonlinear system of equations (QNSE). QNSE appears in all kinds of nonlinear problems, such as quadratic programming Mangasarian (1994), nonlinear element analysis Reddy (2014), and nonlinear differential equations Verhulst (2006). In specific, QNSE often appears when solving quadratic nonlinear differential equations, including the Navier-Stokes equations in fluid dynamics Anderson and Wendt (1995), the logistic equation in biology Hobbie and Roth (2007), and the Lorenz system in atmospheric dynamics Ghil and Childress (2012). QNSE also appears when solving nonlinear differential equations in which the degree of nonlinear polynomials is higher than two because these differential equations can be approximate to quadratic nonlinear differential equations Kerner (1981). Therefore, solving QNSE is a fundamental and important task, and the algorithm for accelerating the solution of QNSE has a wide range of applications.

We propose an effective quantum algorithm for solving $n$ -dimensional QNSE. In our algorithm, based on the homotopy perturbation method and a linearization technique, QNSE is embedded in a finite-dimensional system of linear equations. Then the condition number of the finite-dimensional system is optimized by splitting some subspaces of the finite-dimensional system. Next, we solve the system of linear equations with a quantum linear system solver Harrow et al. (2009); Childs et al. (2017) and obtain a state which is $\epsilon$ -close to the normalized exact solution of the QNSE with success probability $\Omega(1)$ , where $\Omega$ represents an asymptotic notation Cormen et al. (2022), which provides the asymptotic lower bound. The complexity of our algorithm is $O({\rm polylog}(n/\epsilon))$ , which provides an exponential improvement over the optimal classical algorithm in dimension $n$ , and the dependence on $\epsilon$ is almost optimal. Our algorithm places some constraints on the QNSE; it is suitable when the linear component of the QNSE is well conditioned and is dominant in QNSE.

This paper is organized as follows. Sec. II gives the definition of QNSE. The details of our algorithm are introduced in Sec. III. Sec. IV gives the main result of our algorithm. Then we give some applications of our algorithm in Sec. V. Finally, conclusions and discussions of our work are given in Sec. VI.

II Quadratic Nonlinear System of Equations

In this paper, QNSE is defined as

F_{0}+F_{1}{\bm{x}}+F_{2}{\bm{x}^{\otimes 2}}=0,

(1)

where $\bm{x}\in\mathbb{R}^{n}$ and $F_{i}\in\mathbb{R}^{n}\times\mathbb{R}^{n^{i}}$ . We also have the following assumptions and definitions for Eq. (1):

(1)

$F_{1}$ is invertible.
(2)

$F_{1}$ and $F_{2}$ are $s$ -sparse.

(3)

Parameters $\alpha$ , $\beta$ , and $R$ are defined as

		$\displaystyle\alpha:=\\|F_{1}^{-1}\\|\\|F_{0}\\|,$		(2)
		$\displaystyle\beta:=\\|F_{1}^{-1}\\|\\|F_{2}\\|,$
		$\displaystyle R:=\max\{4\alpha\beta,\\|F_{0}\\|\}.$

In this paper, if not specifically noted otherwise, $\|\cdot\|=\|\cdot\|_{2}$ .

(4)

Oracles $O_{F1}$ and $O_{F2}$ extract nonzero elements of $F_{1}$ and $F_{2}$ , respectively. $O_{F1}$ consists of $O_{F11}$ and $O_{F12}$ , and $O_{F2}$ consists of $O_{F21}$ and $O_{F22}$ , which are written as

	$\displaystyle O_{F11}\|i\rangle\|j\rangle=\|i\rangle\|f_{1}(i,j)\rangle,$		(3)
	$\displaystyle O_{F12}\|i\rangle\|j\rangle\|z\rangle=\|i\rangle\|j\rangle\|z\oplus(F_{1})_{i,j}\rangle,$		(4)
	$\displaystyle O_{F21}\|i\rangle\|j\rangle=\|i\rangle\|f_{2}(i,j)\rangle,$		(5)
	$\displaystyle O_{F22}\|i\rangle\|j\rangle\|z\rangle=\|i\rangle\|j\rangle\|z\oplus(F_{2})_{i,j}\rangle,$		(6)

where $f_{1}(i,j)$ and $f_{2}(i,j)$ represent the column index of the $j$ th nonzero entry of the $i$ th row of $F_{1}$ and $F_{2}$ respectively.

(5)

An oracle $O_{F0}$ is used to prepare the amplitude encoding of $F_{0}$ , which is written as

$O_{F0}|0\rangle=\frac{1}{\|F_{0}\|}\sum_{i=0}^{n-1}{F_{0,i}|i\rangle}.$ (7)

Formally, the problem to be solved is defined in Problem 1.

Problem 1.

Consider QNSE defined in Eq. (1); $\bm{x}^{*}$ is an exact solution of Eq. (1), and $|x\rangle$ is the amplitude encoding of $\bm{x}^{*}$ , which is written as

|x\rangle=\frac{1}{\|\bm{x}^{*}\|}\sum_{i=0}^{n-1}{x^{*}_{i}|i\rangle}.

(8)

Given oracles $O_{F0}$ , $O_{F1}$ , and $O_{F2}$ , the goal is to output a state $|\bar{x}\rangle$ such that $\||x\rangle-|\bar{x}\rangle\|\leq\epsilon$ .

III Quantum Homotopy Perturbation Method

In this section, we introduce the overall process of our algorithm. The process contains three steps:

(1)

Transform Eq. (1) into another kind of nonlinear equation with the homotopy perturbation method.
(2)

Embed the transformed nonlinear equations into a finite-dimensional system of linear equations, and solve the linear equations with a quantum linear system solver.
(3)

Measure some qubits of the output state of the quantum linear system solver and obtain the target state which represents a normalized approximate solution of Eq. (1).

The details of the whole process described above are introduced in the following three sections.

III.1 Homotopy perturbation method

The homotopy perturbation method is a classical method for solving nonlinear equations He (1999); Babolian et al. (2009); Chakraverty et al. (2019). The main process of the homotopy perturbation method for solving Eq. (1) is as follows. We construct the homotopy $\nu(p):[0,1]\to\mathbb{R}^{n}$ , which satisfies

H(\nu,p)=F_{0}+F_{1}{\nu}+pF_{2}{\nu^{\otimes 2}}=0.

(9)

With homotopy perturbation method, $\nu$ is written as

\nu=\nu_{0}+p\nu_{1}+p^{2}\nu_{2}+\dots+p^{c}\nu_{c},

(10)

where $c\in\mathbb{N}^{+}$ . Then substituting Eq. (10) into Eq. (9) and equating the terms with identical powers of $p$ , we have

\left\{\begin{aligned} &F_{1}\nu_{i}+F_{0}=0,&i=0,\\ &F_{1}\nu_{i}+F_{2}\sum_{j=0}^{j=i-1}{\nu_{j}\otimes\nu_{i-1-j}}=0,&i=1,2,\dots,c.\end{aligned}\right.

(11)

When $p=1$ , $\nu$ is an approximate solution of Eq. (1), and we define

\bm{\tilde{x}}=\nu(1)=\nu_{0}+\nu_{1}+\dots+\nu_{c}.

(12)

The error bound of $\bm{\tilde{x}}$ is analyzed in Lemma 1.

Lemma 1.

Let $\bm{x}^{*}$ represent an exact solution of Eq. (1); the approximate solution obtained with the homotopy perturbation method is $\tilde{\bm{x}}=\sum_{i=0}^{c}{\nu_{i}}$ . When $R<1$ and $c\geq\log_{1/R}{\frac{\alpha}{\epsilon(1-R)}}$ , where $\alpha$ and $R$ are defined in Eq. (2), $\tilde{\bm{x}}$ satisfies

\|\bm{x}^{*}-\tilde{\bm{x}}\|\leq\epsilon.

(13)

The proof of Lemma 1 is given in Appendix A. From Lemma 1 we see that when $R<1$ , Eq. (12) is convergent; the solution error $\epsilon$ decreases exponentially with $c$ .

III.2 Linear embedding

Then we embed Eq. (11) into a finite-dimensional system of linear equations

A\bm{y}=\bm{b},

(14)

where $\bm{y}=[\bm{y}_{0},\bm{y}_{1},\dots,\bm{y}_{c}]$ . The details of $\bm{y}$ , $A$ , and $\bm{b}$ are explained as follows.

First, ${\bm{y}}_{0}$ is defined as

{\bm{y}}_{0}=[\nu_{0}+\nu_{1}+\dots+\nu_{c}],

(15)

which means ${\bm{y}}_{0}$ contains an $n$ -dimensional vector and ${\bm{y}}_{0,0}$ represents the approximate solution $\tilde{\bm{x}}$ .

Then we substitute ${\bm{y}}_{0}$ into Eq. (11) and get

F_{1}{\bm{y}}_{0,0}+F_{2}\sum_{i=1}^{c}{\sum_{j=0}^{i-1}{\nu_{j}\otimes\nu_{i-1-j}}}=-F_{0}.

(16)

We consider $\nu_{j}\otimes\nu_{i-1-j}$ to be an independent element and define it as a component of ${\bm{y}}_{1}$ . Then ${\bm{y}}_{1}$ contains $c(c+1)/2$ $n^{2}$ -dimensional vectors, which is written as

{\bm{y}}_{1}=[\nu_{0}\otimes\nu_{0},\nu_{0}\otimes\nu_{1},\nu_{1}\otimes\nu_{0},\dots,\nu_{c-1}\otimes\nu_{0}].

(17)

${\bm{y}}_{2}$ is generated from ${\bm{y}}_{1}$ in a similar way. In general, ${\bm{y}}_{i}$ is generated from ${\bm{y}}_{i-1}$ . Repeating this process, we have ${\bm{y}}_{c}=[\nu_{0}^{\otimes{c+1}}]$ , and ${\bm{y}}_{c}$ satisfies

F_{1}^{\otimes c+1}{\bm{y}}_{c,0}=(-F_{0})^{\otimes c+1},

(18)

so ${\bm{y}}_{c}$ does not generate new elements. In summary, ${\bm{y}}_{i}$ can be written as

\bm{y}_{i}=\left\{\begin{array}[]{lll}&[\nu_{0}+\nu_{1}+\dots+\nu_{c}],&i=0,\\ &[\bm{y}_{i,0},\bm{y}_{i,1},\dots,\bm{y}_{i,\beta_{i}-1}],&1\leq i\leq c.\end{array}\right.

(19)

$\beta_{i}$ denotes the number of terms in $\bm{y}_{i}$ , and $\bm{y}_{i,j}$ represents the $j$ th item of $\bm{y}_{i}$ , which is written as $\bm{y}_{i,j}=\otimes_{k=0}^{i}{\nu_{a_{i,j,k}}}$ ; $a_{i,j,k}$ satisfies

\left\{\begin{array}[]{l}a_{i,j,k}\geq 0,\\ i+1\leq\sum_{k=0}^{i}{(a_{i,j,k}+1)}\leq c+1.\end{array}\right.

(20)

By Eq. (20), $\beta_{i}$ satisfies

\beta_{i}=\left\{\begin{array}[]{ll}1,&i=0,\\ \sum_{k=i}^{c}{\binom{k}{i}},&1\leq i\leq c.\end{array}\right.

(21)

We define

\vec{a}_{i,j}=[a_{i,j,0},a_{i,j,1},\dots,a_{i,j,i}].

(22)

The mapping $(i,j)\to\vec{a}_{i,j}$ is a one-to-one mapping, the time complexity to compute this mapping or its reverse is $O({\rm poly}(c))$ .

Next, we discuss the structure of matrix $A$ and vector ${\bm{b}}$ . We set $\vec{a}_{i,0}=[0,0,\cdots,0]$ ; then $\bm{y}_{i,0}=\nu_{0}^{\otimes i+1}$ and satisfies

F_{1}^{\otimes i+1}\bm{y}_{i,0}=(-F_{0})^{\otimes i+1}.

(23)

When $j\geq 1$ , $\vec{a}_{i,j}$ has nonzero elements, and we assume the first nonzero element is $a_{i,j,k}$ ; then we have

		$\displaystyle I_{n}^{\otimes k}\otimes F_{1}\otimes I_{n}^{\otimes i-k}\bm{y}_{i,j}+\nu_{a_{i,j,0}}\otimes\dots\otimes\nu_{a_{i,j,k-1}}$
		$\displaystyle\quad\otimes F_{2}(\sum_{l=0}^{a_{i,j,k}-1}{\nu_{l}\otimes\nu_{a_{i,j,k}-1-l}})\otimes\dots\otimes\nu_{a_{i,j,i}}=0,$		(24)

where $\nu_{a_{i,j,0}}\otimes\dots\otimes\nu_{l}\otimes\nu_{a_{i,j,k}-1-l}\otimes\dots\otimes\nu_{a_{i,j,i}}\in\bm{y}_{i+1}$ . Therefore, Eq. (14) can be expanded in the following form:

\left(\begin{array}[]{cccc}A_{0,0}&A_{0,1}&&\\ &A_{1,1}&\ddots&\\ &&\ddots&A_{c-1,c}\\ &&&A_{c,c}\end{array}\right)\left(\begin{array}[]{c}\bm{y}_{0}\\ \bm{y}_{1}\\ \vdots\\ \bm{y}_{c}\end{array}\right)=\left(\begin{array}[]{c}\bm{b}_{0}\\ \bm{b}_{1}\\ \vdots\\ \bm{b}_{c}\end{array}\right),

(25)

where $A_{i,i}$ is an ( $n^{i+1}\beta_{i}$ )-dimensional square matrix and $A_{i,i+1}$ is an ( $n^{i+1}\beta_{i}\times n^{i+2}\beta_{i+1}$ ) dimensional matrix. The elements of $A_{i,i}$ and $A_{i,i+1}$ are determined by Eqs. (23) and (III.2). With Eqs. (23) and (III.2), the expression of ${\bm{b}}$ can also be obtained; in detail, the $i$ th component of $\bm{b}$ is

\bm{b}_{i}=[(-F_{0})^{\otimes i+1},\bm{0},\dots,\bm{0}].

(26)

From Eq. (23), matrix $A$ contains the block matrix $F_{1}^{\otimes i},i=1,2,\cdots,c+1$ , which causes the condition number $\kappa_{A}$ of $A$ to increase exponentially with $c$ . We optimize $\kappa_{A}$ by splitting $F_{1}^{\otimes i+1}\bm{y}_{i,0}=(-F_{0})^{\otimes i+1}$ in Eq. (25) into

\left(\begin{array}[]{ccccc}B_{i,0}&I&&\\ &B_{i,1}&\ddots&\\ &&\ddots&I\\ &&&B_{i,i}\end{array}\right)\left(\begin{array}[]{c}\nu_{0}^{\otimes i+1}\\ F_{0}\otimes\nu_{0}^{\otimes i}\\ \vdots\\ F_{0}^{\otimes i}\otimes\nu_{0}\end{array}\right)=\left(\begin{array}[]{c}\bm{0}\\ \bm{0}\\ \vdots\\ -F_{0}^{\otimes i+1}\end{array}\right),

(27)

where $B_{i,j}=I_{n}^{\otimes j}\otimes F_{1}\otimes I_{n}^{\otimes i-j}$ . As a result, $\bm{y}_{i,0}$ is redefined as

\bm{y}_{i,0}=[\nu_{0}^{\otimes i+1},F_{0}\otimes\nu_{0}^{\otimes i},\dots,F_{0}^{\otimes i}\otimes\nu_{0}],

(28)

and the linear system $A\bm{y}=\bm{b}$ defined in Eq. (14) is adjusted accordingly. In later sections, Eq. (14) is defaulted to the adjusted linear system. The dimension of the linear system is

$\displaystyle N$	$\displaystyle=\sum_{i=0}^{c}{n^{i+1}(\beta_{i}+i)}$
	$\displaystyle=(n+1)^{c+1}-1-sn+\frac{cn^{c}(n-1)-n^{c}+1}{(n-1)^{2}}n^{2}$
	$\displaystyle\approx(n+1)^{c+1}+cn^{c+1}.$	(29)

Next, we solve Eq. (14) with the quantum linear system solver proposed in Childs et al. (2017); operations $O_{A}$ and $O_{b}$ are required. $O_{A}$ consists of $O_{A1}$ and $O_{A2}$ , which are defined as

\begin{array}[]{c}O_{A1}|i\rangle|j\rangle=|i\rangle|f_{a}(i,j)\rangle,\\ O_{A2}|i\rangle|j\rangle|z\rangle=|i\rangle|j\rangle|z\oplus A_{i,j}\rangle,\end{array}

(30)

where $f_{a}(i,j)$ represents the column index of the $j$ th nonzero entry of the $i$ th row of $A$ . $O_{b}$ is used to prepare the amplitude encoding of $\bm{b}$ , which is defined as

O_{b}|0\rangle=|b\rangle:=\frac{1}{\|\bm{b}\|}\sum_{i=0}^{c}{\sum_{j=0}^{\beta_{i}-1}{\|b_{i,j}\||i,j\rangle|b_{i,j}\rangle}},

(31)

where $|b_{i,j}\rangle$ is the amplitude encoding of $\bm{b}_{i,j}$ ; using Eqs. (26) and (27), $b_{i,j}$ is written as

\bm{b}_{i,j}=\left\{\begin{array}[]{cc}[\bm{0},\bm{0},\dots,-F_{0}^{\otimes i+1}],&0\leq i\leq c,j=0,\\ \bm{0},&0\leq i\leq c,1\leq j\leq\beta_{i}-1.\end{array}\right.

(32)

$O_{A}$ is constructed by querying $O_{F1}$ and $O_{F2}$ , and $O_{b}$ is constructed by querying $O_{F0}$ ; the query complexity is given in Lemma 2, and the proof of Lemma 2 is given in Appendix B.

Lemma 2.

The operations $O_{A1}$ and $O_{A2}$ defined in Eq. (30) can be constructed by querying $O_{F1}$ and $O_{F2}$ $O({\rm poly}(c))$ times; $O_{b}$ defined in Eq. (31) can be constructed by querying $O_{F0}$ $O({\rm poly}(c))$ times.

After running the quantum linear system solver, we obtain the output state $|\bm{y}\rangle$ , which is written as

|\bm{y}\rangle=\sum_{i=0}^{c}{\sum_{j=0}^{\beta_{i}-1}{|i,j\rangle|\bm{y}_{i,j}\rangle}}.

(33)

The query complexity of the quantum linear system solver is $O(s_{A}\kappa_{A}{\rm polylog}(s\kappa/\epsilon))$ , and the sparsity of matrix $A$ is

s_{A}=\frac{c(c+1)}{2}s.

(34)

$\kappa_{A}$ is decided by $c$ , $F_{1}$ , and $F_{2}$ . As shown in Lemma 3, when $c$ , $F_{1}$ , and $F_{2}$ satisfy some conditions, an upper bound of $\kappa_{A}$ is derived.

Lemma 3.

When $c$ , $F_{1}$ , and $F_{2}$ satisfy

\|F_{1}^{-1}\|(1+(c+1)\|F_{2}\|)<1,

(35)

the condition number $\kappa_{A}$ of matrix $A$ satisfies

\kappa_{A}\leq\frac{\kappa_{F_{1}}+1}{1-\|F_{1}^{-1}\|[1+(c+1)\|F_{2}\|]},

(36)

where $\kappa_{F_{1}}$ represents the condition number of $F_{1}$ .

The proof of Lemma 3 is given in Appendix C. From Lemma 3, $\kappa_{A}$ does not increase exponentially with $c$ and mainly depends on $\kappa_{F_{1}}$ .

III.3 Measurement

Finally, by measuring the first qubit register of $|\bm{y}\rangle$ to $|0,0\rangle$ , we obtain the target state $|\bm{\tilde{x}}\rangle=|\bm{y}_{0,0}\rangle$ , a normalized approximate solution of Eq. (1). This step is probabilistic; the success probability is

p=\frac{\|\bm{y}_{0,0}\|^{2}}{\|\bm{y}\|^{2}}.

(37)

A lower bound of $p$ is given in Lemma 4.

Lemma 4.

Consider the linear system $A\bm{y}=\bm{b}$ defined in Eq. (14). Let $\eta^{{}^{\prime}}=\|\bm{y}_{0,0}\|/R$ . When $\|F_{1}^{-1}\|<1$ and $R<\sqrt{2}/2$ , $p$ defined in Eq. (37) satisfies

p\geq\frac{(\eta^{{}^{\prime}})^{2}(1-2R^{2})}{(\eta^{{}^{\prime}})^{2}(1-2R^{2})+2}.

(38)

The proof of Lemma 4 is given in Appendix D. With Lemma 4, $p$ mainly depends on $\eta^{{}^{\prime}}$ and $1-2R^{2}$ .

IV Main Result

In this section, the main result of our work is given in Theorem 1.

Theorem 1.

Given the QNSE $F_{0}+F_{1}\bm{x}+F_{2}\bm{x}^{\otimes 2}=0$ defined in Sec. II with an exact solution $\bm{x}^{*}$ , let $\eta=\|\bm{x}^{*}\|/R$ , $c=\lceil\log_{\frac{1}{R}}{\frac{4\alpha}{\eta R\epsilon(1-R)}}\rceil$ , and $G=\|F_{1}^{-1}\|[1+(c+1)\|F_{2}\|]$ , where $\alpha$ and $R$ are defined in Eq. (2). When

\begin{array}[]{l}G<1,\\ R<\sqrt{2}/2,\end{array}

(39)

there exists a quantum algorithm with success probability $\Omega(1)$ to obtain a normalized quantum state $|\bm{\bar{x}}\rangle$ satisfying $\|\bm{\bar{x}}-\frac{\bm{x}^{*}}{\|\bm{x}^{*}\|}\|\leq\epsilon$ . The query complexity of the algorithm for the oracles of $F_{0}$ , $F_{1}$ , and $F_{2}$ is

O\left(\frac{\kappa_{F_{1}}s{\rm polylog}\left(\frac{\|F_{1}\|}{\epsilon\eta(1-G)(1-2R^{2})\|F_{2}\|}\right)}{\eta(1-G)\sqrt{1-2R^{2}}}\right).

(40)

The gate complexity is the query complexity multiplied by a factor of ${\rm polylog}\left(\frac{n\|F_{1}\|}{\epsilon\eta(1-G)(1-2R^{2})\|F_{2}\|}\right)$ .

The proof of Theorem 1 is given in Appendix E. Here we give some discussion of Theorem 1. The dependence on $n$ and $\epsilon$ of our algorithm is $O({\rm polylog}(n/\epsilon))$ . Compared with the classical algorithm, our algorithm provides exponential acceleration.

Exponential acceleration comes at the cost of stronger constraints on QNSE, which are listed in Eq. (39). The role of $G<1$ is to limit the condition number of matrix $A$ so that it does not increase exponentially with $c$ . $G$ can be written as

G=\kappa_{F_{1}}\frac{1+(c+1)\|F_{2}\|}{\|F_{1}\|}.

(41)

The other constraint $R<\sqrt{2}/2$ is used to bound the success probability of our algorithm. $R<\sqrt{2}/2$ also implies $R<1$ , which is the convergence condition of our algorithm. $R$ mainly depends on $4\alpha\beta$ because when $\|F_{0}\|\geq\sqrt{2}/2$ , $x$ can be rescaled to $\zeta x$ with a suitable constant $\zeta$ which keeps $4\alpha\beta$ unchanged and makes $\|F_{0}\|<\sqrt{2}/2$ . Notice that

4\alpha\beta=4\|F_{1}^{-1}\|^{2}\|F_{2}\|\|F_{0}\|=4\kappa_{F_{1}}^{2}\frac{\|F_{0}\|}{\|F_{1}\|}\frac{\|F_{2}\|}{\|F_{1}\|}.

(42)

From Eqs. (41) and (42), $G$ and $R$ measure the condition number $\kappa_{F_{1}}$ and the dominance of the linear component $F_{1}$ in QNSE from different aspects. When the linear component $F_{1}$ is well conditioned and is dominant in QNSE, $G$ and $R$ are small, and our algorithm is efficient. “Well conditioned” means the condition number $\kappa_{F_{1}}$ is not large. “Dominant” means that the strength of the linear component $\|F_{1}\|$ in QNSE is much larger than that of other components, i.e., $\|F_{1}\|^{2}>>\|F_{0}\|\|F_{2}\|$ and $\|F_{1}\|>>1+(c+1)\|F_{2}\|$ .

V Application

In this section we give two applications of our algorithm.

The first application is a two-dimensional QNSE, which is defined as

\left\{\begin{array}[]{ll}8x_{0}-x_{1}-0.5x_{0}^{2}+0.5x_{0}x_{1}+0.2=0,\\ -x_{0}+8x_{1}-0.5x_{1}^{2}+0.5x_{1}x_{0}-0.2=0.\end{array}\right.

(43)

The corresponding $F_{0},F_{1}$ , and $F_{2}$ are written as

	$\displaystyle F_{0}=\left(\begin{array}[]{c}0.2\\ -0.2\end{array}\right),F_{1}=\left(\begin{array}[]{cc}8&-1\\ -1&8\end{array}\right),$		(48)
	$\displaystyle F_{2}=\left(\begin{array}[]{cccc}-0.5&0.5&0&0\\ 0&0&0.5&-0.5\end{array}\right).$		(51)

We set $c=2$ and compute the parameters $G$ and $R$ ,

	$\displaystyle G\approx 7.24\times 10^{-1}<1,$
	$\displaystyle R\approx 2.82\times 10^{-1}<\sqrt{2}/2.$		(52)

Then $\bm{y}=[\bm{y}_{0},\bm{y}_{1},\bm{y}_{2}]$ , and each component of $\bm{y}$ is

	$\displaystyle\bm{y}_{0}=[\nu_{0}+\nu_{1}+\nu_{2}],$
	$\displaystyle\bm{y}_{1}=[\nu_{0}\otimes\nu_{0},F_{0}\otimes\nu_{0},\nu_{0}\otimes\nu_{1},\nu_{1}\otimes\nu_{0}],$
	$\displaystyle\bm{y}_{2}=[\nu_{0}^{\otimes 3},F_{0}\otimes\nu_{0}\otimes\nu_{0},F_{0}\otimes F_{0}\otimes\nu_{0}].$		(53)

The corresponding matrix $A$ and vector $\bm{b}$ are

\displaystyle A=\left(\begin{array}[]{cccccccc}F_{1}&F_{2}&&F_{2}&F_{2}&&&\\ &F_{1}\otimes I_{2}&I_{4}&&&&&\\ &&I_{2}\otimes F_{1}&&&&&\\ &&&I_{2}\otimes F_{1}&&I_{2}\otimes F_{2}&&\\ &&&&F_{1}\otimes I_{2}&F_{2}\otimes I_{2}&&\\ &&&&&F_{1}\otimes I_{4}&I_{8}&\\ &&&&&&I_{2}\otimes F_{1}\otimes I_{2}&I_{8}\\ &&&&&&&I_{4}\otimes F_{1}\end{array}\right),

(62)

\bm{b}=[-F_{0},\bm{0}_{4},-F_{0}^{\otimes 2},\bm{0}_{4},\bm{0}_{4},\bm{0}_{8},\bm{0}_{8},-F_{0}^{\otimes 3}]^{T},

(63)

where $\bm{0}_{4}=[0,0,0,0]$ and $\bm{0}_{8}$ is similar. Then we solve $A\bm{y}=\bm{b}$ and obtain the component $\bm{y}_{0}$ ,

\bm{\tilde{x}}=\bm{y}_{0}=[-2.2151849674\times 10^{-2},2.2292943149\times 10^{-2}].

(64)

A numerical solution $\bm{x}^{*}$ of Eq. (43) is obtained with the classical algorithm, which is written as

\bm{x}^{*}=[-2.2151848573\times 10^{-2},2.2292944259\times 10^{-2}].

(65)

Then $\bm{\tilde{x}}$ satisfies

\|\bm{x}^{*}-\bm{\tilde{x}}\|\approx 1.56\times 10^{-9}.

(66)

The second application is a nonlinear boundary problem, which is defined as

u_{xx}-u^{2}-\delta x^{2}=0,u(0)=0,\ u(1)=0,

(67)

where $\delta=5\times 10^{-4}$ . Then we discretize the equation in the following way:

x_{i}=(i+1)h,i=0,1,\dots,n-1,h=1/(n+1).

(68)

We define $u_{i}=u(x_{i})$ , and $u_{i,xx}$ is written as

u_{i,xx}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{2h^{2}}.

(69)

Then we have the QNSE

		$\displaystyle-u_{i+1}+2u_{i}-u_{i-1}+2h^{2}u_{i}^{2}+2h^{2}\delta x_{i}^{2}=0,$
		$\displaystyle i=0,1,\dots,n-1,$		(70)

where $u_{-1}=1$ and $u_{n}=0$ . Equation (V) can be represented as

F_{0}+F_{1}{\bm{u}}+F_{2}{\bm{u}}^{\otimes 2}=0,

(71)

where

F_{0}=2\delta h^{2}[x_{0}^{2},x_{1}^{2},\dots,x_{n-1}^{2}]^{T},

(72)

F_{1}=\left(\begin{array}[]{cccc}2&-1&&\\ -1&2&\ddots&\\ &\ddots&\ddots&-1\\ &&-1&2\end{array}\right),

(73)

(F_{2})_{i,j}=\left\{\begin{array}[]{cc}2h^{2},&j=ni,i=0,1,\dots,n-1,\\ 0,&j\neq ni,i=0,1,\dots,n-1.\end{array}\right.

(74)

We set $n=100$ , $c=2$ , and $\zeta=1200$ and rescale ${\bm{u}}$ to ${\bm{w}}=\zeta{\bm{u}}$ ; ${\bm{\omega}}$ satisfies

\zeta^{2}F_{0}+\zeta F_{1}{\bm{w}}+F_{2}{\bm{w}}^{\otimes 2}=0.

(75)

The rescaled QNSE satisfies

	$\displaystyle G\approx 9.01\times 10^{-1}<1,$
	$\displaystyle R\approx 6.27\times 10^{-1}<\sqrt{2}/2.$		(76)

Then we solve Eq. (75) with our algorithm and have $\tilde{\bm{\omega}}$ ; then $\tilde{\bm{u}}=\frac{1}{\zeta}\tilde{\bm{\omega}}$ . Equation (71) is also solved with the classical algorithm, and the solution $\bm{u}^{*}$ is obtained; the error satisfies

\|\bm{u}^{*}-\bm{\tilde{u}}\|\approx 5.41\times 10^{-19}.

(77)

VI Conclusion and Discussion

In this paper, a quantum algorithm for solving QNSE was proposed. When focusing on the equation dimension $n$ and the solution error $\epsilon$ , the complexity of our algorithm is $O(\rm{poly}[\log(n/\epsilon)])$ . The process of solving QNSE with a classical computer is usually transformed into solving a system of linear equations iteratively Rheinboldt (1974). The conjugate gradient method is a widely used linear system solver with complexity $O(ns\kappa)$ Shewchuk et al. (1994), so compared with the classical algorithm, our algorithm provides an exponential improvement in dimension $n$ . The dependence on $\epsilon$ of our algorithm complexity is $O(\rm{poly}(\log[1/\epsilon)])$ , which is almost optimal.

In practice, with the outstanding performance in solving QNSE, our algorithm can be used as a subprogram to accelerate the process of computation in many nonlinear problems, such as quadratic programming Mangasarian (1994) and quadratic nonlinear differential equations Anderson and Wendt (1995); Hobbie and Roth (2007); Ghil and Childress (2012). Furthermore, our algorithm provides a different idea for solving a system of nonlinear equations with quantum computing, which could inspire more quantum algorithms for solving nonlinear equations.

Our algorithm considers only QNSE; it can also be generalized to a higher-order nonlinear system of equations. For example, an $m$ -order nonlinear system of equations $\sum_{j=0}^{m}{F_{j}\bm{x}^{\otimes j}}=0$ can be transformed into a finite-dimensional system of linear equations $Ay=b$ through a process similar to that introduced in Sec. III. The convergence condition, the expression of $Ay=b$ , etc., depend on $F_{j},\ j=0,1,\dots,m$ .

An open question is whether our algorithm can be optimized further. Notice that our algorithm has several restrictions on QNSE. As discussed in Sec. IV, the constraints of our algorithm are that the linear component $F_{1}$ is well conditioned and is dominant in QNSE, which limit the applications of our algorithm. In the future we will consider optimizing the constraints of our algorithm, thereby expanding the applications of our algorithm.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 12034018), and the Innovation Program for Quantum Science and Technology No. 2021ZD0302300.

Appendix A Proof of Lemma 1

Proof.

The exact solution $\bm{x}^{*}$ can be written as

\bm{x}^{*}=\sum_{i=0}^{\infty}{\nu_{i}}.

(78)

Then

\|\bm{x}^{*}-\tilde{\bm{x}}\|=\|\sum_{i=c+1}^{\infty}{\nu_{i}}\|.

(79)

Next, we give an upper bound of $\|\nu_{i}\|$ . From Eq. (11),

\|\nu_{0}\|=\|F_{1}^{-1}F_{0}\|\leq\alpha.

(80)

We assume $\|\nu_{i}\|\leq\gamma_{i}\beta^{i}\alpha^{i+1}$ , where $\gamma_{i}$ is a constant and $\gamma_{0}=1$ . From Eq. (11), $\gamma_{i}$ can be set as

\gamma_{i}=\sum_{j=0}^{i-1}{\gamma_{j}\gamma_{i-1-j}},\gamma_{0}=1.

(81)

Equation (81) represents the Catalan series Koshy (2008), and $\gamma_{i}$ satisfies

\gamma_{i}=\frac{1}{i+1}\binom{2i}{i}\approx\frac{4^{i}}{i^{3/2}\sqrt{\pi}}<4^{i}.

(82)

Thus,

\|\nu_{i}\|\leq\alpha(4\alpha\beta)^{i}\leq\alpha R^{i}.

(83)

Considering $R<1$ , we have

\|\sum_{i=c+1}^{\infty}{\nu_{i}}\|\leq\sum_{i=c+1}^{\infty}{\|\nu_{i}\|}\leq\frac{\alpha R^{c+1}}{1-R}.

(84)

Therefore, when $c\geq\log_{1/R}{\frac{\alpha}{\epsilon(1-R)}}$ , $\tilde{\bm{x}}$ satisfies $\|\bm{x}^{*}-\tilde{\bm{x}}\|\leq\epsilon$ .

Appendix B Proof of Lemma 2

Proof.

(1) Construction process of $O_{A1}$ and $O_{A2}$ .

In fact, the construction process of $O_{A1}$ and $O_{A2}$ is the process to compute nonzero blocks of matrix $A$ . We just need to compute the nonzero blocks of $A_{i,i}$ and $A_{i,i+1}$ , $i\in\{0,1,\dots,c\}$ .

First, we regard $A_{i,i}$ as a $\beta_{i}$ -dimensional diagonal block matrix; the $j$ th block is written as $A_{(i,j),(i,j)}$ . When $j=0$ , $A_{(i,j),(i,j)}$ is the matrix defined in Eq. (27). When $j\in\{1,2,\dots,\beta_{i}-1\}$ , we compute $\vec{a}_{i,j}$ defined in Eq. (22) and find $k$ , the label of $\vec{a}_{i,j}$ ’s first nonzero element. Then $A_{(i,j),(i,j)}=I_{n}^{\otimes k}\otimes F_{1}\otimes I_{n}^{\otimes i-k}$ .

Second, we regard $A_{i,i}$ as a ( $\beta_{i}\times\beta_{i+1}$ )-dimensional block matrix; one block is written as $A_{(i,j),(i+1,j^{{}^{\prime}})}$ . When $j=0$ , $A_{(i,j),(i+1,j^{{}^{\prime}})}=\bm{0}$ for $j^{{}^{\prime}}\in\{0,1,\dots,\beta_{i+1}-1\}$ . When $j\in\{1,2,\dots,\beta_{i}-1\}$ , we compute $\vec{a}_{i,j}$ defined in Eq. (22) and find $a_{i,j,k}$ , the first nonzero element of $\vec{a}_{i,j}$ . Then we compute the related $\vec{a}_{i+1,j^{{}^{\prime}}}$ from Eq. (III.2); the number of $\vec{a}_{i+1,j^{{}^{\prime}}}$ is $a_{i,j,k}$ . Next, we compute $j^{{}^{\prime}}$ from $\vec{a}_{i+1,j^{{}^{\prime}}}$ ; for these $j^{{}^{\prime}}$ , $A_{(i,j),(i+1,j^{{}^{\prime}})}=I_{n}^{\otimes k}\otimes F_{2}\otimes I_{n}^{\otimes i-k}$ .

Therefore, for $i=0,1,\dots,c$ , $j=0,1,\dots,\beta_{i}-1$ , we can compute the nonzero $A_{(i,j),(i,j)}$ and $A_{(i,j),(i+1,j^{{}^{\prime}})}$ , the time complexity is $O({\rm poly}(c))$ . The nonzero elements of $A_{(i,j),(i,j)}$ or $A_{(i,j),(i+1,j^{{}^{\prime}})}$ can be extracted by querying $O_{F1}$ or $O_{F2}$ .

By implementing the above computation process with a quantum circuit, $O_{A1}$ and $O_{A2}$ can be constructed directly. Notice that the expression of $A_{(i,j),(i,j)}$ or $A_{(i,j),(i+1,j^{{}^{\prime}})}$ is related to $i$ and $k$ , where $i,k\leq c$ . $O_{F1}$ and $O_{F2}$ are required for different $i$ and $k$ ; therefore the query complexity of the whole process to $O_{F1}$ and $O_{F2}$ is $O({\rm poly}(c))$ .

(2)Construction process of $O_{b}$ .

Now we give the preparation process of $|b\rangle$ defined in Eq. (31). From Eqs. (31) and (32), $|b\rangle$ can be simplified to

|b\rangle=\frac{1}{\|\bm{b}\|}\sum_{i=0}^{c}{\|F_{0}\|^{i+1}|i,0\rangle|b_{i,0}\rangle}.

(85)

To prepare $|b\rangle$ , we first prepare

|\psi\rangle=\frac{1}{M}\sum_{i=0}^{c}{\|F_{0}\|^{i+1}|i,0\rangle|0\rangle}.

(86)

Notice that $|b_{i,0}\rangle$ can be prepared by querying $O_{F0}$ $i+1$ times; we define

U=\sum_{i=0}^{c}{|i,0\rangle\langle i,0|\otimes U_{i}},

(87)

where $U_{i}|0\rangle=|b_{i,0}\rangle$ . Then we have

|b\rangle=U|\psi\rangle.

(88)

The query complexity of this process to $O_{F0}$ is $O({\rm poly}(c))$ .

Appendix C Proof of Lemma 3

We first give the following lemma.

Lemma 5.

Given an $n$ -dimensional invertible matrix $M$ , $i\in\mathbb{N}^{+}$ , matrix $P$ is defined as

P=\left(\begin{array}[]{ccccc}P_{0,0}&I&&&\\ &P_{1,1}&I&&\\ &&\ddots&\ddots&\\ &&&P_{i-1,i-1}&I\\ &&&&P_{i,i}\end{array}\right),

(89)

where $P_{j,j}=I_{n}^{\otimes j}\otimes M\otimes I_{n}^{\otimes i-j}$ . Then $P$ is invertible, and $P^{-1}$ satisfies

\|P^{-1}\|\leq\frac{\|M^{-1}\|(1-\|M^{-1}\|^{i+1})}{1-\|M^{-1}\|}.

(90)

Proof.

Let $Q=P^{-1}$ and consider $P$ and $Q$ to be $[(i+1)\times(i+1)]$ -dimensional block matrices, $Q_{i,j}$ is written as

$\displaystyle Q_{j,j-k}=\bm{0},$	$\displaystyle 0\leq j\leq i,0\leq j-k\leq i,$	(91)
$\displaystyle Q_{j,j}=P_{j,j}^{-1},$	$\displaystyle 0\leq j\leq i,$
$\displaystyle Q_{j,j+k}=-P_{j,j}^{-1}Q_{j+1,j+k},$	$\displaystyle 0\leq j\leq i,0\leq j+k\leq i.$

From Eq. (91), $Q$ is an upper triangular block matrix. $Q$ is split into $Q=\sum_{k=0}^{i}{\Gamma_{k}}$ , and $\Gamma_{k}$ contains the block $Q_{j,j+k}$ . From Eq. (91) and the definition of $P_{j,j}$ ,

\|\Gamma_{k}\|=\|M^{-1}\|^{k+1},k=0,1,\dots,i.

(92)

Therefore,

\|P^{-1}\|\leq\sum_{k=0}^{i}{\|\Gamma_{k}\|}\leq\frac{\|M^{-1}\|(1-\|M^{-1}\|^{i+1})}{1-\|M^{-1}\|}.

(93)

Then the proof of Lemma 3 is as follows.

Proof.

First, consider the upper bound of $\|A\|$ . By the definition of $A_{i,i}$ ,

\|A_{i,i}\|\leq\|F_{1}\|+1,\quad i=0,1,\dots,c.

(94)

As $A_{0,1}A_{0,1}^{T}=\frac{c(c+1)}{2}F_{2}F_{2}^{T}$ , $\|A_{0,1}\|$ satisfies

\displaystyle\|A_{0,1}\|=\sqrt{\frac{c(c+1)}{2}}\|F_{2}\|<(c+1)\|F_{2}\|.

(95)

When $i\geq 1$ , consider $A_{i,i+1}$ to be $[(\beta_{i}+i)\times(\beta_{i+1}+i+1)]$ -dimensional block matrix; each block has the form $I^{\otimes j}\otimes F_{2}\otimes I^{\otimes i-j}$ , $0\leq j\leq i$ . From the structure of $A_{i,i+1}$ , it has no more than $c$ nonzero block matrices in each row or column, so $A_{i,i+1}$ can be split into at most $c$ matrices with at most one nonzero block matrix in each row or column, which leads to

\|A_{i,i+1}\|\leq c\|F_{2}\|,\quad i\in[1,2,\dots,c-1].

(96)

Combining Eqs. (94), (95) and (96), we have

$\displaystyle\\|A\\|\leq$	$\displaystyle\\|{\rm diag}(A_{0,0},A_{1,1},\dots,A_{c,c})\\|$
	$\displaystyle\ +\\|{\rm diag}(A_{0,1},A_{1,2},\dots,A_{c-1,c})\\|$
$\displaystyle=$	$\displaystyle\max_{i=0}^{c}\{\\|A_{i,i}\\|\}+\max_{i=0}^{c-1}\{\\|A_{i,i+1}\\|\}$
$\displaystyle\leq$	$\displaystyle\\|F_{1}\\|+1+(c+1)\\|F_{2}\\|.$	(97)

Now we analyze the upper bound of $\|A^{-1}\|$ . Let $B=A^{-1}$ , and consider $A$ and $B$ to be $[(c+1)\times(c+1)]$ -dimensional block matrices. The block $B_{i,j}$ satisfies

		$\displaystyle B_{i,i}=A_{i,i}^{-1},\quad 0\leq i\leq c,$		(98)
		$\displaystyle B_{i,i+k}=-A_{i,i}^{-1}A_{i,i+1}B_{i+1,i+k},$
		$\displaystyle\quad\quad 0\leq i\leq c-1,0\leq i+k\leq c.$

$B$ can be split by the distance from the diagonal,

B=\sum_{k=0}^{c}{\Lambda_{k}}.

(99)

$\Lambda_{k}$ contains the block $B_{i,i+k}$ . $\|\Lambda_{0}\|$ satisfies

\|\Lambda_{0}\|\leq\max_{0\leq i\leq c}{\|A_{i,i}^{-1}\|}.

(100)

$A_{i,i}$ can be viewed as a block-diagonal matrix; the first block is the matrix described in Eq. (27), and other blocks are in the form $I_{n}^{\otimes j}\otimes F_{1}\otimes I_{n}^{\otimes i-j}$ . Therefore, from Lemma 5,

\|A_{i,i}^{-1}\|\leq\max\left\{\|F_{1}^{-1}\|,\frac{\|F_{1}^{-1}\|(1-\|F_{1}^{-1}\|^{i+1})}{1-\|F_{1}^{-1}\|}\right\}.

(101)

Notice that $\|F_{1}^{-1}\|(1+(c+1)\|F_{2}\|)<1$ implies $\|F_{1}^{-1}\|<1$ ; we have

\|\Lambda_{0}\|\leq\frac{\|F_{1}^{-1}\|}{1-\|F_{1}^{-1}\|}.

(102)

From Eqs. (95) and (96), $\|A_{i,i+1}\|\leq(c+1)\|F_{2}\|$ . Therefore,

\|A_{i,i}^{-1}A_{i,i+1}\|\leq\frac{\|F_{1}^{-1}\|}{1-\|F_{1}^{-1}\|}(c+1)\|F_{2}\|.

(103)

We define

\zeta=\frac{\|F_{1}^{-1}\|}{1-\|F_{1}^{-1}\|}(c+1)\|F_{2}\|.

(104)

From Eqs. (98), (103), and (104),

\|\Lambda_{k}\|\leq\zeta\|\Lambda_{k-1}\|,

(105)

then

\|B\|\leq\sum_{k=0}^{c}{\|\Lambda_{k}\|}\leq\frac{1-\zeta^{c+1}}{1-\zeta}\|\Lambda_{0}\|.

(106)

With the assumption $\zeta<1$ , $\|A^{-1}\|$ satisfies

	$\displaystyle\\|A^{-1}\\|=\\|B\\|$	$\displaystyle\leq\frac{1}{1-\zeta}\times\frac{\\|F_{1}^{-1}\\|}{1-\\|F_{1}^{-1}\\|}$
		$\displaystyle=\frac{\\|F_{1}^{-1}\\|}{1-\\|F_{1}^{-1}\\|(1+(c+1)\\|F_{2}\\|)}.$		(107)

Therefore, from Eqs. (Proof) and (Proof),

$\displaystyle\kappa_{A}$	$\displaystyle=\\|A\\|\times\\|A^{-1}\\|$
	$\displaystyle\leq\frac{\kappa_{F_{1}}+\\|F_{1}^{-1}\\|[1+(c+1)\\|F_{2}\\|]}{1-\\|F_{1}^{-1}\\|[1+(c+1)\\|F_{2}\\|]}$
	$\displaystyle\leq\frac{\kappa_{F_{1}}+1}{1-\\|F_{1}^{-1}\\|[1+(c+1)\\|F_{2}\\|]}.$	(108)

Appendix D Proof of Lemma 4

Proof.

First, reorder $\bm{y}$ as $\bm{y}=[\bm{y}_{0}^{{}^{\prime}},\bm{y}_{1}^{{}^{\prime}},\dots,\bm{y}_{c}^{{}^{\prime}},\bm{y}_{c+1}^{{}^{\prime}}]$ ; $\bm{y}_{i}^{{}^{\prime}}$ is divided into the following three cases:

(1)

When $i=0$ , $\bm{y}_{0}^{{}^{\prime}}=\bm{y}_{0}=[\sum_{j=0}^{c}{\nu_{j}}]$ .
(2)

When $1\leq i\leq c$ , the element of $\bm{y}_{i}^{{}^{\prime}}$ is denoted as $\nu_{a_{i,0}^{{}^{\prime}}}\otimes\nu_{a_{i,1}^{{}^{\prime}}}\otimes\dots\otimes\nu_{a_{i,k}^{{}^{\prime}}}$ , which satisfies

$\left\{\begin{aligned} &k\geq 1,\\ &a_{i,j}^{{}^{\prime}}\geq 0,\\ &\sum_{j=0}^{k}{(a_{i,j}^{{}^{\prime}}+1)}=i+1.\end{aligned}\right.$ (109)

Therefore, the number of elements in $\bm{y}_{i}^{{}^{\prime}}$ is $2^{i}-1$ . Each item $\bm{y}_{i,j}^{{}^{\prime}}$ in $\bm{y}_{i}^{{}^{\prime}}$ satisfies

$\|\bm{y}_{i,j}^{{}^{\prime}}\|\leq\alpha^{k+1}R^{i-k}\leq\|F_{0}\|^{k+1}R^{i-k}\leq R^{i+1}.$ (110)

Then

$\|\bm{y}_{i}^{{}^{\prime}}\|^{2}\leq(2^{i}-1)R^{2(i+1)}<(2R^{2})^{i}R^{2}.$ (111)

Since $R<\sqrt{2}/2$ ,

$\sum_{i=1}^{c}{\|\bm{y}_{i}^{{}^{\prime}}\|^{2}}\leq\frac{2R^{4}}{1-2R^{2}}.$ (112)

(3)

When $i=c+1$ , $\bm{y}_{c+1}^{{}^{\prime}}$ contains the elements generated in Eq. (28) except $\nu_{0}^{\otimes i+1}$ for $i=1,2,\dots,c$ ; the elements are in the form

F_{0}^{\otimes j}\otimes\nu_{0}^{\otimes i-j},i=2,3,\dots,c+1,j=1,2,\dots,i-1,

(113)

which implies that $\bm{y}_{c+1}^{{}^{\prime}}$ contains $c(c+1)/2$ elements. Each item in $\bm{y}_{c+1}^{{}^{\prime}}$ satisfies $\|F_{0}^{\otimes j}\otimes\nu_{0}^{\otimes i-j}\|\leq\|F_{0}\|^{i}$ . Then $\bm{y}_{c+1}^{{}^{\prime}}$ satisfies

$\displaystyle\\|\bm{y}_{c+1}^{{}^{\prime}}\\|^{2}$	$\displaystyle\leq\\|F_{0}\\|^{4}\sum_{j=1}^{c}{j\\|F_{0}\\|^{2(j-1)}}$
	$\displaystyle<\frac{\\|F_{0}\\|^{4}}{(1-\\|F_{0}\\|^{2})^{2}}$
	$\displaystyle\leq\frac{R^{4}}{(1-R^{2})^{2}}.$	(114)

Therefore, with Eq. (112), Eq. ((3)), and $\|\bm{y}_{0}^{{}^{\prime}}\|=\|\bm{y}_{0}\|=\|\bm{y}_{0,0}\|=\eta^{{}^{\prime}}R$ , $p$ satisfies

$\displaystyle p$	$\displaystyle=\frac{\\|\bm{y}_{0}^{{}^{\prime}}\\|^{2}}{\\|\bm{y}_{0}^{{}^{\prime}}\\|^{2}+\sum_{i=1}^{c}{\\|\bm{y}_{i}^{{}^{\prime}}\\|^{2}}+\\|\bm{y}_{c+1}^{{}^{\prime}}\\|^{2}}$
	$\displaystyle\geq\frac{(\eta^{{}^{\prime}})^{2}R^{2}}{(\eta^{{}^{\prime}})^{2}R^{2}+\frac{2R^{4}}{1-2R^{2}}+\frac{R^{4}}{(1-R^{2})^{2}}}$
	$\displaystyle\geq\frac{(\eta^{{}^{\prime}})^{2}}{(\eta^{{}^{\prime}})^{2}+\frac{2R^{2}}{1-2R^{2}}+2}$
	$\displaystyle\geq\frac{(\eta^{{}^{\prime}})^{2}(1-2R^{2})}{(\eta^{{}^{\prime}})^{2}(1-2R^{2})+2}.$	(115)

Appendix E Proof of Theorem 1

To prove Theorem 1, some lemmas in Berry et al. (2017) are used, which are as follows.

Lemma 6.

Let $|\psi\rangle$ and $|\varphi\rangle$ be two vectors such that $\||\psi\rangle\|\geq\alpha>0$ and $\||\psi\rangle-|\varphi\rangle\|\leq\beta$ . Then

\left\|\frac{|\psi\rangle}{\||\psi\rangle\|}-\frac{|\varphi\rangle}{\||\varphi\rangle\|}\right\|\leq\frac{2\beta}{\alpha}.

(116)

Lemma 7.

Let $|\psi\rangle=\alpha|0\rangle|\psi_{0}\rangle+\sqrt{1-\alpha^{2}}|1\rangle|\psi_{1}\rangle$ and $|\varphi\rangle=\beta|0\rangle|\varphi_{0}\rangle+\sqrt{1-\beta^{2}}|1\rangle|\varphi_{1}\rangle$ , where $|\psi_{0}\rangle$ , $|\psi_{1}\rangle$ , $|\varphi_{0}\rangle$ , and $|\varphi_{1}\rangle$ are unit vectors and $\alpha,\beta\in[0,1]$ . Suppose $\||\psi\rangle-|\varphi\rangle\|\leq\delta<\alpha$ . Then $\||\psi_{0}\rangle-|\varphi_{0}\rangle\|\leq\frac{2\delta}{\alpha-\delta}$ .

Then the proof of Theorem 1 is shown as follows.

Proof.

Construct the system of linear equations $A\bm{y}=\bm{b}$ defined in Eq. (14), and define

\tilde{\bm{x}}=\bm{y}_{0}=\sum_{i=0}^{c}{\nu_{i}}.

(117)

With Lemma 1 and our choice of $c$ ,

\|\bm{x}^{*}-\tilde{\bm{x}}\|\leq\frac{\epsilon\eta R}{4}.

(118)

With $\|\bm{x}^{*}\|=\eta R$ , Eq. (117), and Lemma 6,

\|\frac{\bm{x}^{*}}{\|\bm{x}^{*}\|}-\frac{\tilde{\bm{x}}}{\|\tilde{\bm{x}}\|}\|\leq\epsilon/2.

(119)

The normalized state of $\bm{y}$ is written as

|\bar{\bm{y}}\rangle=\sum_{l=0}^{c}{\alpha_{l}|l\rangle|\bar{\bm{y}}_{l}\rangle}.

(120)

From Eqs. (117) and (120),

\bar{\bm{y}}_{0}=\frac{\tilde{\bm{x}}}{\|\tilde{\bm{x}}\|}.

(121)

Then $A\bm{y}=\bm{b}$ is solved with the quantum linear system solver proposed in Childs et al. (2017), and the output state is written as

|\bm{y}^{{}^{\prime}}\rangle=\sum_{l=0}^{c}{\alpha_{l}^{{}^{\prime}}|l\rangle|y_{l}^{{}^{\prime}}\rangle}.

(122)

Define $\eta^{{}^{\prime}}=\|\tilde{\bm{x}}\|/R$ ; from Eq. (118), $\eta^{{}^{\prime}}$ satisfies

|\eta-\eta^{{}^{\prime}}|\leq\frac{\epsilon\eta}{4}.

(123)

The solution error of the quantum linear system solver is set as

\delta=\frac{\eta^{{}^{\prime}}\sqrt{5(1-2R^{2})}}{30}\epsilon,

(124)

From Theorem $5$ in Childs et al. (2017),

\||\bar{\bm{y}}\rangle-|\bm{y}^{{}^{\prime}}\rangle\|\leq\delta.

(125)

We define $|\bm{\bar{x}}\rangle$ as

|\bm{\bar{x}}\rangle:=|y_{0}^{{}^{\prime}}\rangle.

(126)

From Eqs. (120), (122), (125), and (126) and Lemma 7, $|\bm{\bar{x}}\rangle$ satisfies

\||\bar{y}_{0}\rangle-|\bm{\bar{x}}\rangle\|\leq\frac{2\delta}{\alpha_{0}-\delta}.

(127)

From Lemma 4,

\alpha_{0}\geq\sqrt{\frac{(\eta^{{}^{\prime}})^{2}(1-2R^{2})}{(\eta^{{}^{\prime}})^{2}(1-2R^{2})+4}}>\sqrt{\frac{(\eta^{{}^{\prime}})^{2}(1-2R^{2})}{5}}.

(128)

The second inequality in the above equation defaults to $(\eta^{{}^{\prime}})^{2}(1-2R^{2})<1$ . Combining Eqs. (124), (127), and (128), we have

\||\bar{y}_{0}\rangle-|\bm{\bar{x}}\rangle\|\leq\epsilon/2.

(129)

[Notice that when $(\eta^{{}^{\prime}})^{2}(1-2R^{2})\geq 1$ , $\alpha_{0}\geq 1/\sqrt{5}$ ; then $\delta$ is set as $\delta=C\epsilon$ , where $C<\sqrt{5}/20$ , and Eq. (129) can also be derived.]

From Eqs. (117), (119), (121), and (129), the solution error satisfies

\left\|\frac{|\bm{x}^{*}\rangle}{\|\bm{x}^{*}\|}-|\bm{\bar{x}}\rangle\right\|\leq\epsilon.

(130)

The solution error influences the success probability of our algorithm. We default to $\epsilon<0.1$ ; then from Es. (124) and (128), $\delta$ satisfies

\delta\leq\frac{\epsilon}{6}\alpha_{0}<\frac{1}{60}\alpha_{0}.

(131)

Then

\alpha_{0}^{{}^{\prime}}\geq\alpha_{0}-\delta\geq\frac{59}{60}\alpha_{0}.

(132)

With Eq. (123) and $\eta^{{}^{\prime}}\geq(1-\epsilon/4)\eta$ ,

	$\displaystyle p=$	$\displaystyle(\alpha_{0}^{{}^{\prime}})^{2}\geq\left(\frac{59}{60}\right)^{2}\frac{1}{5}(1-\epsilon/4)^{2}\eta^{2}(1-2R^{2})$
	$\displaystyle>$	$\displaystyle 0.18\eta^{2}(1-2R^{2}).$		(133)

Next, we analyze our algorithm complexity. From Lemma 2, $O_{A}$ defined in Eq. (30) can be constructed by querying $O_{F1}$ and $O_{F2}$ $O(\rm{poly}(c))$ times; $O_{b}$ defined in Eq. (31) can be constructed by querying $O_{F0}$ $O(\rm{poly}(c))$ times. From Eq. (III.2), the dimension of the matrix $A$ is $N\approx(n+1)^{c+1}+cn^{c+1}$ . The sparsity of $A$ is $s_{A}=\frac{c(c+1)}{2}s$ . From Lemma 3, $\kappa_{A}\leq\frac{\kappa_{F_{1}}+1}{1-G}$ . From Theorem $5$ in Childs et al. (2017), when solving the linear system $A\bm{y}=\bm{b}$ , the query complexity of $O_{A}$ and $O_{b}$ is

O(s_{A}\kappa_{A}{\rm polylog}(s_{A}\kappa_{A}/\delta)).

(134)

Substituting the relevant parameters into Eq. (134), the query complexity of $O_{F0}$ , $O_{F1}$ , and $O_{F2}$ is

O\left(\frac{\kappa_{F_{1}}s}{1-G}{\rm polylog}\left(\frac{\|F_{1}\|}{\epsilon\eta(1-G)(1-2R^{2})\|F_{2}\|}\right)\right).

(135)

From Eq. (Proof), the success probability of the algorithm is $O(\eta^{2}(1-2R^{2}))$ . Using amplitude amplification Brassard et al. (2002), we repeat the above procedure $O(\frac{1}{\eta\sqrt{1-2R^{2}}})$ times and obtain the state $|\bm{\bar{x}}\rangle$ with probability $\Omega(1)$ . Therefore, the query complexity of our algorithm is

O\left(\frac{\kappa_{F_{1}}s\cdot{\rm polylog}\left(\frac{\|F_{1}\|}{\epsilon\eta(1-G)(1-2R^{2})\|F_{2}\|}\right)}{\eta(1-G)\sqrt{1-2R^{2}}}\right),

(136)

and the gate complexity is the query complexity multiplied by a factor of ${\rm polylog}\left(\frac{n\|F_{1}\|}{\epsilon\eta(1-G)(1-2R^{2})\|F_{2}\|}\right)$ .

References

Anderson and Wendt (1995) John David Anderson and J Wendt, Computational fluid dynamics (Springer, New York, 1995).
Hobbie and Roth (2007) Russell K Hobbie and Bradley J Roth, Intermediate physics for medicine and biology (Springer, New York, 2007).
Ghil and Childress (2012) Michael Ghil and Stephen Childress, Topics in geophysical fluid dynamics: atmospheric dynamics, dynamo theory, and climate dynamics (Springer, New York, 2012).
Bishop and Johnson (2011) Richard Evelyn Donohue Bishop and Daniel Cowan Johnson, The mechanics of vibration (Cambridge University Press, Cambridge, 2011).
Wilcox et al. (2006) David C Wilcox et al., Turbulence modeling for CFD (DCW industries La Canada, CA, 2006).
Schuster (1984) H. G. Schuster, Deterministic Chaos. An Introduction (PhysikVerlag, Weinheim, 1984).
Vicsek (1992) Tam Vicsek, Fractal growth phenomena (World scientific, Singapore, 1992).
Dennis Jr and Schnabel (1996) John E Dennis Jr and Robert B Schnabel, Numerical methods for unconstrained optimization and nonlinear equations (Society for Industrial and Applied Mathematics, Philadelphia, 1996).
Shor (1999) Peter W Shor, “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,” SIAM review 41, 303–332 (1999).
Grover (1996) Lov K Grover, “A fast quantum mechanical algorithm for database search,” in Proceedings of the twenty-eighth annual ACM symposium on Theory of computing (ACM, New York, 1996) pp. 212–219.
Harrow et al. (2009) Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd, “Quantum Algorithm for Linear Systems of Equations,” Physical Review Letters 103, 150502 (2009).
Childs et al. (2017) Andrew M. Childs, Robin Kothari, and Rolando D. Somma, “Quantum Algorithm for Systems of Linear Equations with Exponentially Improved Dependence on Precision,” SIAM Journal on Computing 46, 1920–1950 (2017).
Subaş ı et al. (2019) Yi ğit Subaş ı, Rolando D. Somma, and Davide Orsucci, “Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing,” Phys. Rev. Lett. 122, 060504 (2019).
Xu et al. (2021) Xiaosi Xu, Jinzhao Sun, Suguru Endo, Ying Li, Simon C. Benjamin, and Xiao Yuan, “Variational algorithms for linear algebra,” Science Bulletin 66, 2181–2188 (2021).
Clader et al. (2013) B David Clader, Bryan C Jacobs, and Chad R Sprouse, “Preconditioned quantum linear system algorithm,” Physical review letters 110, 250504 (2013).
Berry (2014) Dominic W Berry, “High-order quantum algorithm for solving linear differential equations,” Journal of Physics A 47, 105301 (2014).
Montanaro and Pallister (2016) Ashley Montanaro and Sam Pallister, “Quantum algorithms and the finite element method,” Physical Review A 93, 032324 (2016).
Berry et al. (2017) Dominic W. Berry, Andrew M. Childs, Aaron Ostrander, and Guoming Wang, “Quantum algorithm for linear differential equations with exponentially improved dependence on precision,” Communications in Mathematical Physics 356, 1057–1081 (2017).
Xin et al. (2020) Tao Xin, Shijie Wei, Jianlian Cui, Junxiang Xiao, Iñigo Arrazola, Lucas Lamata, Xiangyu Kong, Dawei Lu, Enrique Solano, and Guilu Long, “Quantum algorithm for solving linear differential equations: Theory and experiment,” Physical Review A 101, 032307 (2020).
Cao et al. (2013) Yudong Cao, Anargyros Papageorgiou, Iasonas Petras, Joseph Traub, and Sabre Kais, “Quantum algorithm and circuit design solving the poisson equation,” New Journal of Physics 15, 013021 (2013).
Costa et al. (2019) Pedro CS Costa, Stephen Jordan, and Aaron Ostrander, “Quantum algorithm for simulating the wave equation,” Physical Review A 99, 012323 (2019).
Fillion-Gourdeau et al. (2017) François Fillion-Gourdeau, Steve MacLean, and Raymond Laflamme, “Algorithm for the solution of the Dirac equation on digital quantum computers,” Physical Review A 95, 042343 (2017).
Engel et al. (2019) Alexander Engel, Graeme Smith, and Scott E Parker, “Quantum algorithm for the Vlasov equation,” Physical Review A 100, 062315 (2019).
Arrazola et al. (2019) Juan Miguel Arrazola, Timjan Kalajdzievski, Christian Weedbrook, and Seth Lloyd, “Quantum algorithm for nonhomogeneous linear partial differential equations,” Physical Review A 100, 032306 (2019).
Linden et al. (2022) Noah Linden, Ashley Montanaro, and Changpeng Shao, “Quantum vs. classical algorithms for solving the heat equation,” Communications in Mathematical Physics (2022), 10.1007/s00220-022-04442-6.
Childs and Liu (2020) Andrew M. Childs and Jin-Peng Liu, “Quantum spectral methods for differential equations,” Communications in Mathematical Physics 375, 1427–1457 (2020).
Childs et al. (2021) Andrew M. Childs, Jin-Peng Liu, and Aaron Ostrander, “High-precision quantum algorithms for partial differential equations,” Quantum 5, 574 (2021).
Nielsen and Chuang (2002) Michael A. Nielsen and Isaac Chuang, “Quantum computation and quantum information,” American Journal of Physics 70, 558–559 (2002).
Leyton and Osborne (2008) Sarah K. Leyton and Tobias J. Osborne, “A quantum algorithm to solve nonlinear differential equations,” (2008), arXiv:0812.4423 [quant-ph] .
(30) Seth Lloyd, Giacomo De Palma, Can Gokler, Bobak Kiani, Zi-Wen Liu, Milad Marvian, Felix Tennie, and Tim Palmer, “Quantum algorithm for nonlinear differential equations,” arXiv:2011.06571 [quant-ph] .
Liu et al. (2021) Jin-Peng Liu, Herman Øie Kolden, Hari K. Krovi, Nuno F. Loureiro, Konstantina Trivisa, and Andrew M. Childs, “Efficient quantum algorithm for dissipative nonlinear differential equations,” Proceedings of the National Academy of Sciences 118, e2026805118 (2021).
Kyriienko et al. (2021) Oleksandr Kyriienko, Annie E. Paine, and Vincent E. Elfving, “Solving nonlinear differential equations with differentiable quantum circuits,” Phys. Rev. A 103, 052416 (2021).
Xue et al. (2021a) Cheng Xue, Yu-Chun Wu, and Guo-Ping Guo, “Quantum homotopy perturbation method for nonlinear dissipative ordinary differential equations,” New Journal of Physics 23, 123035 (2021a).
(34) Hari Krovi, “Improved quantum algorithms for linear and nonlinear differential equations,” arXiv:2202.01054 .
(35) Yen Ting Lin, Robert B Lowrie, Denis Aslangil, Yiğit Subaşı, and Andrew T Sornborger, “Koopman von Neumann mechanics and the Koopman representation: A perspective on solving nonlinear dynamical systems with quantum computers,” arXiv:2202.02188 .
(36) Shi Jin and Nana Liu, “Quantum algorithms for computing observables of nonlinear partial differential equations,” arXiv:2202.07834 .
Lubasch et al. (2020) Michael Lubasch, Jaewoo Joo, Pierre Moinier, Martin Kiffner, and Dieter Jaksch, “Variational quantum algorithms for nonlinear problems,” Physical Review A 101, 010301 (2020).
Chen et al. (2022) Zhao-Yun Chen, Cheng Xue, Si-Ming Chen, Bing-Han Lu, Yu-Chun Wu, Ju-Chun Ding, Sheng-Hong Huang, and Guo-Ping Guo, “Quantum approach to accelerate finite volume method on steady computational fluid dynamics problems,” Quantum Inf. Process. 21, 137 (2022).
Ljubomir (2022) Budinski Ljubomir, “Quantum algorithm for the Navier–Stokes equations by using the streamfunction-vorticity formulation and the lattice Boltzmann method,” Int. J. Quantum. Inf. 20, 2150039 (2022).
(40) Peng Qian, Wei-Cong Huang, and Gui-Lu Long, “A quantum algorithm for solving systems of nonlinear algebraic equations,” arXiv:1903.05608 .
Xue et al. (2021b) Cheng Xue, Yu-Chun Wu, and Guo-Ping Guo, “Quantum Newton’s method for solving the system of nonlinear equations,” SPIN 11, 2140004 (2021b).
Giovannetti et al. (2008a) Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone, “Quantum random access memory,” Phys. Rev. Lett. 100, 160501 (2008a).
Giovannetti et al. (2008b) Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone, “Architectures for a quantum random access memory,” Phys. Rev. A 78, 052310 (2008b).
(44) Iordanis Kerenidis and Anupam Prakash, “Quantum recommendation systems,” arXiv:1603.08675 .
Kerenidis et al. (2020) Iordanis Kerenidis, Jonas Landman, and Anupam Prakash, “Quantum algorithms for deep convolutional neural networks,” in International Conference on Learning Representations (ICLR, Addis Ababa, Ethiopia, 2020).
Mangasarian (1994) Olvi L Mangasarian, Nonlinear programming (Society for Industrial and Applied Mathematics, Philadelphia, 1994).
Reddy (2014) Junuthula Narasimha Reddy, An Introduction to Nonlinear Finite Element Analysis: with applications to heat transfer, fluid mechanics, and solid mechanics, 2nd ed. (Oxford University Press, Oxford, 2014).
Verhulst (2006) Ferdinand Verhulst, Nonlinear differential equations and dynamical systems (Springer, Heidelberg, 2006).
Kerner (1981) Edward H. Kerner, “Universal formats for nonlinear ordinary differential systems,” Journal of Mathematical Physics 22, 1366–1371 (1981).
Cormen et al. (2022) Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein, Introduction to algorithms (MIT press, Cambridge, MA, 2022).
He (1999) Ji-Huan He, “Homotopy perturbation technique,” Computer Methods in Applied Mechanics and Engineering 178, 257–262 (1999).
Babolian et al. (2009) E. Babolian, A. Azizi, and J. Saeidian, “Some notes on using the homotopy perturbation method for solving time-dependent differential equations,” Math. Comput. Modell. 50, 213–224 (2009).
Chakraverty et al. (2019) Snehashish Chakraverty, Nisha Mahato, Perumandla Karunakar, and Tharasi Dilleswar Rao, Advanced numerical and semi-analytical methods for differential equations (Wiley, Hoboken, NJ, 2019).
Rheinboldt (1974) Werner C. Rheinboldt, Methods for Solving Systems of Nonlinear Equations (Society for Industrial and Applied Mathematics, Philadelphia, 1974).
Shewchuk et al. (1994) Jonathan Richard Shewchuk et al., “An introduction to the conjugate gradient method without the agonizing pain,” (1994).
Koshy (2008) Thomas Koshy, Catalan numbers with applications (Oxford University Press, Oxford, 2008).
Brassard et al. (2002) Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp, “Quantum amplitude amplification and estimation,” Contemp. Math. 305, 53–74 (2002).

		$\displaystyle\alpha:=\\|F_{1}^{-1}\\|\\|F_{0}\\|,$		(2)
		$\displaystyle\beta:=\\|F_{1}^{-1}\\|\\|F_{2}\\|,$
		$\displaystyle R:=\max\{4\alpha\beta,\\|F_{0}\\|\}.$

$\displaystyle\\|A\\|\leq$	$\displaystyle\\|{\rm diag}(A_{0,0},A_{1,1},\dots,A_{c,c})\\|$
	$\displaystyle\ +\\|{\rm diag}(A_{0,1},A_{1,2},\dots,A_{c-1,c})\\|$
$\displaystyle=$	$\displaystyle\max_{i=0}^{c}\{\\|A_{i,i}\\|\}+\max_{i=0}^{c-1}\{\\|A_{i,i+1}\\|\}$
$\displaystyle\leq$	$\displaystyle\\|F_{1}\\|+1+(c+1)\\|F_{2}\\|.$	(97)

	$\displaystyle\\|A^{-1}\\|=\\|B\\|$	$\displaystyle\leq\frac{1}{1-\zeta}\times\frac{\\|F_{1}^{-1}\\|}{1-\\|F_{1}^{-1}\\|}$
		$\displaystyle=\frac{\\|F_{1}^{-1}\\|}{1-\\|F_{1}^{-1}\\|(1+(c+1)\\|F_{2}\\|)}.$		(107)

$\displaystyle\kappa_{A}$	$\displaystyle=\\|A\\|\times\\|A^{-1}\\|$
	$\displaystyle\leq\frac{\kappa_{F_{1}}+\\|F_{1}^{-1}\\|[1+(c+1)\\|F_{2}\\|]}{1-\\|F_{1}^{-1}\\|[1+(c+1)\\|F_{2}\\|]}$
	$\displaystyle\leq\frac{\kappa_{F_{1}}+1}{1-\\|F_{1}^{-1}\\|[1+(c+1)\\|F_{2}\\|]}.$	(108)

$\displaystyle\\|\bm{y}_{c+1}^{{}^{\prime}}\\|^{2}$	$\displaystyle\leq\\|F_{0}\\|^{4}\sum_{j=1}^{c}{j\\|F_{0}\\|^{2(j-1)}}$
	$\displaystyle<\frac{\\|F_{0}\\|^{4}}{(1-\\|F_{0}\\|^{2})^{2}}$
	$\displaystyle\leq\frac{R^{4}}{(1-R^{2})^{2}}.$	(114)