A Hybrid Quantum Solver for the Lorenz System

Abstract

We develop a hybrid classical–quantum method for solving the Lorenz system. We use the forward Euler method to discretize the system in time, transforming it into a system of equations. This set of equations is solved by using the Variational Quantum Linear Solver (VQLS) algorithm. We present numerical results comparing the hybrid method with the classical approach for solving the Lorenz system. The simulation results demonstrate that the VQLS method can effectively compute solutions comparable to classical methods. The method is easily extended to solving similar nonlinear differential equations.

1. Introduction

Dynamical systems theory, along with its numerical and simulated counterparts, has been employed to study various phenomena, such as weather patterns, population dynamics, economic trends, and the flow of chemical and biological elements [1]. Whether it is the macroscopic dynamics of a temperature–pressure system or the microscopic complexities of a nonlinear system with bifurcations and irreversible physics, mathematical modeling, based on dynamical systems theory, offers a valuable tool for analysis.

Dynamical systems may develop deterministic chaos, with the Lorenz system being a noticeable example. It was introduced by Lorenz in 1963 [2] as a simple nonlinear model of heat convection, and it stands out as one of the earliest attempts to capture atmospheric physics through a model consisting of three differential equations. It is a chaotic three-dimensional system that has been studied extensively, mostly numerically. For a known set of parameter values, the three-dimensional motion converges to the well-recognized butterfly-shaped attractor, which can be observed by solving the system numerically. However, no analytic proof of this fact is known.

In this paper, we focus on how quantum computing might help to improve the numerical solutions of this important nonlinear dynamical system. Unlike classical physics, quantum mechanics generally only provides probabilities for the outcomes of measurements. Quantum systems evolve in a linear fashion, and when a quantum system in superposition is measured, the act of measurement forces the system into a particular state, and the theory predicts the probabilities of obtaining various possible outcomes. This inherent probabilistic measurement of outcomes is a crucial feature of quantum uncertainty as opposed to uncertainty in a weather system which arise due to nonlinear behavior [3].

When we linearize the system by using the forward Euler method and solve it on a classical computer, the subsequent point is uniquely determined by the current point and is obtained by solving a system of linear equations. Classical chaos arises in systems governed by nonlinear equations with a positive Lyapunov coefficient, where small variations in the initial conditions can lead to drastically different outcomes over time. This sensitivity to the initial conditions is visible in solutions computed on a classical computer by using the forward Euler method.

When solving a single iteration of the Euler method by using a quantum algorithm, there are two sources of uncertainty. The initial state can only be prepared with a certain level of precision, and the solution computed by a quantum algorithm is inherently probabilistic, resulting in a random sample rather than a deterministic value. However, it is important to note that this randomness does not introduce new features or chaotic behavior into the system. The differences in trajectories between classical and quantum solutions arise from the inherent probabilistic nature of quantum algorithms, not from any additional characteristics of the Lorenz system itself. Both methods—classical and quantum—can provide approximate solutions, but the quantum method may do so more efficiently by exploiting quantum randomness, similar to how stochastic classical methods use randomness.

Suppose we replace the equation-solving step in the Euler method with a quantum algorithm; given the new source of uncertainty, the following questions arise: (i) Does the chaos in the system increase when a quantum subroutine is used in each iteration? (ii) Does the system still have the butterfly attractor (for the choice of parameters given above)? We answer these questions in the affirmative under the assumption that it is possible to prepare the state exactly and recover the solution exactly. These assumptions can only be justified with sufficient advances in quantum random-access memory (QRAM) and quantum tomography techniques [4,5].

One of the earliest attempts to solve nonlinear differential equations by using a quantum approach is due to Leyton [6]. They used the forward Euler method and multiple copies of an initial state, which were evolved according to the linear Euler system. This algorithm scaled poorly as a function of the time step, as shown in Berry [7], who went on to show that a far more efficient quantum procedure can be obtained for linear differential equations. For this, they computed the bounds on the condition number of the matrix. Liu [8] gave an efficient quantum algorithm for dissipative nonlinear differential equations.

Berry [9] gave an exponentially improved quantum algorithm for solving linear ordinary differential equations with constant coefficients. Although the paper primarily focuses on diagonalizable matrices, it can be extended to approximate solutions for non-diagonalizable matrices by using nearby diagonalizable matrices, as discussed in [10].

Krovi [10] presented a quantum algorithm for solving linear, inhomogeneous ordinary differential equations (ODEs). The algorithm shows improved gate and query complexity for specific diagonalizable classes of matrices and is extended to handle non-diagonalizable and singular matrices.

Lloyd [11] introduced a quantum algorithm called the “quantum nonlinear solver” which can solve nonlinear differential equations—a critical component of weather prediction models. Applications of these algorithms to weather prediction have been studied by Tennie [12]. This solver leverages multiple copies of a quantum state to simulate nonlinearity, potentially providing an advantage over classical methods. Initial tests of this algorithm on a simplified model have shown promising results, demonstrating its agreement with classical methods and ensemble averages for a single-particle system, $dx/dt=x-\alpha x^3$. The algorithm discretizes the time domain, similar to the forward Euler method. The set of equations can be represented in matrix form as $Ax=b$, with the state vector x obtained by using the Harrow–Hassidim–Lloyd (HHL) algorithm. The nonlinear forward Euler method is utilized in the case of nonlinear differential equations.

For the potential applications of quantum computing to climate change studies, see the recent review by Rahman [13]. The authors identified quantum principal component analysis (qPCA) and HHL as algorithms with quantum advantages that are crucial to analyzing climate models. It is difficult to leverage the advantages of quantum computing for large weather and climate datasets due to limited readout capacity and data accessibility challenges [12].

