Euler Lagrange Equation Solver

[Aiko Bartels]

Thankyou for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Recently, the advent of quantum computers and the invention of new quantum algorithms have provided a novel paradigm for solving PDEs. A cornerstone of many of these quantum algorithms is the seminal Harrow-Hassidim-Lloyd (HHL) algorithm1 for solving linear systems, which can be utilized to solve PDEs by discretizing the PDE and mapping it to a system of linear equations. Compared to classical algorithms, the HHL algorithm can be shown to exhibit an exponential speedup. Unfortunately, attractive as it may sound, the HHL algorithm works only in an idealized setting, and a list of caveats must be addressed before it can be used to realize a quantum advantage2. Moreover, implementing HHL and many other quantum algorithms would require the use of a fault-tolerant quantum computer, which may not be available in the near future3. Instead, the machines we have today are imperfect, noisy intermediate-scale quantum (NISQ) devices4 with both coherent and incoherent errors limiting practical circuit depths.

Over the last few years, variational quantum algorithms (VQAs) have emerged as a leading strategy to realize a quantum advantage on NISQ devices. Specifically, VQAs employ shallow circuit depths to optimize a cost function, expressed in terms of an Ansatz with tunable parameters, through iterative evaluations of expectation values5. Applications of VQAs include the variational quantum eigensolver (VQE) for finding the ground or excited states of a system Hamiltonian6,7,8, the quantum approximate optimization algorithm (QAOA) for solving combinatorial optimization problems9, and solvers for linear10,11,12 and non-linear13 systems of equations.

Consider the second-order homogeneous evolution equation defined on the set \(\Omega \times J\), where \(\Omega \subset \mathbb R^d\) denotes a d-dimensional bounded spatial domain and \(J = [0,T]\), where \(T>0\) denotes a bounded temporal domain, as

where \(\mathcal I\) is the identity matrix of the same size, \(u^k=\left[ u_i j^k\right] _0 \le i \le n_x, 0 \le j \le n_y\) and \(f^k=\left[ f_i j^k\right] _0 \le i \le n_x, 0 \le j \le n_y\).

where \(f^k+1/2=(f^k+1+f^k )/2\). However, the CN method may introduce spurious oscillations to the numerical solution for non-smooth data unless the algorithm parameters satisfy the maximum principle25.

Here, we explore a variational quantum approach towards the solution of the evolution equation (1). In addition to potential quantum speedup, a variational quantum algorithm could also benefit from data compression, where a matrix of dimension N can be expressed by a quantum system with only \(\log _2 N\) qubits, where N is the size of the problem. Consider the Poisson equation, which is a time-independent form of Eq. (1), expressed as

Since expectation values of the identity operator are equal to 1, i.e. \(\langle \phi I^\otimes n \phi \rangle = \langle \phi ' I^\otimes n \phi ' \rangle = 1\), evaluating the expectation value of the operator \(A_x, \beta \) requires only the evaluation of expectation values of the simple Hamiltonians \(H_1-4\) (\(H_1-3\) for Dirichlet boundary condition). The required number of quantum circuits is therefore limited to a constant \(O(n^0)\)18. Similar decomposition expressions apply to problems of higher dimensions, including \(A_y,\beta \) in the y direction19.

Once the matrix A is decomposed, a parameterized quantum state \( \psi (\theta ) \rangle\) is prepared using an Ansatz represented by a sequence of quantum gates \(U(\theta )\) parameterized by \(\theta\) applied to a basis state \( 0 \rangle ^\otimes n\), such that \( \psi (\theta ) \rangle = U(\theta ) 0 \rangle ^\otimes n\). Here, we use a hardware-efficient Ansatz consisting of multiple layers of \(R_Y\) gates across n qubits entangled by controlled-X gates (see Fig. 1). For the source term f in (10), a quantum state \( b \rangle\) is prepared by encoding a real vector with the unitary \(U_b\), such that \( b \rangle = U_b 0 \rangle ^\otimes n\). Depending on the actual input, conventional amplitude encoding methods26,27 may introduce a global phase that must be corrected by its complex argument for computing in the real space.

In this study, we propose to solve the evolution equation (1) through successive time-stepping of the quasi-steady Poisson equation using a variational quantum algorithm. Using a parameter set \(\theta ^k\) obtained at time-step k, we encode a normalized source state \( \hatb^k \rangle := b \rangle /\sqrt\langle b\) from \( b(\theta ^k) \rangle\) and seek an implicit solution to

where \(\theta ^k+1=\theta _\mathrm opt(t^k+1 )\) is the parameter set and \(r^k+1=r_\mathrm opt(\theta ^k+1)\) is the norm at next time-step \(k+1\). This process is then iterated in time up to \(n_t\) number of time-steps as desired (see Algorithm 1).