Armaos [14] explored the use of the Variational Quantum Eigensolver (VQE), a quantum computing algorithm, to analyze the Lorenz system. The main focus of their analysis was to determine the eigenvalues of the modified Jacobian matrix by using the VQE. This identified stable and unstable points of the Lorenz system.

Contributions: The main contribution of this work is a method, described in Section 2.2, for solving the Lorenz system of atmospheric convection. This method uses a hybrid quantum–classical algorithm. The technique involves discretizing the Lorenz system in time by using the forward Euler method. The resulting system of equations is solved by using the Variational Quantum Linear Solver (VQLS) described in Section 3. Section 4 is dedicated to the algorithm. The results of the simulations are reported in Section 5. The simulation of the VQLS method is computationally intensive on classical computers compared with classical algorithms that require inverting a matrix. However, it requires significantly fewer qubits (three) compared with other quantum algorithms, like the HHL algorithm, which requires a minimum of nine qubits for one iteration, making it potentially suitable for implementation on near-term noisy intermediate-scale quantum (NISQ) computers. Simulations demonstrate that the VQLS method effectively computes solutions comparable to classical methods. A detailed error analysis of the two methods, classical and quantum, based on Richardson extrapolation is in Section 5.5. Unlike previous quantum approaches [12] to computing atmospheric dynamics that used multiple copies of the state, this method simulates only a single time step, avoiding measurement errors that accumulate over time. The method developed in this paper is easily extended to solve other similar nonlinear differential equations. Results and limitations of the approach are discussed in Section 5.7.

2. Lorenz System

The Lorenz system consists of three coupled, nonlinear ordinary differential equations. The system was first introduced by Edward Lorenz in 1963 in a paper where he explored the underlying mechanisms of long-range weather prediction [2]. The equations describe the rate of change of three quantities over time, often interpreted physically as the movement of a fluid cell in a larger circulation. The Lorenz equations are as follows:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dx}{dt}}& =& \sigma (y-x),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dy}{dt}}& =& x(\rho -z)-y,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dz}{dt}}& =& xy-\beta z,\hfill \end{array}$$

(3)

where the following apply:

• x, y, and z represent the state variables of the system. These variables can be thought of as proportional to the intensity of the convective motion, the temperature difference between ascending and descending currents, and the deviation of the vertical temperature profile from linearity, respectively.
• $\sigma$ is the Prandtl number, representing the ratio of viscous diffusivity to thermal diffusivity.
• $\rho$ represents the Rayleigh number, which measures the thermal buoyancy force relative to viscous damping in the fluid.
• $\beta$ is a geometric factor associated with the problem.

For certain values of the parameters, $\sigma =10,$ $\rho =28,$ and $\beta =8/3$, the three-dimensional motion converges to the well-known butterfly-shaped attractor, which can be observed by solving the system numerically.

It is customary to use Lyapunov exponents to measure how quickly two initially close points in a system diverge from each other, as time progresses, indicating the rate of separation of nearby trajectories. Positive Lyapunov exponents suggest chaotic behavior and sensitivity to initial conditions, while negative or zero exponents indicate stability or convergence of trajectories over time. The existence of positive Lyapunov exponents in the Lorenz system shows the presence of chaos [1,15], and indeed, the Lorenz system can be chaotic, as shown in Figure 1.

These dynamical systems are also studied for the stability behavior, the existence of attractor points, saddle points, etc., under the time flow. As systems become increasingly complex, they require more computational power. Quantum computers, a technology of the future, promise more efficient and faster ways to solve numerical systems and speed up numerical methods for solving such differential equations.

A chaotic system exhibits complex and unpredictable behavior. Figure 2 shows two very different trajectories for two starting points that are extremely close to each other: ($1\times {10}^{-16}$, $1\times {10}^{-16}$, $1\times {10}^{-16}$) and ($1\times {10}^{-16}$, $-1\times {10}^{-16}$, $1\times {10}^{-16}$). This means that small changes in the starting point can lead to significantly different outcomes, making long-term prediction impossible. The Lorenz system is a classic example of such a system, displaying chaotic solutions for specific parameter values and initial conditions [2,16]. This system has been extensively analyzed in the literature and is a fundamental example in the study of chaotic systems [1,15].

2.1. Linearization

Consider the following approximation of the Lorenz system, (1)–(3), obtained by dropping the nonlinear terms in Equation (2), $xz$, and in Equation (3), $xy$:

$$\begin{array}{cc}\hfill \dot{x}& =\sigma (y-x),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \dot{y}& =\rho x-y,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \dot{z}& =-\beta z,\hfill \end{array}$$

(6)

where we denote $dx/dt$ by $\dot{x}$, etc.

Now, we can use the forward Euler method to discretize this system in time. If $\dot{x}=(x_{n+1}-x_n)/h$, etc., for small h, the linearized system becomes

$$\begin{array}{cc}\hfill x_{n+1}& =x_n+h\sigma (y_n-x_n),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill y_{n+1}& =y_n+h\rho x_n-h{y}_n,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill z_{n+1}& =z_n-h\beta z_n.\hfill \end{array}$$

(9)

By denoting with $w_n=(x_n,y_n,z_n)$ the three-dimensional vector, then for all n, the linearized system can be written as

$$w_{n+1}=A_L{w}_n,$$

(10)

where “L” stands for linear, and

$$A_L=\left[\begin{array}{ccc}1-h\sigma & h\sigma & 0\\ h\rho & 1-h& 0\\ 0& 0& 1-h\beta \end{array}\right],$$

(11)

is the $3\times 3$ matrix encoding the single time step, from n to $n+1$, of the linear evolution.

If the system evolves for T time steps, the evolution is described by the following matrix equation:

$$B_{L}w\equiv \left[\begin{array}{ccccc}I& \mathcal{O}& \dots & \mathcal{O}& \mathcal{O}\\ A_L& -I& \dots & \mathcal{O}& \mathcal{O}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \mathcal{O}& \mathcal{O}& \dots & -I& \mathcal{O}\\ \mathcal{O}& \mathcal{O}& \cdots & A_L& -I\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{T-1}\\ w_T\end{array}\right]=\left[\begin{array}{c}{w}_1\\ O\\ ⋮\\ O\\ O\end{array}\right]\equiv b,$$

(12)

where $B_L$ is a $3T\times 3T$ matrix; $w=(w_1,\dots ,w_T)$ is a $3T$-dimensional vector; I is the $3\times 3$ identity matrix, whereas $\mathcal{O}$ is a $3\times 3$ matrix with all 0 entries; O is a 3-dimensional vector with all 0 entries; hence, $b=(w_1,O,\dots ,O)$ is also a $3T$-dimensional vector.

Indeed, from the first line, we obtain $w_1=w_1$, but then, from the second one, we obtain $A_L{w}_1-w_2=O$, i.e., $w_2=A_L{w}_1$, that is, the matrix of Equation (10) for $n=1$, and so on until the last line, which gives $A_L{w}_{T-1}-w_T=O$, i.e., $w_T=A_L{w}_{T-1}$, that is, the matrix of Equation (10) for $n=T-1$.

The evolution of the system is given by $w=(w_1,w_2,\dots ,w_T)$, which can be obtained by inverting the block diagonal matrix $B_L$ and multiplying it by b, i.e., $w=B_L^{-1}b$.

We can also use HHL to compute the solution of system (12). This linearization on a collection of qubits along with HHL is used in [12] for solving atmospheric dynamics by using quantum approaches. The advantage offered is the exponential speedup due to HHL subject to limitations posed by the process of initial state preparation, solution readout, and noisy quantum bits constraints. A similar linearization is used in [14] along with a Variational Quantum Eigensolver to study the stability of a simplified model of the atmospheric system (Lorenz system) at equilibrium points.

2.2. Inclusion of Nonlinear Terms

When the nonlinear terms are disregarded, the resulting system yields to a simpler and more tractable analysis. However, the linear system has important limitations:

• It is unable to capture global behavior, especially in regions far from equilibrium.
• Chaotic behavior, bifurcations, and other complex dynamics are completely overlooked by linear analysis. In particular, it does not encompass chaos, which is a defining feature of the Lorenz system.
• The effectiveness of the linear approximation can vary significantly with changes in parameters ($\sigma$, $\rho$, $\beta$). Some dynamics observable in one set of parameter values may be absent in others.

To overcome these challenges, we propose here a new method incorporating nonlinear terms. The discretized version of the nonlinear Lorenz system in (1)–(3), obtained through the very same steps as for the linear case, is simply

$$\begin{array}{ccc}\hfill x_{n+1}& =& x_n+h\sigma (y_n-x_n),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill y_{n+1}& =& y_n+h[x_n(\rho -z_n)-y_n],\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill z_{n+1}& =& z_n+h(x_n{y}_n-\beta z_n).\hfill \end{array}$$

(15)

The value h is the timestep that controls the resolution of the simulation, e.g., $h={10}^{-5}$. It is important to handle the timestep h carefully to ensure stability and accuracy.

We rewrite this system as

$$\begin{array}{ccc}\hfill x_{n+1}-h\sigma y_n-x_n(1-h\sigma )& =& 0\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill y_{n+1}-h{x}_n(\rho -z_n)-y_n(1-h)& =& 0\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill z_{n+1}-h{x}_n{y}_n-z_n(1-\beta h)& =& 0.\hfill \end{array}$$

(18)

In matrix form, this is

$$A_{NL}W=b_{NL}$$

(19)

where “$NL$” stands for nonlinear; $A_{NL}$ is an $8\times 8$ matrix, which we shall explicitly write in a moment; and W and $b_{NL}$ are 8-dimensional vectors given by

$$W=(x_n,y_n,z_n,x_{n+1},y_{n+1},z_{n+1},x_n{z}_n,x_n{y}_n)$$

(20)

and

$$b_{NL}=(x_n,y_n,z_n,0,0,0,x_n{z}_n,x_n{y}_n).$$

(21)

The vector W contains the solutions.

It is important to notice that the matrix in Equation (19) refers to a single time step; hence, it is the nonlinear generalization of Equation (10), and not that of Equation (12), that refers to multiple time steps. With the latter, it shares the higher dimensionality, compared with the linear system of Equation (10) (eight dimensions vs. three dimensions), but that is due to the inclusion of the nonlinear terms into a single time step.

Thus, system (16)–(18), in its matrix form (19), can be written as

$$\begin{array}{c}\hfill A_{NL}W\equiv \left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ -(1-h\sigma )& -h\sigma & 0& 1& 0& 0& 0& 0\\ -h\rho & -(1-h)& 0& 0& 1& 0& h& 0\\ 0& 0& -(1-\beta h)& 0& 0& 1& 0& -h\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_n\\ y_n\\ z_n\\ x_{n+1}\\ y_{n+1}\\ z_{n+1}\\ x_n{z}_n\\ x_n{y}_n\end{array}\right]=\left[\begin{array}{c}{x}_n\\ y_n\\ z_n\\ 0\\ 0\\ 0\\ x_n{z}_n\\ x_n{y}_n\end{array}\right]\equiv b_{NL}.\end{array}$$

(22)

The solution $x_{n+1},y_{n+1},z_{n+1}$ is obtained by inverting the matrix $A_{NL}$ and multiplying it by $b_{NL}$.