The cost function \(E(\theta ^k)\) may be evaluated on a quantum computer by computing each of the inner products in the expression in Eq. (17) separately. Using the decomposition provided in Eq. (16), these inner products may be expressed in terms of expectation values \(\langle \varphi O_i \varphi \rangle\) for preparable states \( \varphi \rangle\) and simple Hermitian operators \(O_i\). Each expectation value \(\langle \varphi O_i \varphi \rangle\) is evaluated on a quantum computer by preparing the state \( \varphi \rangle\) using the quantum circuits described above and then measuring the operator \(O_i\) in the state \( \varphi \rangle\)18.

In this study, the variational quantum algorithm is implemented in Pennylane (Xanadu)28 using a statevector simulator with the Qulacs29 plugin as a backend for quantum simulations, and the L-BFGS-B optimizer for parametric updates. Amplitude encoding is carried out via the standard Mortonnen state preparation template30 with custom global phase correction. For hardware emulation via the QASM simulator (Qiskit), we refer the reader to the excellent cost-sampling analysis of Sato et al. 18.

To solve Eq. (22), the variational quantum evolution algorithm (Algorithm 1) can be employed with a suitable time-stepping scheme (7). For the implicit Euler (IE) method (8), the matrix A and source state \( b(\theta ^k) \rangle\) can be decomposed into

For the Crank-Nicolson (CN) method (9), it follows that the \(\mathcal Au^k\) term, which carries a small but non-trivial evaluation cost, can be eliminated using the source state of the previous time-step \(k-1\), leading to

where \(2\left \baru_D^k+1/2\right\rangle :=\left( \left u_D^k+1\right\rangle +\left u_D^k\right\rangle \right)\). Here, the presence of a \(k-1\) term in \(\left b^k-1\right\rangle\) is not unexpected due to temporal finite differencing at second-order accuracy.

For a space-time domain \(\Omega \times J\in [0,1] \times [0,1]\), let the number of time-steps be \(n_t=20\) and the spatial intervals be \(n_x=2^n+1\), where n is the number of qubits, and \(\delta _x=1\) is the diffusion parameter. We employ the Dirichlet boundary condition with boundary values \((g_L,g_R )=(1,0)\) and initial values \(u_0= \mathbf 0\). With initial random parameters \(\theta ^0 \in [0,2 \pi ]\), we run a limited-memory Broyden-Fletcher-Goldfarb-Shanno boxed (L-BFGS-B) optimizer31,32,33,34 to optimize \(\theta\) with absolute and gradient tolerances set at \(10^-8\).

where \(\left \hatu^k\right\rangle :=\left u^k\right\rangle / \sqrt\left\langle u^k \mid u^k\right\rangle \) is the normalized classical solution vector at time k. The trace errors of solutions shown in Fig. 2a are 0.0008 and 0.0025 for n- l sets of 3-3 and 4-4 respectively.

Figure 2b shows how the cost function E depends on the number of optimization steps for n-l of 3-3 and 4-4 (10 sampled runs each). Each distinct step in E represents sequential optimization from solution \( \psi (\theta ^k) \rangle\) at time-step k towards the solution \( \psi (\theta ^k+1) \rangle\) at \(k+1\). For small time-step \(\Delta t\), \(\theta ^k\) provides a good initial parameter set for solving optimization step \(k+1\). If the Ansatz parameters were re-initialized randomly \(\theta ^k\in [0,2\pi ]\) before each time-step, significantly more optimization steps would be required on average for convergence for each run (see Fig. 2b, inset).

Here we briefly examine the time complexity of the quantum algorithm excluding the classical computing components. Following the analysis of the variational Poisson solver18, the time complexity of the proposed variational evolution equation solver per time-step reads

where the terms within the inner parentheses indicate the time complexity of the state preparation scaling as \(\mathcal O(l+e+n^2 )\), which consists of the Ansatz depth l, the encoding depth \(e = \mathcal O(n^2)\)35, the depth of the circuit needed to implement the n-qubit cyclic shift operator \(O(n^2)\), and that of the number of shots \(\mathcal O(\varepsilon ^-2)\) required for estimation of expectation values up to a mean squared error of \(\varepsilon ^2\). The required number of quantum circuits depends on the boundary conditions applied (3 for periodic, 4 for Dirichlet and 5 for Neumann conditions), scaling only as \(\mathcal O(n^0 )\). \(\barT_\mathrm eval\) is the time-averaged number of function evaluations,

For deep and wide quantum circuits, the increase in optimization time is exacerbated by the presence of barren plateaus, or vanishingly small gradients in the energy landscape, where re-initialization can leave one trapped at a position far removed from the minimum37,38,39. Conversely, short time-steps lead to efficient solution trajectories that remain close to the local cost minima, leading to significant reduction in optimization times. To verify this, we conduct numerical simulations varying the diffusion parameter \(\delta\) with \(l=n\) up to time \(T=1\). Figure 3c shows that the number of iterations, or required optimization steps, per time-step increases linearly with \(\delta\). Close inspection of the time-averaged trace distance shows bimodal distributions at higher \(\delta\), which separates success and failure during convergence towards the global minimum (see Fig. 3d, dotted lines), resembling local minima traps due to poor optimization or expressivity of Anstze19,40.

3a8082e126