This inversion can be performed classically or by using HHL [11], as in [12]. HHL assumes a Hermitian matrix; therefore, one has to solve for ${\mathcal{A}}_{NL}{\mathcal{W}}_{NL}=B_{NL}$, where ${\mathcal{A}}_{NL}=\left[\begin{array}{cc}{\mathcal{O}}_8& A_{NL}\\ A_{NL}^{†}& {\mathcal{O}}_8\end{array}\right]$, ${\mathcal{W}}_{NL}=({\mathbf{O}}_8,W)$, and $B_{NL}=(b_{NL},{\mathbf{O}}_8)$, with ${\mathcal{O}}_8$ being an $8\times 8$ matrix with all 0 entries and ${\mathbf{O}}_8$ being an eight-dimensional column vector of zeros. We need 4 qubits to represent the input $b_{NL}$, 1 qubit for controlled rotation, and, say, 4 qubits for the clock register (needed for the simulation of the Hamiltonian). So, a total of 9 qubits are needed if HHL is used to solve $A_{NL}W=b_{NL}$. The precision increases with the increase in the number of clock qubits. The condition number for ${\mathcal{A}}_{NL}$ is 3.03 for an h, step size, of $0.01$.

Although the HHL and VQLS methods are both prominent algorithms for solving linear systems, they present distinct advantages and limitations. The HHL algorithm is known for its theoretical exponential speedup over classical methods for certain classes of linear systems, particularly those matrices A that are sparse and well-conditioned [17,18]. However, HHL is highly susceptible to quantum gate errors and qubit decoherence, which significantly impairs its performance on current quantum hardware [19,20]. Additionally, the simulation of HHL is computationally intensive.

In contrast, the VQLS method employs simpler and often shorter quantum circuits, making it more feasible on available noisy intermediate-scale quantum (NISQ) hardware [21,22]. While the VQLS may not achieve the same exponential speedup as HHL, its variational approach facilitates noise mitigation, resulting in greater robustness and reliability in noisy quantum environments [23,24]. Additionally, the VQLS requires fewer quantum resources by combining quantum subroutines with classical optimization, making it particularly suited for NISQ devices, where quantum coherence time is limited [25,26].

Here, we should emphasize again that while this work and the work [6,12] use the forward Euler method to discretize the Lorenz system, there are two main differences. First, in this work, we simulate only a single time step, whereas [12] simulates several time steps (T) in a single iteration. The second big difference is the use of a collection of qubits in [6,12] to give a mean-field approximation. Since we simulate only a single time step, we do not need a collection of qubits to avoid measurement errors that compound with time steps.

On the other hand, Berry [7] has an efficient quantum procedure for linear differential equations, while in the present work, we consider a nonlinear system, so the approach here is not directly comparable with the results there.

System (16)–(18) has nonlinear terms, specifically $xz$ and $xy$, which distinguish it significantly from the linear system presented in [14]. However, the present system is not as complex as the nonlinear system presented in [11].

To show the efficiency of the proposed method, we compare the results of our proposed method and the Euler method, which is a nonlinear approach, to solve the Lorenz system. We use the same number of iterations, step size, and initial point. We plot this compression in Figure 3.

Encouragingly, the solutions obtained from our proposed linear system approach are similar to the those from nonlinear Euler method. This result demonstrates that although our method is based on solving a linear system, it effectively captures the dynamics of the nonlinear Lorenz system. To solve (22), we begin from an initial point and use iterative methods to find the optimal solution. While numerous classical methods are available [27], we will use a variational quantum algorithm called the VQLS.

3. Variational Quantum Linear Solver

3.1. Overview

The VQE is a sophisticated hybrid quantum–classical algorithm designed to determine the ground state energy of quantum systems, such as molecular structures. The VQE operates by utilizing a parameterized quantum circuit to prepare an approximate quantum state, which is then measured to estimate the expectation value of the Hamiltonian, representing the system’s energy. This energy estimation is subsequently fed into a classical optimization process that iteratively adjusts the circuit’s parameters to minimize the energy, ultimately converging on the ground state [28]. The VQLS [21] is a hybrid quantum algorithm for solving linear systems of equations by using the variational principle, which minimizes the expectation value of the Hamiltonian of the system, with respect to a parameterized quantum circuit. Given A and b, the VQLS aims to prepare variationally a state $|w\rangle$ such that $A|w\rangle \sim |b\rangle .$ One of the main appeals of the VQLS is that it can be implemented on near-term NISQ computers. Experimental studies on a limited set of test instances on Rigetti machines seem to offer evidence that the VQLS scales linearly in the condition number $\kappa$ and $log{\displaystyle \frac{1}{ϵ}}$ where $ϵ$ is the desired precision [21].

The VQLS is a parameterized circuit. To find a good set of values for the parameters, classical optimization is used, which is computationally challenging. Estimating the values with a precision of $\pm 1/poly\left(n\right)$ is DQC1-hard, where n represents the number of qubits. It is believed that classical algorithms cannot efficiently find precise values of the parameters because the efficient simulation of DQC1 would lead to the collapse of the polynomial hierarchy to the second level. Consequently, there are doubts about the efficient classical simulation of the VQLS.

In order to effectively use the VQLS algorithm, the input matrix A must satisfy certain requirements. A should be representable as a linear combination of unitaries, similar to how the Hamiltonian is represented in the Variational Quantum Eigensolver as a linear combination of Pauli operators. The method provided by [21] is based on Szegedy walks and efficiently decomposes a sparse matrix into a linear combination of unitaries. A should also be sparse and well-conditioned, a finite $\kappa$. The norm of A should be bounded, with $\left|\right|A\left|\right|\le 1$. Lastly, the unitaries in the decomposition must be efficiently implementable. All these assumptions are satisfied by the matrix in (22).

We assume that A can be expressed as a linear combination of unitary operators, such that

$$A=\sum _i{c}_i{A}_i,$$

(23)

where $A_i$ are the unitaries and $c_i$ are complex coefficients. This representation effectively models a system Hamiltonian. Typically, the decomposition involves a linear combination of tensor products of the identity and Pauli matrices. These gates are widely used due to their well-known properties and ease of implementation. The matrix representation of the gates are below.

$$\begin{array}{c}\hfill I=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),X=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),Y=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),Z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).\end{array}$$

(24)

We use a recently proposed algorithm to decompose a square real symmetric matrix of any size into a tensor product of Pauli spin matrices for all the application matrices discussed. This algorithm, which is detailed in [29], is available in Pennylane and was utilized to generate decompositions for stiffness matrices of general sizes commonly encountered in discrete finite-element methods. We do not use the decomposition algorithm from [21].

3.2. Cost Function

Two types of cost functions have been introduced for the VQLS method: local cost functions and global cost functions. We describe the cost functions and highlight their features.

The residual-based cost function is given by

$$C\left(\theta \right)=\underset{\theta }{min}\parallel A|\psi \left(\theta \right)\rangle -{|b\rangle \parallel }^2=\underset{\theta }{min}{\langle \left(A\right|\psi \left(\theta \right)\rangle -|b\rangle )}^{†}\left(A\right|\psi \left(\theta \right)\rangle -|b\rangle )\rangle .$$

(25)

This cost function can be viewed as the expectation value of an effective Hamiltonian, which is defined as follows [21]:

$$H_G=A^{†}\left(I-|b\rangle \langle b|\right)A.$$

(26)

Therefore, we can write the cost function associated with Hamiltonian $H_G$ as

$$C_G\left(\theta \right)=\langle \psi \left(\theta \right)|H_G|\psi \left(\theta \right)\rangle .$$

(27)

3.3. Ansatz

In variational algorithms, an Ansatz refers to an assumed initial form for the quantum state. This Ansatz is typically represented by a parameterized quantum circuit, which is used to prepare a trial state that can be optimized. The Strongly Entangling Layer (SEL) Ansatz is a type of parameterized quantum circuit used in variational quantum algorithms, such as the VQE and Quantum Approximate Optimization Algorithm (QAOA). This Ansatz is designed to introduce a high degree of entanglement between qubits while maintaining a relatively simple and regular structure, making it a popular choice for many quantum applications. Let us describe the structure of the SEL Ansatz.

Each SEL Ansatz consists of layered circuits composed of rotation gates and entanglement operations:

• Parameterized rotation gate: Each qubit undergoes a parameterized rotation, typically represented as $R(\alpha ,\beta ,\gamma )$. This gate can be a combination of rotations around different axes, such as $R_X$, $R_Y$, and $R_Z$.
• Entangling operations: Following the rotation gates, to entangle the qubits, a series of controlled-NOT (CNOT) gates are applied. The pattern of these CNOT gates may vary, but they generally ensure that every qubit is entangled with at least one other qubit in the layer.
• Layered structure: To enhance the expressiveness of the Ansatz, multiple layers of the above combination are stacked. Each layer applies a new set of parameterized rotation gates $R(\alpha ,\beta ,\gamma )$,

  $$\begin{array}{c}\hfill R(\alpha ,\beta ,\gamma )=R_Z\left(\gamma \right)R_Y\left(\beta \right)R_Z\left(\alpha \right)=\left[\begin{array}{cc}{e}^{-i(\alpha +\gamma )/2}cos(\beta /2)& -e^{i(\alpha -\gamma )/2}sin(\beta /2)\\ e^{-i(\alpha -\gamma )/2}sin(\beta /2)& e^{i(\alpha +\gamma )/2}cos(\beta /2)\end{array}\right],\end{array}$$

  (28)

Followed by entangling CNOT gates, as shown in Figure 4.

4. Algorithm

To ease the notation, since the focus will be fully on nonlinear system (22), from here on, we indicate with A the matrix $A_{NL}$, with w the vector W, and with b the vector $b_{NL}$; see Equations (19)–(22).

Algorithm 1 outlines the VQLS method for solving a system of linear equations. The algorithm starts with the following inputs: a specified matrix A, a vector b, the number of layers, the maximum number of iterations, the convergence tolerance, and a fixed step size for the optimization process.

Algorithm 1 Variational quantum linear solver [21]

Require: Matrix A, vector b, number of layers $num\_layers$, maximum iterations $max\_iterations$, convergence tolerance $conv\_tol$, step size $stepsize$

Ensure: Optimized parameters

1: $A^{†}\leftarrow$ Hermitian conjugate of A

2: $b\_norm\leftarrow b/\parallel b\parallel$

3: $P_b\leftarrow b\_norm\otimes b\_nor{m}^T$

4: $I\leftarrow$ Identity matrix of size A

5: $H_G\leftarrow A^{†}·(I-P_b)·A$

6: Define cost function $C_G$

7: $num\_qubits\leftarrow \lceil {log}_2\left(\mathrm{size}\mathrm{of}A\right)\rceil$

8: Define quantum device with $num\_qubits$

9: Define Ansatz with $num\_layers$

10: Initialize optimizer with $stepsize$

11: Initialize random parameters for Ansatz

12: for $it\leftarrow 1$ to $max\_iterations$ do

13:       Compute gradient

14:       Update parameters and compute the cost function

15:       Check stop condition

16: end for

17: return optimized parameters

18: Extract the solution

In lines 1–6 of the algorithm, the cost function is defined based on the matrix A and the vector b. After defining $H_G$, the cost function is evaluated by using the expectation value of $H_G$. The number of qubits and the Ansatz are specified. Initially, the parameter $\theta$ is assigned a random value. The rest of the algorithm is the classic gradient descent loop. The algorithm then computes the cost function and gradient within the for loop and updates the parameters $\theta$. This process repeats until the stopping condition is met. Once the optimal value for the parameter $\theta$ is found, the solution to the system $Aw=b$ is determined. It is important to note that the solution must satisfy $Aw-b=0$ to be considered valid.

5. Numerical Results

We present the results of implementing quantum Algorithm 1 and classical approaches to solving nonlinear Lorenz system (22). We use the following values of the initial parameters:

• Number of layers: The number of layers is a crucial parameter for achieving optimal parameter values. We ran the algorithm with the same initialization and plotted these behaviors of the cost functions as shown in Figure 5 as a function of the number of layers. Three or more layers appear to be necessary to compute a good set of values for the parameters.
  Figure 5. The expectation value as a function of the number of layers.

  Figure 5. The expectation value as a function of the number of layers.
• Stop condition: We used two different criteria for the stop condition, i.e., the maximum number of iterations and convergence tolerance, given by

  $$\text{Max-iter}=200,ϵ=1\times {10}^{-8},$$

  respectively.
• Initialization: The initial value of the parameter $\theta$ was chosen randomly.
• Cost function: We considered the following cost function:

  $$C_G\left(\theta \right)=\langle \psi \left(\theta \right)|H_G|\psi \left(\theta \right)\rangle$$
• Number of qubits: Error-free qubits represent an important quantum resource. In our implementation, three qubits are needed, which is a very low resource requirement compared with algorithms based on HHL in which an extra register is needed to store the eigenvalue computed during the phase inversion process. For greater precision, a larger-sized quantum register is needed to store the eigenvalues. In general, the VQLS algorithm needs $logn$ qubits, where n is the size of the matrix A.

5.1. Solution

Quantum states have unit norms. Thus, to find the exact solution that satisfies the system of equations $Aw=b$, where b is not of unit length, we need to scale the state vector. Once the optimal parameters are obtained by using the classical optimization method, we execute the Ansatz to obtain the final state. Now, this state can only be measured, and we make an assumption that we can reconstruct w from $|w\rangle$ efficiently. Next, we seek a coefficient C such that the solution of the VQLS can be scaled by C to match the exact solution $w=A^{-1}b$. To determine the value of C, we divide $\parallel b\parallel$ by $\parallel A\tilde{w}\parallel$. By multiplying this factor by the solution, we obtain the exact solution.

5.2. Handling the Sign of the Solution

We know that if w is an eigenvector of the matrix A, then $-w$ is also an eigenvector of A. To obtain the correct solution, we check the distance between $Aw$ and b for the i-th element in the solution. If this distance is less than the absolute value of $Aw\left(i\right)$, we do not need to multiply by $-1$. Otherwise, we must multiply by $-1$ to ensure the correct solution.

5.3. Condition Number $\kappa$

A well-known measure of how ill-conditioned a non-singular matrix A is is the condition number $\kappa \left(A\right)=\left|\right|A^{-1}\left|\right|\ast \left|\right|A\left|\right|$, where $\left|\right|\circ \left|\right|$ is any norm. When 2-norm is used, the condition number is the ratio of the largest to the smallest eigenvalue. The condition number of the input matrix A significantly affects the performance of the VQLS algorithm; the running time is linear in $\kappa$. Specifically, for the matrix A defined by (22), the condition number depends largely on the step size parameter h. As h gets closer to 1, the condition number of A increases exponentially in general. On the other hand, the lowest condition number is achieved as h tends to zero. As mentioned in [21], the VQLS algorithm performs optimally when the condition number is small. In our study, we use $h=0.01$ to achieve better performance. Figure 6 shows the relationship between the condition number of the matrix A and the value of h. Even for large values of the step size ($\simeq 0.1$), the condition number is bounded by 70.

5.4. Starting Point

To solve system (22) by using the VQLS, we begin by initializing the algorithm with the following values:

• $b=(1,-2,4,0,0,0,4,-2)$, $h=5\times {10}^{-3}$, $\sigma =10$, $\rho =28$, $\beta ={\displaystyle \frac{8}{3}}$, and $T=2000$.
• $b=({10}^{(-16)},-{10}^{(-16)},{10}^{(-16)},0,0,0,0,0)$, $h={10}^{-3}$, $\sigma =10$, $\rho =13.92655742$, $\beta ={\displaystyle \frac{8}{3}}$, and $T=10,000$.

5.5. Error Analysis

To compare the trajectories calculated by using the classical and quantum algorithms, as outlined in Section 2.2, we can analyze them point by point. For each time step, we calculate the difference in the coordinate positions obtained from the two methods and normalize them by using the formula provided next. Let $w_n^c=(x_n^c,y_n^c,z_n^c,)$ and $w_n^q=(x_n^q,y_n^q,z_n^q,)$ represent the points computed in the classical and quantum algorithms, respectively, at iteration $I=n$. The relative error is defined as

$${\displaystyle \frac{|x_n^c-x_n^q|+|y_n^c-y_n^q|+|z_n^c-z_n^q|}{1+|x_n^c|+|y_n^c|+|z_n^c|}}.$$

(29)

The relative error for the two algorithms after 500 iterations is illustrated in Figure 7. This shows that the trajectories calculated by the quantum method proposed in this paper closely match the trajectory computed by using a classical computer.

However, analyzing the relative error alone does not provide insights into the impact of step size on error. Therefore, we conducted a more comprehensive analysis by using Richardson’s extrapolation method.

For a deeper analysis of errors, we consider the error in computing the gradients $\dot{x},\dot{y},$ and $\dot{z}$. This error serves as a good proxy for the error in computing the trajectories, as the only error that occurs in each iterative step is an error in the gradient computation. We make the strong assumption that errors arising from quantum uncertainty are negligible.

Let $\dot{x_h}$ be the numerically computed gradient as a function of the step size h, where $h>0$ represents a positive number close to zero. The numerical gradient $\dot{x_h}$ can be expressed as the actual gradient $\dot{x}$ plus some error $E\left(h\right)$, which is a function of h:

$${\dot{x}}_h=\dot{x}+E\left(h\right),$$

(30)

where $E\left(h\right)$ is the error associated with the forward difference method. The error term $E\left(h\right)$ for the forward difference method is

$$\begin{array}{c}\hfill E\left(h\right)=c·h+O\left(h^2\right),\end{array}$$

(31)

where c is a constant. This relationship is true for all step sizes h. Therefore, the numerically computed derivatives for step sizes h and $2h$ can be written as

$$\begin{array}{c}\hfill {\dot{x}}_h=\dot{x}+ch+O\left(h^2\right)\end{array}$$

(32)

$$\begin{array}{c}\hfill {\dot{x}}_{2h}=\dot{x}+c\left(2h\right)+O\left(h^2\right).\end{array}$$

(33)

We obtain the error for the step size h by subtracting the second equation from the first one:

$${\dot{x}}_{2h}-{\dot{x}}_h=ch+O\left(h^2\right)=E\left(h\right).$$

(34)

This intuitive calculation can be formalized, as shown in [30]. At each time step n, we use the numerical gradients computed at step sizes $h,2h$ to bound the error in the computation. The initial point for both the step sizes is point $(x_{n-1},y_{n-1},z_{n-1})$. Similarly, we can obtain the error in computation of $\dot{y}$ and $\dot{z}$. The numerical derivatives are computed in two ways: classical and using a quantum algorithm. We compute the error for both methods.

5.6. Error for VQLS Method

Here, we briefly review some key concepts related to errors in the VQLS method [21]. For the VQLS algorithm, the deviation between observable expectation values for the approximate solution $\left|x\right({\alpha }_{\mathrm{opt}})\rangle$ and the true solution $|x_0\rangle$ can be upper-bounded based on the value of the cost function. Specifically, the error tolerance $ϵ$ can be set before running the algorithm, where $ϵ$ is defined as the trace distance between the exact and approximate solutions:

$$ϵ={\displaystyle \frac{1}{2}}\mathrm{Tr}\left(∥|x_0\rangle \langle x_0|-|x\left({\alpha }_{\mathrm{opt}}\right)\rangle \langle x\left({\alpha }_{\mathrm{opt}}\right)|∥\right).$$

(35)

Moreover, it has been proven that the cost function satisfies the relation [21]

$$C_G\ge {\displaystyle \frac{{ϵ}^2}{{\kappa }^2}},$$

(36)

where $\kappa$ is the condition number of the input matrix, and $C_G$ is the cost function defined by (27). By using the above formulas, we can observe that the difference between the exact solution and the quantum solution is related to ${\kappa }^2$ multiplied by the value of the objective function. This shows that the error between the quantum solution and the exact solution is inherently connected to the condition number of the matrix and the chosen cost function, allowing us to control and quantify the error before running the algorithm.

Several noise mitigation strategies are particularly well suited to enhance the robustness of the VQLS on NISQ devices due to the method’s variational structure, and we can apply these methods. Techniques such as Zero-Noise Extrapolation (ZNE) and Measurement Error Mitigation (MEM) allow the VQLS to reduce the impact of noise without increasing the quantum circuit depth, making them practical for the shallow circuits typically used in variational algorithms [31,32]. Additionally, the Variational Error Suppression (VES) technique can adaptively adjust the parameters within the VQLS to minimize noise effects, enhancing solution accuracy even on noisy devices [22,25]. Together, these approaches demonstrate VQLS’s practical suitability for today’s quantum hardware while mitigating the effects of common noise sources.

5.7. Results and Discussion

The VQLS circuit has a short depth (20). Short-depth circuits can be implemented in NISQ computers because these circuits require fewer quantum gates and operations to execute. NISQ computers are currently limited by noise and errors, which makes long computations and complex circuits prone to errors. Short-depth circuits require fewer operations, reducing the likelihood of errors and increasing the chances of successfully executing the computation on current NISQ hardware. Short-depth circuits also offer better error mitigation [31].

One of the critical resources used by any quantum algorithm is the number of qubits, which significantly impacts the simulation time. The algorithm proposed here to solve the Lorenz system using the VQLS method requires only three fault-tolerant logical qubits. However, since the simulations were performed by using a classical computer, the execution time is slightly longer than that of classical methods. However, it should be noted that this simulation time is insignificant compared with the simulation of the HHL method. We coded the VQLS and a classical algorithm.

An important observation is that although we solved a linear system of equations, the results resemble those of a nonlinear system. The new system has nonlinear product terms, so it simulates the nonlinear system of differential equations well in low dimensions.

For two different starting points, the trajectories computed by using the classical and the VQLS methods are shown in Figure 8. Both approaches compute similar trajectories. The chaotic system appears to have more error compared with the attractor.

We use Equation (34) to estimate the error for algorithms (both classical and quantum) as a function of the step size. We use the error measure to quantify the error of the two approaches individually. We also use the error to perform a comparative analysis. The VQLS algorithm takes considerable time to simulate; therefore, only a limited comparative study is conducted.

First, we examine the computation shown on the left side of Figure 1. The starting point is $(1,-2,4)$. Figure 9a shows the error in the computation of $\dot{x},\dot{y},$ and $\dot{z}$ in each iteration and the total error (defined as the sum of individual errors). It is evident that equation system (22) solved by using the classical approach has good accuracy. The average total error is the same for a step size of ${10}^{-3}$. The error decreases as the number of iterations increases.

Next, we study the classical error in the computation of the blue curve shown in Figure 2. Figure 9b plots the individual and the total error as a function of the iteration for a step size of ${10}^{-3}$. For the first 200 iterations, the total error is less than ${10}^{-16}$. Initially, the error is insignificant, but for the last half of the computation, the error is significant, close to 100 times greater than the step size. The average error is still comparable to the step size.

The total error in the two figures is of the order of the step size, so we infer that the system given by (22) is a good model.

Next, we compare the error of the quantum method with the classical method. Since the VQLS take a long time to simulate, we limit the experiments to two starting points and the number of iterations to 200 in the preliminary results reported here. The individual errors and the total error in the computation of $\dot{x},\dot{y},$ and $\dot{z}$ for the trajectories shown in Figure 8a,b are shown in Figure 10a,b, respectively. It appears that the error of the hybrid quantum method is comparable.

5.8. Limitations

The VQLS-based method is capable of handling nonlinearities. However, it does not offer the same exponential advantage as HHL when dealing with a “giant” linear transformation. The effectiveness of the VQLS method depends upon parameter values computed by using a classical optimization approach, which are also sensitive to the initial points. The method relies on strong assumptions, such as the requirement for the exact preparation of initial states, noise-free evolution, and the ability to recover the final answer with the desired precision. The impact of quantum uncertainty on trajectories is not well understood; while we have shown that the trajectories computed by classical and quantum algorithms are “close”, it is conceivable that quantum uncertainty can lead to completely different trajectories from the classical ones, which requires further examination. The VQLS is a hybrid algorithm, and simulating it on a classical computer is more time-consuming than computing the inverse of a small matrix. At present, there is no computational advantage for nonlinear difference equations involving few variables.

Figure 10. Classical error vs. quantum error. (a) Errors for the trajectories computed in Figure 8a. (b) The error for the trajectory computed in Figure 8b.

Figure 10. Classical error vs. quantum error. (a) Errors for the trajectories computed in Figure 8a. (b) The error for the trajectory computed in Figure 8b.

It is worth noting that the VQLS algorithm proposed for solving the Lorenz system aligns with existing implementations of quantum linear solvers for near-term quantum devices [21]. In that study, a 10-qubit implementation of the VQLS was performed on Rigetti’s Aspen-4 quantum computer to solve a Quantum Linear System Problem (QLSP) with a 1024 × 1024 matrix. This experiment utilized a specifically tailored Ansatz incorporating $R_y\left({\alpha }_i\right)$ gates and computed the cost function by expanding the Hamiltonian in terms of Pauli operators to meet hardware constraints. The findings indicate that the cost function was minimized effectively on quantum hardware, closely matching simulated results and verifying that the solution to the linear system was obtained. This research supports the feasibility of implementing our proposed method on current quantum devices.

5.9. Matrix Formulation for High-Dimensional Systems

To extend the proposed method to solve larger nonlinear systems of differential equations, we first need to ensure that the matrix $A_{NL}$ defined by (19) is compatible with quantum algorithms, which often require matrix dimensions to be powers of 2. If the original matrix $A_{NL}$ does not meet this requirement, we can pad it by adding an identity matrix or zero elements to reach the nearest $2^n\times 2^n$ dimension. This allows the system to be compatible with the VQLS and enables efficient processing on quantum hardware.

Consider a general nonlinear system of ordinary differential equations given by

$$\begin{array}{c}\hfill {\displaystyle \frac{d\mathbf{u}}{dt}}=\mathbf{f}\left(\mathbf{u}\left(t\right)\right),\end{array}$$

(37)

where $\mathbf{u}\left(t\right)={[u_1\left(t\right),u_2\left(t\right),\dots ,u_N\left(t\right)]}^{\top }$ is an N-dimensional vector of state variables and $\mathbf{f}\left(\mathbf{u}\right)={[f_1\left(\mathbf{u}\right),f_2\left(\mathbf{u}\right),\dots ,f_N\left(\mathbf{u}\right)]}^{\top }$ represents the nonlinear terms.

By using a timestep h, the discretized form at each time step n can be written as

$$\begin{array}{c}\hfill {\mathbf{u}}_{n+1}={\mathbf{u}}_n+h\mathbf{f}\left({\mathbf{u}}_n\right),\end{array}$$

(38)

which we rearrange as

$$\begin{array}{c}\hfill {\mathbf{u}}_{n+1}-{\mathbf{u}}_n-h\mathbf{f}\left({\mathbf{u}}_n\right)=0.\end{array}$$

(39)

Now, we can rewrite (39) into matrix form as

$$A_{NL}W=b_{NL},$$

(40)

where $A_{NL}$ is an $(N+M)\times (N+M)$ matrix (or the nearest $2^n\times 2^n$ size if padding is applied), W is an $(N+M)$-dimensional vector containing state variables at the current and next time steps, and $b_{NL}$ denotes an $(N+M)$ -dimensional vector with known values. Now, we can apply the VQLS to solve the linear system given by (40).

6. Conclusions

This is an exploratory study on the use of quantum algorithms for studying complex systems. We show how a variational quantum algorithm can be used to solve the Lorenz system. We perform a comparative error analysis of the quantum method and the classical method. The quantum method is found to be reliable in the simulations. The method has limitations, such as dependence on classical optimization, the readout–read-in capabilities of quantum systems, and noise in quantum computation.