# A New Variational Quantum Algorithm Based on Lagrange Polynomial Encoding to
Solve Partial Differential Equations

Josephine Hunout  Department of Aeronautics, Imperial College of London,
London SW7 2AZ, United Kingdom Corresponding author<EMAIL_ADDRESS>Sylvain Laizet  Department of Aeronautics, Imperial College of London, London
SW7 2AZ, United Kingdom Lorenzo Iannucci  Department of Aeronautics, Imperial
College of London, London SW7 2AZ, United Kingdom

###### Abstract

Partial Differential Equations (PDEs) serve as the cornerstone for a wide
range of scientific endeavours, their solutions weaving through the core of
diverse fields such as structural engineering, fluid dynamics, and financial
modelling. PDEs are notoriously hard to solve, due to their the intricate
nature, and finding solutions to PDEs often exceeds the capabilities of
traditional computational approaches. Recent advances in quantum computing
have triggered a growing interest from researchers for the design of quantum
algorithms for solving PDEs. In this work, we introduce two different
architectures of a novel variational quantum algorithm (VQA) with Lagrange
polynomial encoding in combination with derivative quantum circuits using the
Hadamard test differentiation to approximate the solution of PDEs. To
demonstrate the potential of our new VQA, two well-known PDEs are used: the
damped mass-spring system from a given initial value and the Poisson equation
for periodic, Dirichlet and Neumann boundary conditions. It is shown that the
proposed new VQA has a reduced gate complexity compared to previous
variational quantum algorithms, for a similar or better quality of the
solution.

## 1\. Introduction

Solving equations and particularly partial differential equations (PDEs) is a
key element in various fields of science. In aerospace engineering, they are
essential to model and understand the dynamic behaviour of physical systems in
structural analysis or fluid dynamics. In classical computation, PDEs can be
solved using differentiation methods such as finite-difference, finite
element/volume methods, or using spectral methods like polynomial
approximations. With these methods, reaching a high accuracy often requires a
fine spatial discretisation or a large basis set, and in both cases, it can
lead to a significant computational cost. In engineering applications, despite
PDEs typically having low dimensions from 1D to 3D, the need for a very fine
mesh in all spatial directions is common, which can result in high memory and
computational requirements.

Quantum computing is a transformative new paradigm which takes advantage of
the quantum phenomena seen at the microscopic physical scale. While
significantly more challenging to design, quantum computers can run
specialised algorithms that can scale better than their classical
counterparts, sometimes exponentially faster. It is therefore natural to
investigate the potential of quantum computers for solving PDEs. Quantum
algorithms for tackling PDEs fall into two distinct groups: fully quantum
algorithms, which utilise quantum circuits to manipulate the quantum state in
accordance with the relevant PDE, and hybrid quantum-classical algorithms, in
which a quantum computer plays a role in a broader classical computing
process.

Fully quantum approaches do offer certain advantages in solving linear PDEs
because quantum operators act linearly on quantum states. However, the
evaluation of the performance and quantum advantage of these algorithms is
still unclear, as these metrics highly depend on the complexity of the problem
and it is often challenging to test the algorithms on real quantum machines.
Fully quantum algorithms face inherent limitations due to the necessity of
encoding and retrieving extensive classical data in a quantum superposition.
The advancement of this technology, known as quantum Random Access Memory
(qRAM) [1], is crucial for a significant application of such algorithms. The
current state of quantum hardware and the challenges associated with fault-
tolerant error-correction make the implementation of a reliable and large-
scale qRAM unfeasible at present [2, 3]. Although the quick evolution of
quantum computers and the symbolic limit of one thousand qubits reached with
the release of the IBM Condor quantum computer [4], a large number of physical
qubits do not guarantee a large number of logical qubits due to actual error-
correction techniques [5].

For these reasons, hybrid quantum-classical algorithms have been developed in
recent years. In particular, Variational Quantum Algorithms (VQAs) are notable
for their hybrid quantum-classical approach, making them suitable for
implementation on near-term quantum devices that are not yet fault-tolerant or
fully error-corrected. Such systems are often referred to as Noisy
Intermediate-Scale Quantum (NISQ) devices. The ”variational” aspect in VQAs
refers to the use of variational principles, which are a mathematical approach
used to find approximations to the lowest energy states of a system. This is
akin to finding the most stable configuration of a physical system. In the
context of quantum computing, these principles are used to find the state of a
quantum system that minimises a certain objective function, which is often
related to the problem being solved. The hybrid nature of VQAs, leveraging
both quantum and classical computing resources, makes them versatile and
adaptable to various types of problems and quantum hardwares. This
adaptability positions VQAs as a promising approach for exploring the
potential of quantum computing in the near term.

Among the earliest algorithms that demonstrated the potential of quantum
computers to solve specific problems more efficiently than classical
counterparts were the Shor’s algorithm [6] and the Grover’s algorithm [7].
Both of these algorithms were foundational in quantum computing and
highlighted the potential for quantum solvers to address specific
computational challenges more efficiently than classical methods. They spurred
significant interest in quantum computing research and development, leading to
the exploration of various quantum algorithms for solving a wide range of
problems, including those involving PDEs. Many of the first quantum solvers
developed in the 2000s were encoded by vectors of amplitudes and built on the
Quantum Amplitude Amplification and the Quantum Amplitude Estimation
Algorithms (QAAA and QAEA) [8]. Hence, the first quantum solvers for PDEs used
differentiation methods over discretised spaces of variables. In 2006, the
Kacewicz quantum algorithm [9] was developed to solve initial value problems
(IVP), i.e., ordinary differential equations (ODEs) with initial conditions.
In this quantum algorithm, the QAAA and QAEA [8] are used to estimate the mean
value of the ODE over numerous sets of points. Later, the Kacewicz quantum
algorithm has been improved to be applied to partial derivative equations
(PDEs) and simulate the solution of the steady state of the Navier-Stokes
equations [10], the Burgers’ equation [11], and the heat equation [12].
Another important quantum solver for PDEs is the famous HHL Algorithm, a
quantum solver developed in 2009 [13] which has demonstrated an exponential
speed-up over classical solvers in certain cases. Since then, several solvers
for linear PDEs have been proposed, with algorithms based on the finite
difference method [14, 15, 16, 17], the finite element method [18, 19], high-
order methods [20, 21], or spectral methods [22]. Quantum solvers have also
been used to solve non-linear PDEs, often by applying linearisation techniques
[23, 24, 25].

Regarding VQAs specifically, they have been studied for many applications,
including equation solvers, and have shown potential near-term quantum
advantages. The aim of VQAs, often referred as quantum neural networks, is to
classically train one or multiple parametrised quantum circuits to solve a
given problem. For example, a well-known VQA is the Variational Quantum
Eigensolver [26, 27, 28] designed to find the ground state energy of a given
system, which is particularly used in quantum chemistry and materials science.
VQAs have already been designed and tested to solve PDEs [29, 30, 31] or ODEs
[32].

One of the most challenging aspects of VQAs is their trainability. Training
variational circuits using gradient-based optimisers often prove daunting,
given the exponential disappearance of the gradient with the system’s size,
referred as the Barren plateaus phenomenon. Nevertheless, strategic
approaches, such as wise choices of observables, gradient determination
techniques and limiting the depth of the system, offer potential solutions to
overcome the barren-plateau obstacle [33, 34].

In this work, two architectures for a new VQA based on Lagrange polynomials
encoding are presented and compared to existing VQAs. To the best of our
knowledge, this is the first time that such an encoding is combined with a VQA
approach. The proposed strategy is based on differentiable quantum circuits
employing Lagrange polynomials and the Hadamard test differentiation method.
The comparison is conducted in terms of solution quality and algorithm gate
and circuit complexity (number of gates and circuits). Furthermore, it should
be noted that all VQAs discussed here are simulated, and the comparison is
carried out in a noise-free simulation environment. In order to assess the
potential of our new VQA approach, the focus of the present study is on two
PDEs: the damped mass-spring system from a given initial value and the Poisson
equation for periodic, Dirichlet and Neumann boundary conditions.

The paper is organised as follows: first, a background section defines VQAs
and highlights key elements such as quantum circuits differentiation methods,
and existing VQAs for PDEs, which serve as a basis for comparison with the VQA
architectures proposed in this article. Then, methods used for the development
of the novel VQA are presented, and finally results for both applications are
discussed. The paper is ended with a conclusion.

## 2\. Background

###  2.1 Variational Quantum Algorithms

VQAs have been developed to obtain a quantum advantage using NISQ devices,
i.e. with limited qubit resources, connectivity, and circuit depth, setting
them apart from fault-tolerant quantum algorithms. The objective of a VQA is
to ally the strengths of both classical and quantum computing paradigms, by
training a parameterised quantum circuit with a classical optimiser. While
VQAs have been explored across diverse applications, they may feature
different algorithmic approaches and quantum circuit architectures, but they
all share a general definition using common key elements [35]: a quantum
circuit built on a set of variational parameters ${\\{\theta_{i}\\}}_{i}$, a
cost operator $\hat{C}$, a loss function $\mathcal{L}$ and a classical
optimiser.

####  2.1.1 Variational Quantum Circuit

In the VQAs studied in this paper, the Variational Quantum Circuits (VQCs) are
composed of two parts: an encoding block, sometimes referred to as a quantum
feature map circuit and a variational ansatz, sometimes directly referred to
as the variational quantum circuit itself. While the encoding block is used to
initialise the circuit with fixed parameters, the variational ansatz is built
on the set of variational parameters $\\{{\theta_{i}\\}}_{i}$ later optimised
throughout the algorithm.

####  2.1.2 Cost Operator and Loss Function

One of the core features of a VQA is to select the architecture of the VQC and
a suitable cost operator $\hat{C}$ for a given problem. In this paper, the
cost operator is distinguished from the loss function $\mathcal{L}$. The cost
operator maps observables to the VQCs, while the loss function guides the
optimisation process by quantifying the distance between the algorithm’s
outputs and the desired solution.

####  2.1.3 Classical optimisers

Finally, the VQC is trained by a classical optimiser. As the loss function may
have multiple local extrema, training a VQC is recognised as NP-hard [36]
while facing new challenges due to the stochastic nature of quantum
computation through measurements, noise, and presence of the barren plateau
(detailed in 2.1.4) [35]. Consequently, the choice of a classical optimiser
remains an active area of research.

One of the initial classical optimisers used for VQAs is the Adam optimiser,
widely used in classical machine learning in the training of neural networks
due to its ease of implementation, computational efficiency, and little memory
requirements [37]. The Adam optimiser is an advanced version of the stochastic
gradient descent (SGD) which uses first-order gradients to estimate first and
second moments, thereby adapting the learning rate for each parameter
iteratively throughout the optimisation process. This optimiser finds
application in the VQA by Kyrrienko et al. [30], and is chosen for the VQA
proposed in this paper, detailed in Section 3. The second classical optimiser
mentioned in this paper through the VQA proposed by Sato et al. [32], is the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimiser. While the Adam optimiser is
designed to handle sparse gradients and noise, the BFGS aims to solve
unconstraint nonlinear optimisation problems. This optimiser is a quasi-Newton
approach which approximates the Hessian matrix of the loss function from the
first-order gradient or an approximation of the gradient.

The performance of optimisers for VQAs depends on the VQC structure and the
loss function. A recent comparison between classical gradient-based and
gradient-free optimisers [38] for the Variational Linear Quantum Solver [39]
over 3 to 5 qubits with realistic noise did not highlight an optimiser type
that outperformed the other. Among gradient-based optimisers, the best
performances were obtained with the BFGS optimiser, while the overall best
performances were obtained with the Simultaneous Perturbation Stochastic
Approximation (SPSA) [40], a gradient-free optimiser. This optimiser is
considered as gradient-free as it approximates gradients from a single partial
derivative obtained with finite differencing along a random direction. The
SPSA is expected to be an efficient method for VQAs due to its reduced
complexity and better performance compared to a SGD approach. However, other
gradient-free optimisers such as COBYLA [41] performed poorly in the presence
of noise compared to other optimisers [38].

Since these latter classical optimisers were initially proposed for classical
computation, other methods have been developed recently for quantum hybrid
computation, with, for example, the meta-learner [42], a SGD-based optimiser,
which adapts the learning rate of the optimisation process by training a
neural network based on the optimisation history of similar problems. Another
adaptive and SGD-based approach for VQA is to variate the number of shots
required to determine the loss function gradient, that is, its precision,
throughout the optimisation process to significantly reduce the complexity of
the optimisation [43, 44]. Another gradient-based approach for VQAs is the
quantum natural gradient descent method, which instead evaluates the gradient
of the loss function in terms of the Euclidean norms. It uses a metric tensor
that quantifies the sensitivity of the quantum circuit to variation of the
parameters [45]. Finally, in the case of a loss function defined as a sum of
trigonometric functions, a gradient-free approach can be used where all the
local parameters are sequentially updated [46, 47].

####  2.1.4 Barren Plateau

The Barren plateau phenomenon refers to the problem of vanishing gradients in
VQAs. It severely limits the trainability of these algorithms and poses a
challenge for optimisation problems in quantum systems. This phenomenon
continues to receive a significant attention from researchers, to
comprehensively understand the phenomenon, ascertain the reasons behind its
occurrence, and explore effective strategies to surmount its challenges.

The first cause identified for the Barren plateau phenomenon is the use of
deep randomly initialised hardware efficient ansatzes [48]. Consequently,
attention to overcome the barren plateau was initially directed toward a
clever initialisation of the VQC parameters [49, 50, 51, 52] or a clever VQC
structure [53]. Later, Cerezo et al. [33] studied the correlation between the
depth of the circuit as well as the use of a local or global measurement and
the occurrence of barren plateau. Their findings revealed that even shallow
VQCs could lead to a barren plateau, highlighting the selection of the cost
function as a main cause of the phenomenon. Specifically, opting for local
cost functions over global cost functions was shown to mitigate the barren
plateau, resulting in a polynomial gradient decay with circuit size instead of
an exponential decay. Recently, another approach proposed counteracting the
emergence of the barren plateau by computing the gradient and variance of the
cost function using Clifford approximant quantum circuits [34].

###  2.2 Derivative Quantum Circuits

Analytical derivatives of quantum circuits with respect to Pauli gates
parameters can be obtained easily by adding a corresponding rotational gate,
depending on the circuit structure.

####  2.2.1 Parameter Shift Rule Method

The Parameter Shift Rule method is widely used to differentiate quantum
circuits. This method requires to build two modified circuits to obtain the
derivative with respect to one Pauli gate parameter [29]. In the example of a
circuit of two qubits as described in Figure 1.a, the derivative with respect
to the parameter gate $\alpha$ is obtained by the difference of two modified
versions where the parameter gate $\alpha$ is shifted to $\pm\pi/2$, as
described in Figure 1.b.

Figure 1: a) Original quantum circuit, b) Derivative quantum circuits

Thus, the derivative of the original quantum circuit in figure 1.a with
respect to $\alpha$ is:

$\frac{\partial}{\partial\alpha}\expectationvalue{\hat{C}}{f_{\alpha}}=\frac{1}{2}\left(\expectationvalue{\hat{C}}{f_{\alpha^{+}}}-\expectationvalue{\hat{C}}{f_{\alpha^{-}}}\right),$
(1)

where $\ket{f_{\alpha}}$ is the state vector of the original quantum circuit,
$\ket{f_{\alpha^{+}}}$ and $\ket{f_{\alpha^{-}}}$ are the state vectors of the
quantum circuit where the parameter $\alpha$ has been shifted respectively of
$\pi/2$ and $-\pi/2$. $\hat{C}$ is a chosen cost function.

####  2.2.2 Hadamard Test Differentiation

A second method to derive quantum circuits is based on the Hadamard test,
originally used to extract the real part of an expectation value by measuring
the magnetisation of the ancilla qubit $a$, as shown in Figure 2.a:
$\bra{a}\hat{Z}\ket{a}=Re\bra{\psi}\hat{\mathcal{U}}\ket{\psi}$.

Figure 2: a) Hadamard Test, b) Hadamard Test Differentiation

This structure can be reused to obtain the derivative of the expectation value
real part with respect to a Pauli gate parameter $\theta_{j}$ (Figure 2.b)
[54]:

$\bra{a}\hat{Z}\ket{a}=\frac{\partial}{\partial\theta_{j}}Re\bra{\psi}\hat{\mathcal{U}}\ket{\psi},$
(2)

where $\hat{R}_{j}=\exp(\frac{-i}{2}\theta_{j}\hat{P}_{j})$, and $\hat{P}$ is
a tensor of Pauli matrices. With this differentiation method, only one circuit
is required to obtain the derivative of the expectation value, but only the
real value of the expectation is extracted.

###  2.3 Kyriienko et al. Algorithm

The first quantum algorithm studied in this work for comparison purposes is a
VQA based on the algorithm developed by Kyriienko et al. [30], which trains
VQCs to approximate the solution of PDEs by polynomial fitting. The Kyriienko
et al. algorithm has been originally designed to solve first-order
differential equations, but in this work, its application is expanded to
second-order differential equations.

####  2.3.1 Circuit Structure

Figure 3: Circuit Structure of the Kyriienko et al. algorithm.

The Kyriienko et al. algorithm [30] aims to represent an approximation of the
solution of a differential equation by a parameterised quantum circuit divided
into two parts : a quantum feature map $\mathcal{\hat{U}}_{\varphi}(x)$ and a
variational ansatz $\mathcal{\hat{U}}_{\theta}$. The first part encodes the
equation variable $x$ through a fitting set of nonlinear Chebyshev polynomials
of the first and second kind over $]0,1[$. Encoding the variable $x$ through
functions allows a high expressivity of the quantum circuit while ensuring its
derivability. The quantum feature map $\mathcal{\hat{U}}_{\varphi}(x)$ is
composed by Y-rotational gates $\hat{R}_{Y}$, parameterised by the functions
$\varphi_{j}(x)=2j\arccos{(x)}$, where the degree $j$ of the Chebyshev
polynomials grows with the number of qubits $n$ [30].
$\hat{\mathcal{U}}_{\varphi}(x)$ can be expressed as

$\hat{\mathcal{U}}_{\varphi}(x)=\bigotimes_{j=1}^{n}\hat{R}_{Y,j}\left[2j\arccos{(x)}\right]$
(3)

This quantum feature map encodes the Chebyshev polynomials of the first and
second kind as follow:

$\hat{R}_{Y,n}\big{(}\varphi(x)\big{)}\ket{0}=T_{n}(x)\ket{0}+\sqrt{1-x^{2}}\
U_{n-1}(x)\ket{1},$ (4)

where the Chebyshev polynomials of the first kind is defined by
$T_{n}=\cos(n\arccos(x))$ and of the second kind by
$U_{n-1}(x)=\sin(n\arccos(x))/\sin(\arccos(x))$ [30]. The chosen variational
ansatz $\hat{\mathcal{U}}_{\theta}$ is the hardware efficient ansatz [55],
composed of multiple layers of a sequence of rotational gates and controlled
NOT (CNOT) gates.

####  2.3.2 Algorithm

In Kyriienko et al. [30], the resulting quantum state
$\ket{f_{\varphi,\theta}(x)}=\hat{\mathcal{U}}_{\varphi}(x)\hat{\mathcal{U}}_{\theta}\ket{\text{\O}}$
is evaluated by the total magnetisation of the system
$\hat{C}=\sum_{j=1}^{n}\hat{Z}_{j}$. The derivatives with respect to $x$ and
$\theta$ are determined by the parameter shift rule method and the chosen
classical optimiser is Adam [37]. Finally, the loss function is defined as the
mean square of the PDE with all terms on one side, as expressed in Equation 6.

###  2.4 Sato et al. Algorithm

The second VQA used in this work for comparison purposes is the Sato et al.
algorithm [32], which was designed to solve the Poisson equation. This
algorithm is based on a discretised approach in which the VQC is used to
approximate the unitary vector of the solution of the 1D Poisson equation.

####  2.4.1 Circuit Structure

The circuit structure has been built to solve a 1D second order differential
equation of the form : $d^{2}f(x)/dx^{2}+s(x)=0$, where $s$ is the source term
encoded through its amplitudes by $\hat{\mathcal{U}}_{s}$ and the direction
vector of the solution $f$ is approximated by the variational ansatz
$\hat{\mathcal{U}}_{\theta}$. This structure using an ancilla gives the
following state :

$\ket{s,\psi_{\theta}}=\sqrt{2}(\ket{0}\ket{s}+\ket{1}\ket{\psi_{\theta}})$
(5)

where the state vector $\ket{\psi_{\theta}}$ encodes the unitary vector of the
discretised vector solution $f$.

Figure 4: Circuit Structure of the Sato et al. algorithm.

####  2.4.2 Algorithm

This approach evaluates the VQC described in Figure 4 with multiple cost
operators depending on the boundary conditions to quantify the loss function,
which here is the potential energy of the PDE. As this algorithm is a
discretised approach, derivatives with respect to the variable $x$ are
obtained with a finite-difference scheme. In contrast, the gradient of the
loss function $\nabla_{\theta}\mathcal{L}_{\theta}$ requires differentiating
the VQC, using a method based on the Hadamard test Differentiation [32]. The
classical optimiser chosen by Sato et al. is the BFGS optimiser, described in
Section 2.1.3.

## 3\. Methods

###  3.1 Algorithms Workflow

For both algorithms studied here, the quantum circuits are aimed to
approximate the solution of a PDE. To do so, the algorithms follow a workflow
inspired by Kyriienko et al. [30], summarised in Figure 5. First, an initial
set of quantum circuits and their derivatives with respect to $x$ is built for
each point $x_{i}$ in the chosen interval of interest $I$, with a common
random vector of variational parameters $\theta$. The value of the resulting
function is read as an expectation value through a cost operator $\hat{C}$
applied to the quantum state $\ket{f_{\varphi,\theta}(x)}$ of the circuit:
$f_{\theta}(x)=\expectationvalue{\hat{C}}{f_{\varphi,\theta}(x)}$. The
classical optimiser uses a loss function, which evaluates if the resulting
function satisfies the PDE as well as other conditions like the boundary
conditions or regularisation points. In this work, the selected classical
optimiser is Adam [37], which demonstrated satisfactory performances in
previous studies [30]. As this optimiser relies on gradient descent
optimisation, derivative quantum circuits with respect to the variational
parameters $\theta$ must be calculated. Finally, the vector of variational
parameters $\theta$ is updated and all circuits are recalculated in another
loop until the resulting function converges to the desired function.

This work specifically focuses on enhancing the VQC structure to not only
encode other polynomials but also to effectively determine the gradient and
derivatives at a lower computation cost by comparison to existing approach
while not compromising on the quality of the solution, thereby making this new
algorithm attractive.

Figure 5: Proposed algorithm workflow.

###  3.2 Loss Function

With a similar approach to the Kyrrienko et al. algorithm [30], the loss
function $\mathcal{L}_{\theta}^{DE}$ quantifies the distance between the PDE -
written with all terms on one side - and zero, over a set of $M$ points within
the chosen interval $I$. This loss function is defined as

$\mathcal{L}_{\theta}^{DE}\left[\frac{\partial^{2}f_{\theta}}{\partial
x^{2}},\frac{\partial f_{\theta}}{\partial
x},f_{\theta},x\right]=\frac{1}{M}\sum_{i=1}^{M}L\left(DE\left[\frac{\partial^{2}f_{\theta}}{\partial
x^{2}}\left(x_{i}\right),\frac{\partial f_{\theta}}{\partial
x}\left(x_{i}\right),f_{\theta}\left(x_{i}\right),x_{i}\right],0\right),$ (6)

where the distance $L$ is evaluated with the mean square error
$L(a,b)=(a-b)^{2}$. Other distance evaluations can be used, such as the mean
absolute error $L(a,b)=|a-b|$. This loss function can be completed with
additional functions that quantify the satisfaction of boundary conditions
$\mathcal{L}_{\theta}^{BC}$ or regularisation points
$\mathcal{L}_{\theta}^{R}$ [30]. Arbitrary weights $\eta$ can also be used to
prioritise a condition over the others, as expressed in the equation

$\mathcal{L}_{\theta}=\eta_{DE}\ \mathcal{L}_{\theta}^{DE}+\eta_{BC}\
\mathcal{L}_{\theta}^{BC}+\eta_{R}\ \mathcal{L}_{\theta}^{R}.$ (7)

###  3.3 Boundary Conditions

The boundary conditions can be taken into account in the optimisation loop by
contributing to the loss function, or can also be handled separately as an
automatic shift of the function read out from the quantum circuit as proposed
by Kyriienko et al. [30].

####  3.3.1 Contribution to the loss function

For a second-order PDE, the contribution of the boundary conditions to the
loss function can be described as follows:

$\mathcal{L}_{\theta}^{BC}=L\left(f_{\theta}(x_{0}),u_{0}\right)+L\left(\frac{\partial
f_{\theta}}{\partial x}(x_{0}),\dot{u}_{0}\right),$ (8)

where $u_{0}$ and $\dot{u}_{0}$ are the value of the sought function and its
derivative at $x=x_{0}$. For a lower- or higher-order PDE, the corresponding
terms must be removed or added.

####  3.3.2 Floating boundary handling

As an alternative approach, a shift term $f_{BC}$ can be calculated at each
iteration to automatically match the boundary condition regardless of the
value of the variable parameters $\theta$ [30]. Hence, the contribution to the
loss function is null ($\mathcal{L}_{\theta}^{BC}=0$), and we have

$f_{\theta}(x)=\expectationvalue{\hat{C}}{f_{\varphi,\theta}(x)}+f_{BC}\ ,$
(9) $f_{BC}=u_{0}-\expectationvalue{\hat{C}}{f_{\varphi,\theta}(x_{0})}.$ (10)

This definition can be extended for the derivatives of the sought function for
all PDEs. However, to ensure the flexibility of the algorithm and its
convergence toward the sought function, the boundary condition can be handled
with a shift term for the function and a loss contribution for the
derivatives.

###  3.4 Regularisation

For equations where some values of the desired function $\\{u_{r}\\}$ are
known, an additional condition can be created to ensure the convergence of the
desired function toward those values. This contribution to the loss function
can be defined as

$\mathcal{L}_{\theta}^{R}=\frac{1}{M_{R}}\sum_{r=1}^{M_{R}}\left(L(f_{\theta}(x_{r}),u_{r})\right),$
(11)

where $M_{R}$ is the number of regularisation points and $u_{r}$ the value of
the sought function at those points. A regularisation procedure can also be
implemented for derivatives values [30]. To minimise the complexity of the
algorithm, regularisation can be added to the loop by dividing the set of
points $\\{x_{i}\\}_{i=1}^{M}$ into subsets of points that will, one by one,
be used as training points for the differential equation and then as
regularisation points once the algorithm converged on this subset. This
approach can save an important amount of computing resources but can in
certain cases decrease the overall accuracy of the algorithm output. In the
first application presented in Results section (4.1), regularisation have been
used for the Hadamard-Lagrange algorithm, to reduce its overall gate
complexity.

###  3.5 Cost Operator

The cost operator $\hat{C}$ is used to read out the value of the approximated
function and its derivatives from the different quantum circuits built
throughout the algorithm, with

$f(x)=\bra{f_{\varphi,\theta}(x)}\hat{C}\ket{f_{\varphi,\theta}(x)},$ (12)

where $\ket{f_{\varphi,\theta}(x)}$ is the quantum state of the quantum
circuit evaluated at $x$ within the chosen interval. The choice of the cost
operator impacts the flexibility and trainability of the algorithm. Most
common cost functions are the magnetisation of a single qubit $j$ :
$\hat{C}=\hat{Z}_{j}$, which can be used for the measurement of an ancilla
qubit for example, or the total magnetisation of the system, where all qubits
are measured: $\hat{C}=\sum_{j=1}^{n}\hat{Z}_{j}$.

In this study, the cost functions employed are centred around the total
magnetisation of the system, or variations of it. The decision to opt for
projective measurement along the Z-axis is motivated by its straightforward
implementation and the potential to simplify the variational ansatz, reducing
the overall gate complexity of VQAs as detailed in appendix A.

###  3.6 Derivatives and Gradients

The two new VQA architectures proposed in this paper aim to represent the
function solution of the PDE through a quantum circuit. To do so, the
derivatives with respect to the variable $x$ are necessary to evaluate the
loss function. Both quantum feature maps encode the variable $x$ using a set
of functions $\\{\varphi_{i}\\}_{i=0}^{n}$ where $n$ is the number of gates
used for the quantum feature map encoding. Hence, the first- and second-order
derivatives are defined as follow:

$\frac{df}{dx}=\sum_{i=0}^{n}{\frac{\partial}{\partial\varphi_{i}}\bra{f_{\varphi,\theta,i}}\hat{C}\ket{f_{\varphi,\theta,i}}\
\frac{d\varphi_{i}}{dx_{i}}},$ (13)
$\frac{d^{2}f}{dx^{2}}=\sum_{i=0}^{n}{\frac{\partial}{\partial\varphi_{i}}\bra{f_{\varphi,\theta,i}}\hat{C}\ket{f_{\varphi,\theta,i}}\
\frac{d^{2}\varphi_{i}}{dx_{i}^{2}}}+\\\
\sum_{i=0}^{n}{\sum_{j=0}^{n}{\frac{\partial^{2}}{\partial\varphi_{i}\varphi_{j}}\bra{f_{\varphi,\theta,i,j}}\hat{C}\ket{f_{\varphi,\theta,i,j}}\
\frac{d\varphi_{i}}{dx_{i}}\frac{d\varphi_{j}}{dx_{j}}}}.$ (14)

However, in the variational ansatz, the variational parameter $\theta$ is
encoded directly into the parameterised rotational gates. In this work, all
gradients $\nabla_{\theta}$ are determined using the Parameter Shift Rule
method as defined in equation 1.

###  3.7 Quantum Circuits Structures

For the two new VQA architectures, the variable function $x$ is encoded
through a quantum feature map circuit $\hat{\mathcal{U}}_{\varphi}(x)$ and the
vector of variational parameters $\theta$ is encoded through a variable ansatz
$\hat{\mathcal{U}}_{\theta}$. The VQA inspired by the work of Kyriienko et al.
[30] is based on two circuits acting on the same register, one after the
other, and all qubits in the register are measured, as shown in Figure 3. For
the newly designed VQA, two different structures have been studied (Figures 6
and 7) and in both cases the circuit is inspired by the Hadamard test
structure (Figure 2), over the ancilla register of one qubit or multiple
qubits.

The development of the Hadamard-Lagrange algorithm has been driven by two
primary factors. Firstly, due to the computational demands of the parameter
shift rule method, the gate complexity involved in determining the first and
second derivatives with respect to $x$, along with their gradients with
respect of $\theta$, was substantial. Secondly, observations from the initial
application of the Kyriienko-inspired algorithm to a second-order DE, as
outlined in Results section (4.1), revealed oscillations which may be caused
by the use of excessively high-degree polynomials for interpolation tasks.
Consequently, the present study focused on a VQC structure that not only
reduces the computational overhead associated with evaluating derivatives and
gradients but also implementing the Lagrange interpolation, while not
impacting the quality and accuracy of the solution. One of the main objectives
was also to perform the interpolations with the lowest possible degree of
polynomial.

The first approach was initially designed across two quantum registers of the
same sizes to ensure the VQC’s trainability. As discussed in Background
Section (2.1.4), local measurements can mitigate the barren plateau
phenomenon. However, for circuits with very shallow depths, this expanded
structure is suboptimal as it necessitates a high number of qubits. Therefore,
the Hadamard-Lagrange structure was re-designed to minimise the required
number of qubits for encoding, while retaining the same cost function and
differentiation method (second approach). Intermediate structures where the
number of qubits required for the encoding is between 1 and $n$ could also
have been considered, but the study of such structure is left for a future
work.

###  3.8 Kyriienko-Inspired Algorithm for comparison

In order to assess the performance of our new VQA, the algorithm developed by
Kyriienko et al.[30] is expanded following Equation 14. The hardware efficient
ansatz originally used for the variational ansatz is replaced by two layers of
a single X-rotational gate $\hat{R}_{X}$ per qubit and a linear network of
CNOT gates to minimise the complexity of the classical optimiser and ensure
the convergence of the algorithm as a complete hardware efficient ansatz [55]
structure is not required with a projective measurement onto to Z-axis (see
Appendix A for more details).

###  3.9 Proposed new Hadamard-Lagrange Algorithm

The new algorithm designed and tested in this work has been developed based on
the Hadamard test structure (Figure 2) but extended to a multiple-qubit
register, to approximate the solution of the differential equation by Lagrange
interpolation. In this work, two different structures of the quantum feature
map have been studied.

####  3.9.1 Quantum Feature map

First, an extended version has been developed where an encoding block
$\hat{\mathcal{U}}_{\varphi}$, with
$\varphi_{i}(x)=\arccos(\frac{x-x_{j}}{2})$, is applied in between two
networks of CNOT gates, over a quantum register of $n$ qubits, leading to

$\hat{\mathcal{U}}_{\varphi}(x)=\bigotimes_{j=1}^{n}\hat{R}_{Y,j}\Big{(}\arccos(\frac{x-x_{j}}{2})\Big{)}.$
(15)

Figure 6: Circuit Structure of the Hadamard Lagrange Algorithm with an
extended structure.

A second version with a simplified structure has been developed, where only
one ancilla qubit is used for the encoding. Here, the variable $x$ is encoded
through Y-rotational gates parameterised by the same encoding function
$\varphi_{i}$ defined in the extended structure but within the networks of
CNOT gates.

Figure 7: Circuit Structure of the Hadamard Lagrange Algorithm with a
simplified structure.

In both structures, the set of points $\\{x_{j}\\}_{j=1}^{n}$ included in the
interval $I$ used for the polynomial interpolation is directly linked to the
size of the first quantum register. While the second register can be built
over $n$ qubits (Figure 6) or only one qubit (Figure 7), for both structures,
only the first register is measured, leading to local measurements instead of
global measurements. As discussed in Background section 2.1.4, this allows to
delay the appearance of the barren plateau phenomenon.

Additionally, determining whether to use a simplified or expanded circuit
structure, depends on the quantum device being employed. The expanded circuit
structure requires more qubits but fewer gates in comparison to its simplified
counterpart. Consequently, while it reduces noise from gate errors or
decoherence, this structure may be prone to errors from qubits or cross-talk,
and can possibly be limited by the existing quantum resources available in
NISQ devices. Conversely, the simplified structure shows promise, particularly
with the development of efficient mitigation techniques and error gate
correction. In the realm of quantum circuit simulations, the preference leans
towards the simplified circuit structure due to its lower computational
resource demands (as both circuit provides the same solution on noise-free
qubits).

####  3.9.2 Variational Quantum Circuit

In our new VQA, the variational ansatz $\hat{\mathcal{U}}_{\theta}$ is
composed of one layer, following the same structure as presented in the
Kyriienko-inspired Algorithm which is used for comparison (see section 3.8).

####  3.9.3 Cost Operator

For both circuit structures to be read as Lagrange interpolating polynomial,
the first register is evaluated by the same cost operator:

$\hat{C}=\sum_{j=1}^{n}\frac{1}{\rho_{j}}\hat{Z}_{j},\text{ with }\
\rho_{j}=\bra{f_{\varphi,0}(x_{j})}\hat{Z}_{j}\ket{f_{\varphi,0}(x_{j})}.$
(16)

Hence, the resulting function is:

$f(x)=\sum_{j=1}^{n}\left(\alpha_{\theta}\prod_{i\neq
j}^{n}\frac{x-x_{i}}{x_{j}-x_{i}}\right),$ (17)

where $\alpha_{\theta}$ is a coefficient determined by the variational ansatz.

## 4\. Results

In order to assess the performance of our new Hadamard-Lagrange algorithm, two
PDEs will be used, as well as existing VQAs and analytical solutions for
comparisons, which will be performed with a focus on two metrics: the quality
of the solution and the gate complexity of the VQAs.

###  4.1 Damped Mass Spring System

A Damped Mass-Spring System (DMSS) is a physical model used to describe the
behaviour of a mass attached to a spring when it is subject to a damping
force. This system is widely studied in mechanics and physics to understand
oscillatory motion and damping effects. In this work, we focus on an initial
value problem (IVP) of a DMSS solved over a time interval of $I=[0,10]$
seconds, using the Kyriienko-inspired algorithm and the proposed new Hadamard-
Lagrange algorithm with an expanded circuit structure. The PDE can be
expressed as

$m\frac{d^{2}f}{dt^{2}}+b\frac{df}{dt}+kf=0,$ (18)

with $m=1kg$; $b=1kg.s^{-1}$; $k=1N.m^{-1}$; $u_{0}=1m$ and
$\dot{u}_{0}=0m.s^{-1}$. The results obtained by the two quantum methods are
compared with the analytical solution, detailed in the appendix C.

In this study, two quantum approaches are applied to the Damped Mass Spring
System differential equation, and analysed by looking at their performance
using two distinct sets of training points. These sets are based on Chebyshev
nodes, recognised as optimal for avoiding the Runge phenomenon [56] when
interpolating over equispaced points induces oscillations at the edges of the
interval [56]. The first set $\\{x_{i}^{(1)}\\}_{i}$ corresponds to the
original Chebyshev nodes (originally between $[-1,1]$) rescaled over an
interval $[a,b]$ (with $0\leq a<b\leq 1$). The second set
$\\{x_{i}^{(2)}\\}_{i}$ corresponds only to the original positive Chebyshev
nodes scaled over an interval $[0,b]$. Here, $a=0$ and $b=0.9$ as imposed by
qubits encoding and for stability reason, for both sets of points. These two
sets can be expressed as:

$x_{k}^{(1)}=\frac{a+b}{2}+\frac{b-a}{2}\cos\left({\frac{2k-1}{2N_{points}}\pi}\right),0\leq
k<N_{points}$ (19)

$x_{k}^{(2)}=b\cos\left({\frac{2k-1}{4(N_{points}-1)}\pi}\right),0\leq
k<N_{points}$ (20)

To assess the quality of the result and facilitate a comparative analysis of
the two quantum algorithms, we will differentiate the DE loss contribution at
each individual point, represented by the square of the DE evaluated at that
point, and, the total loss as the average of DE loss contributions across all
points considered, as defined in Equation 7. Here, the loss is assessed over
50 equispaced points within the specified interval.

To estimate the gate complexity between the two quantum approaches, the number
of circuits and basic quantum gates required for one iteration and the entire
calculation are compared. The detailed methodology for estimating this gate
complexity is provided in Appendix B. However, these results can be highly
dependent on the initial random values of the vector $\theta$, providing only
an estimate of the order of magnitude of the quantum circuits and quantum
gates required to solve the equation.

####  4.1.1 Kyriienko-inspired Algorithm

The application of the Kyriienko-inspired (KI) algorithm to the DMSS problem
has been investigated using a range of VQCs with 3 to 5 qubits and 2 to 7
layers, across 7, 9, or 12 points for both sets defined in Equations 19 and
20. Each VQC is initialised with a random vector of variational parameters
$\theta$. To ensure a fair comparison of the influence of the size and type of
training set, the same initial vector $\theta$ has been used for a given size
and depth of the VQC. Furthermore, the quantum feature map used in the KI
algorithm is based on high-degree Chebyshev polynomials, which can cause
oscillations at the extremities of the algorithm. To mitigate this issue, the
KI algorithm has been trained over an interval of $[0,12]$ seconds, ensuring a
smooth solution over the required interval of $[0,10]$ seconds.

Additionally, the quantum feature map of this algorithm’s structure is
independent of the training set of points, and increasing the number of
training points throughout the algorithm does not guarantee a more efficient
optimisation. In fact, more iterations may be necessary to achieve
convergence. Hence, the VQCs have been trained over all points in the chosen
set, using the Adam classical optimiser [37] with a learning rate of
$\alpha=0.01$. The algorithm is considered to have converged when the absolute
value of each component of the loss function gradient falls below a specified
threshold of $\epsilon_{\partial\mathcal{L}}=10^{-4}$, with:

$|\frac{\partial\mathcal{L}}{\partial\theta_{i}}|\leq\epsilon_{\partial\mathcal{L}}.$
(21)

For this application, the boundary and initial conditions are implemented with
a floating boundary handling (3.3.2) for $u_{0}$, in combination with a
contribution to the loss function (3.3.1) for $\dot{u}_{0}$. This ensures the
satisfaction of the first condition while avoiding to over-constraint the
system. In this case, the loss function can be expressed by Equation 7, with
$\eta_{DE}=1$ ; $\eta_{BC}=1$ and $\eta_{R}=0$.

Among the different VQCs studied for the KI algorithm, the lowest average DE
losses were observed with a VQC of 5 qubits and 2 layers over a set of 12
points of the second kind and with a VQC of 4 qubits and 3 layers over a set
of 12 points of the first kind. According to the results obtained from this
algorithm, an increased number of points in the training set does not always
correlate with a higher gate complexity or highly improved accuracy. However,
the selection of the training points themselves and the structure of the VQC
seems to have the most significant impact on both the quality of the solution
and the gate complexity, as shown in the Figure 8.

Figure 8: Average total loss and gate complexity for the Kyriienko-inspired
algorithm applied to the most efficient VQC structures with varying training
set lengths over 5 tries. $N5L2I2$ stands for the structure of 5 qubits and 2
layers across a training set of the second kind, while $N4L3I1$ stands for the
structure of 4 qubits and 3 layers across a training set of the first kind.

Additionally, applying this algorithm with different random initial vectors of
variational parameters $\theta$ has demonstrated a high sensitivity to these
initial parameters. Across five different attempts to solve this DMSS problem
with the structures presented above, a coefficient of variation of $20\%$ in
the total loss was observed for $N5L2I2$, along with a coefficient of
variation of $57\%$ in the number of iterations required for the algorithm to
converge. Similarly, for $N4L3I1$, a $50\%$ coefficient of variation in total
loss and a $40\%$ coefficient of variation in iterations were observed.
Consequently, the accuracy and gate complexity shown in Figure 8 only reflect
the average outcome of the algorithm, not its worst or best results. This
variability in the accuracy of the results is illustrated in Figure 9, which
shows the best result achieved for each VQC structure.

For the two VQC structures that produced the best results, some results
presented important oscillations toward the end of the interval, despite
placing some training points outside the interval to mitigate this issue.
Additionally, for both structures, the boundary conditions for the first- and
second-order derivatives were not satisfied. Finally, the corresponding total
DE loss obtained for the best results highlighted in Figure 9 is,
respectively, $2.50\cdot 10^{-3}$ for $N5L2I2$ and $2.77\cdot 10^{-3}$ for
$N4L3I1$. However, the boundary conditions for the first- and second-order
derivatives are not fully met, resulting in a BC loss of $0.447$ and $0.485$
for each VQC structure, respectively, as illustrated in Figure 12.

Figure 9: Results are presented for each of the most efficient VQC structures,
illustrating variability and the best outcome achieved across five attempts
for: a) the function $f$, b) the first derivative $\dot{f}$, c) the second
derivative $\ddot{f}$ and d) the DE loss contribution across 50 equispaced
points in the interval, plotted on a logarithmic scale. $N5L2I2$ stands for
the structure of 5 qubits and 2 layers across a training set of the second
kind, while $N4L3I1$ stands for the structure of 4 qubits and 3 layers across
a training set of the first kind.

####  4.1.2 Hadamard-Lagrange algorithm

The new algorithm proposed in this work introduces the use of interpolation
points in the encoding part of the VQC. These interpolation points are
selected to be similar to the training set of points. This approach avoids
over-constraining the optimisation process and helps reduce the overall gate
complexity. This new method has been studied using sets of 5 to 8 points.
Among these, the optimal VQC size was found to be based on 7 points,
corresponding to the use of 8 qubits, with the use of the simplified structure
as illustrated in Figure 7. As already mentioned, the focus of this work is on
the theoretical aspects of variational quantum algorithms, hence why the
quantum circuits examined in this work are simulated in a noise-free
environment. Consequently, the use of the simplified structure is preferred
because of its lower computational resource demands.

It should be noted that, for this algorithm, training points in pairs with the
number of qubits can be modified during the calculation. To accelerate the
optimisation, the process is divided into two parts, using an evolving set of
training points and regularisation points. Furthermore, the boundary
conditions are implemented following the same methodology as for the
Kyriienko-inspired algorithm, with this time $\eta_{DE}=1$ ; $\eta_{BC}=0.6$
and $\eta_{R}=1$.

Figure 10: Results of the new Hadamard-Lagrange VQA for both parts of the
optimisation process, achieved across five attempts. a) the function $f$, b)
the first derivative $\dot{f}$, c) the second derivative $\ddot{f}$ and d) the
DE loss contribution across 50 equispaced points in the interval, plotted on a
logarithmic scale. $I1$ stands for the first kind of training set of points
while $I2$ stands for the second.

For the first part of the optimisation, the set of training points initially
contains only the three Chebyshev nodes on the left of the domain (using 4
qubits). The first node contributes to the BC loss function
$\mathcal{L}_{\theta}^{BC}$ and the other two points contribute to the DE loss
function $\mathcal{L}_{\theta}^{DE}$. Once the algorithm converged for these
two nodes, which means that all components of the loss function gradient
reached a threshold value of $10^{-4}$ (defined in Equation 21) the next node
is added and will contribute to the DE loss function
$\mathcal{L}_{\theta}^{DE}$, while the second point will become a
regularisation point contributing to the regularisation loss function
$\mathcal{L}_{\theta}^{R}$. This process is repeated until convergence is
achieved for all nodes (using 8 qubits in this example), using the Adam
classical optimiser [37] with an evolving learning rate from $\alpha=0.04$ to
$\alpha=0.01$ depending on the total loss function value
$\mathcal{L}_{\theta}$. In the second part of the optimisation, the size of
the VQC do not evolve and the nodes contribute three by three to the DE loss,
whilst other points act as regularisation points, starting again from the left
side of the interval. During this second part, the Adam classical optimiser
[37] is applied with a fixed learning rate of $\alpha=0.01$.

The application of this new VQA, using this two-part optimisation, has been
observed to be insensitive to the initial random vector of variational
parameters $\theta$. For each type of training set, the solutions obtained
after the first part of the optimisation over 5 trials with different initial
random values of the vector of variational parameters $\theta$, were extremely
consistent, with a mean total loss of $1.24\cdot 10^{-2}$ and $2.76\cdot
10^{-2}$ with a coefficient of variation of $1.0\%$ and $1.5\%$ for the first
and second kind of training set, respectively. Similarly, the average number
of iterations to achieve the convergence of the first part is, respectively,
$670$ and $620$ iterations with a coefficient of variation of $11\%$ in both
cases. The second part of the optimisation process consistently starts from a
similar state, resulting in comparable deviation in terms the solution
accuracy and the number of iterations.

As illustrated in Figure 10, the first part of the optimisation yielded a
close approximation of the solution, achieving slightly better accuracy for
the initial set of training points. Similarly, the second part of the
optimisation process produced improved results for the first training set,
with a total DE loss of respectively $1.51\cdot 10^{-3}$ and $2.87\cdot
10^{-3}$. Additionally, the boundary conditions are more effectively met this
time, with a BC loss of respectively $1.18\cdot 10^{-3}$ and $1.07\cdot
10^{-4}$, as depicted in the Figure 12.

####  4.1.3 Hadamard-Lagrange algorithm vs Kyriienko-inspired algorithm

For both quantum algorithms, the number of qubits and Chebyshev nodes is
different. The idea here was to select, by trial and error, the optimal set of
training points and structures to obtain the most accurate solution for each
approach.

As shown in Figure 11, both quantum algorithms provide a satisfactory
approximation of the desired function $f_{\theta_{opt}}$ over the chosen
interval, with better accuracy with the Hadamard-Lagrange algorithm.
Noticeable oscillations in the general approximation of the derivatives with
respect to $x$ can be seen for the Kyriienko-inspired algorithm, despite
training being performed over Chebyshev nodes (different locations for the
nodes have been tested to confirm that the oscillations were still present).
Additionally, for this algorithm, the initial conditions are not met for the
derivatives, suggesting that the contribution of the loss function for the
initial values may conflict with the overall loss function. It suggest that
this quantum algorithm may lack flexibility to take into account all types of
boundary conditions, potentially due to constraints in the encoding or
variational ansatz.

The Hadamard-Lagrange algorithm reproduces the Lagrange interpolation and
allows the determination of the fitting polynomial with the smallest degree
over the given set of points. Hence, the function as well as their derivatives
are smooth and do not present much oscillations in comparison with the
Kyriienko-inspired algorithm. Moreover, the algorithm better adheres to the
boundary conditions. In general, the solutions obtained via our new Hadamard-
Lagrange algorithm appear to be of better quality compared to the solutions
obtained with the Kyriienko-inspired algorithm, especially at the boundaries,
and are in reasonably good agreement with the analytical solution as seen in
Figure 11. Furthermore, regarding the different loss contribution compared in
Figure 12, both algorithms converge toward a solution with a similar range of
DE loss, but the contribution of the boundary conditions to the loss function
considerably reduces the quality of the solutions obtained with the Kyriienko-
inspired VQAs.

Figure 11: Results comparison between the new Hadamard-Lagrange VQA and the
best outcomes of the Kyriienko-inspired VQA for both kinds of training sets of
points, for : a) the function $f$, b) the first derivative $\dot{f}$, c) the
second derivative $\ddot{f}$ and d) the DE loss contribution across 50
equispaced points in the interval, plotted on a logarithmic scale.  Figure 12:
DE and BC loss comparison between the new Hadamard-Lagrange VQA and the best
outcomes of the Kyriienko-inspired VQA for both kinds of training sets of
points, plotted on a logarithmic scale.

A direct comparison of the number of gates required between the solutions
obtained with the new VQA and the best outcomes of the Kyriienko-inspired VQA
reveals a significant difference, as illustrated in Figure 13. This comparison
positions the new Hadamard-Lagrange algorithms presented in this work as the
most precise and efficient for solving, at least, the damped mass spring
system problem. Finally, when considering both the quality of the solution and
its gate complexity, the optimal solution is derived from the new Hadamard-
Lagrange VQA, trained on a set of 7 training points of the first kind.

Figure 13: gate complexity comparison between the new Hadamard-Lagrange VQA
and the best outcomes of the Kyriienko-inspired VQA for both kinds of training
sets of points.

###  4.2 Poisson Equation

The Poisson equation is one of the most computationally intensive PDEs to
solve, requiring very expensive iterative methods. Previous work from Sato et
al. [32] solved the Poisson equation using a quantum approach for different
boundary conditions and a step function $s$ from $1/2^{n/2}$ to $-1/2^{n/2}$
for $n=5$, using a VQA based on the minimal potential energy, with $n+1=6$
qubits over $2^{n}=32$ points. In the present study, we solved the same
equation using the Hadamard-Lagrange algorithm with a simplified structure on
4 qubits trained over 3 points. The Poisson equation can be expressed as

$\frac{d^{2}f}{dx^{2}}+s=0,$ (22)

where $f$ is the state field and $s$ the given step function as a source term,
plotted in Figure 14. Different types of boundary conditions can be tested for
this equation, namely periodic, Neumann and Dirichlet boundary conditions. It
is a great opportunity to assess the versatility of the proposed new VQA. This
Poisson equation can be analytically solved by splitting the interval in half.
The solution over each half interval is a second degree polynomial where the
coefficients are real numbers and depend on the source function and the
boundary conditions. The definition of the analytical solutions (depending on
the boundary conditions) is further detailed in appendix C.

Figure 14: a) Source term $s$ used for the Poisson equation application. b)
Iterative evolution of the Hadamard-Lagrange solution for Neumann boundary
conditions.

####  4.2.1 Hadamard-Lagrange Algorithm vs Sato et al. Algorithm

The two quantum approaches being compared here for the solution of a Poisson
equation are quite different. The Sato et al. algorithm uses a discretised
approach, where the output is the unitary vector of amplitudes representing
the approximated solution function across the specified interval. In contrast,
the new algorithm using Lagrange polynomials encoding outputs the approximated
solution function, which can be computed for any point/node within the chosen
interval. Given its discretised nature, the Sato et al. algorithm encodes both
the source function $s$ and the algorithm output $\psi_{\theta}$ into quantum
states as vectors of amplitudes (Equation 5). However, this implementation
assumes the existence of an efficient unitary operator $\hat{\mathcal{U}}_{s}$
to prepare the quantum state $\ket{s}=\hat{\mathcal{U}}_{s}\ket{0}^{\otimes
n}$, which may limit its application to other source functions or other PDEs.

Figure 15: a) Iterative evolution of the DE loss $\mathcal{L}_{\theta}^{DE}$
obtained with the Hadamard-Lagrange algorithm; equal loss
$\mathcal{L}_{\theta}^{DE}$ to the Sato et al. [32] algorithm is obtained
after 300 iterations. b) Detailed comparison of the DE loss contribution over
the interval, $\mathcal{L}_{\theta,x}^{DE}$. Both plots are obtained with
Neumann boundary conditions, plotted on a logarithmic scale.

Our proposed Hadamard-Lagrange algorithm operates with continuous functions.
In cases where the source function is discontinuous, as in this particular
application, at least the second-order derivative will exhibit a
discontinuity. To address this challenge, the algorithm can be applied to sets
of sub-intervals. For example, the algorithm can be initially applied to the
first half of the interval, and the solution over the second half can be
determined by either symmetry or applying the algorithm again to the second
half. Here, the Hadamard-Lagrange algorithm was applied to half of the
interval and iterated until it achieved a precision comparable to the Sato et
al. solution, allowing a fair comparison under similar controlled conditions,
as shown in Figure 14.b. Figure 14.a illustrates the shape of the source term,
with a discontinuity. This figure illustrates how our VQA is eventually
converging to the analytical solution. It should also be noted that for our
VQA, the DE loss contribution over the half interval is constant, which is not
the case for the Sato et al. algorithm (dotted line in the Figure 15.b), for
which the loss is much larger at the boundaries of the computational domain.
The red line in Figure 15.b, corresponding to 300 iterations, is the loss
obtained with our new VQA, which was stopped when reaching the same averaged-
in-space loss as the Sato et al. solution. The green line, corresponding to
399, is the loss obtained when reaching the local minimum loss of the Sato et
al. solution. Note that the data obtained in this section are averaged over 5
tries to minimise the influence of the random initialisation.

The loss function of the proposed Hadamard-Lagrange algorithm accounts for
both the evaluation of the differential equation and the boundary conditions.
Figure 16 illustrates the contribution in the total loss function from the
differential equation and the constraints of the boundary conditions. It can
be seen that the contribution from the boundary conditions constraint is
different depending on the type of boundary conditions, with a significant
contribution for Neumann boundary conditions and a relatively small
contribution for periodic and Dirichlet boundary conditions. This is due to
the small difference in the analytical solutions depending on the boundary
conditions as explained in Appendix C.

Figure 16: Overall loss value of the Hadamard-Lagrange algorithm applied to
the different boundary conditions, plotted on a logarithmic scale.

For all boundary conditions (i.e., Periodic, Dirichlet, and Neumann), a direct
comparison to the analytical solution revealed an overall good agreement for
both VQA, even at the boundaries of the computational domain, as seen in
Figure 17. As a reminder, this is achieved with the use of a simplified
structure for the circuit of our new VQA. For each boundary condition, the
mean absolute error of the solution obtained using the Sato et al. algorithm
is 4, 3, and 6, larger than the error obtained for the proposed algorithm, for
periodic, Dirichlet and Neumann boundary conditions, respectively.

Figure 17: Comparison of the Hadamard-Lagrange solution with the one obtained
with the Sato et al. algorithm [32] for (a) periodic, (b) Dirichlet and (c)
Neumann boundary conditions

.

To further investigate the performance of both algorithms, it can be helpful
to look at the first- and second-order derivatives of the solution of the
Poisson equation. Focusing on the case with the Neumann boundary condition, it
can be seen in Figure 18, that the Sato et al. solution deviates from the
analytical solution at the boundary of the domain for the second order
derivative and near the discontinuity for the first order derivative. It is
important to note that the accuracy of the Sato et al. solution could be
improved by increasing the spatial discretisation, but this would dramatically
increase the size of the VQC and consequently the complexity.

Figure 18: Comparison of the first (a) and second (b) order derivatives of the
solution of the Poisson equation for the Neumann boundary condition, obtained
with the novel Hadamard-Lagrange algorithm and the Sato et al. algorithm [32].

####  4.2.2 Gate complexity estimation

Similarly to the previous application, the gate complexity of both VQAs has
been estimated by the number of circuits as well as the number of basic gates
required to complete the calculation. However, this estimation can only be
regarded as an order of magnitude estimation due to the fundamentally
different approaches of the two VQAs being compared. Moreover, for both
algorithms, a different classical optimiser had been used, which significantly
impacts the number of iterations and, consequently, the overall gate
complexity. As previously mentioned in section 2.1.3, the BFGS optimiser used
by Sato et al. [32] has better performances in training linear quantum solvers
over few qubits [38] than the Adam optimiser used for the Hadamard-Lagrange
algorithm.

For one iteration of the Sato et al. algorithm, the main circuit (Figure 4) is
evaluated with 3, 4, or 5 observables for periodic, Dirichlet, and Neumann
boundary conditions, respectively, to obtain the loss function and its
gradient with respect to the variational vector $\theta$, which contains 45
variational parameters. This requires numerous derivative circuits to compute
the gradient and a main circuit with a large variational ansatz. Additionally,
the main circuit (Figure 4) includes an encoding part and, in some cases, a
shift part depending on the observable.

Figure 19: Gate complexity estimation of the Hadamard-Lagrange algorithm,
averaged over five attempts per boundary conditions, compared with the Sato et
al. algorithm [32], plotted on a logarithmic scale.

In contrast, the proposed Hadamard-Lagrange algorithm requires circuits to
evaluate the loss function (Equation 7) and its gradient with respect to the
variational vector $\theta$ for each iteration. This includes the differential
equation loss over the control points $\\{x_{i}\\}$ and for the boundary
conditions at $x_{0}$. Here, the variational vector contains only 3
parameters. This structural difference leads to a significant disparity in
gate complexity. Although the Sato et al. algorithm [32] requires, on average,
half as many iterations to converge to the solution of the Poisson equation,
its circuit structure is much more complex, requiring 200 times more gates per
iteration compared to the Hadamard-Lagrange algorithm, as illustrated in
Figure 19. This is in line with the observation made for the first
application, for which the proposed Hadamard-Lagrange algorithm was already
very competitive in terms of gate complexity by comparison to the Kyriienko-
inspired algorithm.

Moreover, while being flexible due to its numerous variational parameters
$\theta$, the HEA [55] used in the Sato et al. algorithm faces complexities
arising from its multilayered structure and the sequence of rotational gates.
This large number of parametrised and non-parameterised gates also exacerbates
classical optimisation challenges, potentially leading to overfitting and
encountering barren plateaus. This complexity may explain why the Sato et al.
algorithm exhibits a fairly high gate complexity.

## 5\. Conclusions & Future work

PDEs are critical in various scientific fields such as structural engineering,
fluid dynamics, and financial modelling, but can be difficult to solve with
traditional computational methods due to their complexity. Recent advancements
in quantum computing spur researchers’ interest in designing quantum
algorithms for solving PDEs. This work introduced two distinct architectures
of a novel variational quantum algorithm (VQA) based on Lagrange polynomial
encoding alongside derivative quantum circuits with Hadamard test
differentiation to approximate PDE solutions. To demonstrate the potential of
this new VQA, it was applied to two well-known PDEs: the damped mass-spring
system from a given initial value and the Poisson equation under periodic,
Dirichlet, and Neumann boundary conditions. The results showed that the new
VQA with its two different circuit architectures has a reduced gate complexity
compared to existing variational quantum algorithms while providing similar or
better solutions.

Among the algorithms compared in this study, both the Kyriienko-inspired
algorithm and the proposed new Hadamard-Lagrange algorithm employ meshless
approaches and rely on polynomial fitting, demanding continuity in the
definition of the problem. Such feature is not required by the Sato et al.
algorithm. However, while the latter algorithm needs to discretise the problem
and hence can solve only linear or linearised PDEs, polynomial fitting
algorithms do not require any linearisation procedure, which makes them good
candidates to solve nonlinear PDEs.

Furthermore, these meshless approaches allow to distinguish up to three types
of points within the interval $I$: the training points, used to evaluate the
gradient of the loss function and optimise the variational parameters; the
evaluation points, used to evaluate the optimised function; and finally, the
encoding points, used in the quantum feature map of the Hadamard-Lagrange
algorithm. These different sets of points are independent, enabling the
training of VQCs on various sets or evolving sets with algorithm iterations.
This independence also allows to introduce additional losses at specific
points, different from the training set or not, such as boundary conditions.

On the other hand, where the Kyriienko-inspired algorithm encodes the variable
$x$ through a nonlinear set of polynomials, which causes some oscillations at
the extremity of the interval, encoding the variable with linear polynomials
such as with the proposed Hadamard-Lagrange algorithm may not be optimal in an
application to a nonlinear problem. In both cases the degree of the fitting
polynomial scales with the size of the circuit, which might, at some point,
limit their implementation to NISQ devices.

Based on the encouraging results showcased in this study and these comments,
we intend to further investigate the potential of our VQA and improve its
performance by looking at the following:

  * •

High dimensions and different PDEs: The work presented focused on 1D PDEs
(Damped Mass-Spring System (DMSS) and the Poisson equation), but most problems
of interest in engineering are in 2D and in 3D, and using a wide range of
PDEs. We therefore need to assess the performance of our VQA with more than
one spatial dimension and for a wide range of PDEs (advection-diffusion
equation, Burgers’s equation, and possibly the Euler equations).

  * •

Non-linear PDEs: Nonlinear PDEs and their corresponding dynamics are
frequently encountered in engineering. Due to the linear evolution of quantum
systems, constructing quantum algorithms for solving non-linear PDEs is quite
challenging. The next natural step is therefore to look at how our VQA could
deal with non-linear PDEs.

  * •

Tests on quantum computers: As discussed in the manuscript determining whether
to use the simplified or expanded circuit structure for our new VQA depends on
the quantum device being employed. The expanded circuit structure requires
more qubits but fewer gates. Hence, while it reduces noise from gate errors or
decoherence, this structure may be prone to errors from qubits or cross-talk.
Conversely, the simplified structure shows promise, particularly with the
development of efficient mitigation techniques and error gate correction. It
would therefore be of interest to investigate the performance of both
structures on quantum computers.

## Acknowledgments

This paper is dedicated to the memory of Professor Lorenzo Iannucci, whose
guidance and support shaped this project. The first author would like to
acknowledge the Department of Aeronautics, Imperial College London, for
supporting this work with a fully funded doctoral studentship. This research
was funded by the Engineering and Physical Sciences Research Council in the
United Kingdom, grant number EP/W032643/1.

## References

  * [1] V. Giovannetti, S. Lloyd, and L. MacCone, “Quantum random access memory,” Physical Review Letters, vol. 100, no. 16, pp. 1–4, 2008.
  * [2] O. D. Matteo, “A primer on quantum RAM.” 2020.
  * [3] O. D. Matteo, V. Gheorghiu, and M. Mosca, “Fault-Tolerant Resource Estimation of Quantum Random-Access Memories,” IEEE Transactions on Quantum Engineering, vol. 1, no. 1, pp. 1–13, 2021.
  * [4] D. Castelvecchi, “IBM releases first-ever 1,000-qubit quantum chip,” Nature, vol. 624, no. 7991, p. 238, 2023.
  * [5] J. M. Gambetta, J. M. Chow, and M. Steffen, “Building logical qubits in a superconducting quantum computing system,” npj Quantum Information, vol. 3, no. 1, pp. 0–1, 2017.
  * [6] P. W. Shor, “Algorithms for quantum computation: discrete logarithms and factoring,” in Proceedings 35th annual symposium on foundations of computer science, pp. 124–134, Ieee, 1994.
  * [7] L. K. Grover, “A fast quantum mechanical algorithm for database search,” in Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, pp. 212–219, 1996.
  * [8] G. Brassard, P. Høyer, M. Mosca, and A. Tapp, “Quantum amplitude amplification and estimation.” 5 2002.
  * [9] B. Kacewicz, “Almost optimal solution of initial-value problems by randomized and quantum algorithms,” Journal of Complexity, vol. 22, pp. 676–690, 10 2006.
  * [10] F. Gaitan, “Finding flows of a Navier–Stokes fluid through quantum computing,” npj Quantum Information, vol. 6, no. 1, 2020.
  * [11] F. Oz, R. K. Vuppala, K. Kara, and F. Gaitan, “Solving Burgers’ equation with quantum computing,” Quantum Information Processing, vol. 21, no. 1, pp. 1–13, 2022.
  * [12] F. Oz, O. San, and K. Kara, “An efficient quantum partial differential equation solver with chebyshev points,” Scientific Reports, vol. 13, no. 1, p. 7767, 2023.
  * [13] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” Physical Review Letters, vol. 103, no. 15, pp. 1–4, 2009.
  * [14] Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, “Quantum algorithm and circuit design solving the Poisson equation,” New Journal of Physics, vol. 15, no. 1, 2013.  013021\.
  * [15] P. C. Costa, S. Jordan, and A. Ostrander, “Quantum algorithm for simulating the wave equation,” Physical Review A, vol. 99, no. 1, 2019.  012323\.
  * [16] S. Wang, Z. Wang, W. Li, L. Fan, Z. Wei, and Y. Gu, “Quantum fast Poisson solver: the algorithm and complete and modular circuit design,” Quantum Information Processing, vol. 19, pp. 1–25, 2020.
  * [17] A. M. Childs, J.-P. Liu, and A. Ostrander, “High-precision quantum algorithms for partial differential equations,” Quantum, vol. 5, 2021.  574\.
  * [18] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, “Preconditioned quantum linear system algorithm,” Physical Review Letters, vol. 110, no. 25, 2013.  250504\.
  * [19] A. Montanaro and S. Pallister, “Quantum algorithms and the finite element method,” Physical Review A, vol. 93, p. 032324, 3 2016.
  * [20] D. W. Berry, “High-order quantum algorithm for solving linear differential equations,” Journal of Physics A: Mathematical and Theoretical, vol. 47, no. 10, 2014.
  * [21] D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, “Quantum Algorithm for Linear Differential Equations with Exponentially Improved Dependence on Precision,” Communications in Mathematical Physics, vol. 356, no. 3, pp. 1057–1081, 2017.
  * [22] A. M. Childs and J. P. Liu, “Quantum Spectral Methods for Differential Equations,” Communications in Mathematical Physics, vol. 375, no. 2, pp. 1427–1457, 2020.
  * [23] S. K. Leyton and T. J. Osborne, “A quantum algorithm to solve nonlinear differential equations.” 12 2008.
  * [24] M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, “Variational quantum algorithms for nonlinear problems,” Phys. Rev. A, vol. 101, p. 010301, 2020.
  * [25] A. Sarma, T. W. Watts, M. Moosa, Y. Liu, and P. L. McMahon, “Quantum variational solving of nonlinear and multi-dimensional partial differential equations,” arXiv preprint arXiv:2311.01531, 2023.
  * [26] A. Peruzzo, J. McClean, P. Shadbolt, M. H. Yung, X. Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications, vol. 5, no. May, 2014.
  * [27] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, “The theory of variational hybrid quantum-classical algorithms,” New Journal of Physics, vol. 18, no. 2, 2016.
  * [28] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, The Variational Quantum Eigensolver: A review of methods and best practices, vol. 986.  Physics reports, 2022.
  * [29] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, “Quantum circuit learning,” Physical Review A, vol. 98, no. 3, pp. 1–3, 2018.
  * [30] O. Kyriienko, A. E. Paine, and V. E. Elfving, “Solving nonlinear differential equations with differentiable quantum circuits,” Physical Review A, vol. 103, pp. 1–22, 5 2021.
  * [31] A. J. Pool, A. D. Somoza, M. Lubasch, and B. Horstmann, “Solving Partial Differential Equations using a Quantum Computer,” Proceedings - 2022 IEEE International Conference on Quantum Computing and Engineering, QCE 2022, pp. 864–866, 2022.
  * [32] Y. Sato, R. Kondo, S. Koide, H. Takamatsu, and N. Imoto, “Variational quantum algorithm based on the minimum potential energy for solving the Poisson equation,” Physical Review A, vol. 104, pp. 1–4, 6 2021.
  * [33] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits,” Nature Communications, vol. 12, no. 1, pp. 1–39, 2021.
  * [34] V. Heyraud, Z. Li, K. Donatella, A. Le Boité, and C. Ciuti, “Efficient Estimation of Trainability for Variational Quantum Circuits,” PRX Quantum, vol. 4, no. 4, p. 1, 2023.
  * [35] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, “Variational quantum algorithms,” Nature Reviews Physics, vol. 3, no. 9, pp. 625–644, 2021.
  * [36] L. Bittel and M. Kliesch, “Training Variational Quantum Algorithms Is NP-Hard,” Physical Review Letters, vol. 127, no. 12, 2021.
  * [37] D. P. Kingma and J. L. Ba, “Adam: A method for stochastic optimization,” 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings, pp. 1–15, 2015.
  * [38] A. Pellow-Jarman, I. Sinayskiy, A. Pillay, and F. Petruccione, “A comparison of various classical optimizers for a variational quantum linear solver,” Quantum Information Processing, vol. 20, no. 6, pp. 1–14, 2021.
  * [39] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subaşı, L. Cincio, and P. J. Coles, “Variational Quantum Linear Solver,” Quantum, vol. 7, 2023.
  * [40] J. Spall, “Multivariate stochastic approximation using a simultaneous perturbation gradient approximation,” IEEE Transactions on Automatic Control, vol. 37, pp. 332–341, 3 1992.
  * [41] M. J. D. Powell, “A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation,” in Advances in Optimization and Numerical Analysis, pp. 51–67, Dordrecht: Springer Netherlands, 1994.
  * [42] M. Wilson, R. Stromswold, F. Wudarski, S. Hadfield, N. M. Tubman, and E. G. Rieffel, “Optimizing quantum heuristics with meta-learning,” Quantum Machine Intelligence, vol. 3, no. 1, pp. 1–13, 2021.
  * [43] J. M. Kübler, A. Arrasmith, L. Cincio, and P. J. Coles, “An Adaptive Optimizer for Measurement-Frugal Variational Algorithms,” Quantum, vol. 4, p. 263, 5 2020.
  * [44] R. Sweke, F. Wilde, J. J. Meyer, M. Schuld, P. K. Faehrmann, B. Meynard-Piganeau, J. Eisert, P. K. Fährmann, B. Meynard-Piganeau, and J. Eisert, “Stochastic gradient descent for hybrid quantum-classical optimization,” Quantum, vol. 4, p. 314, 8 2020.
  * [45] K. M. Nakanishi, K. Fujii, and S. Todo, “Sequential minimal optimization for quantum-classical hybrid algorithms,” Physical Review Research, vol. 2, p. 043158, 10 2020.
  * [46] J. Stokes, J. Izaac, N. Killoran, and G. Carleo, “Quantum Natural Gradient,” Quantum, vol. 4, p. 269, 5 2020.
  * [47] B. Koczor and S. C. Benjamin, “Quantum natural gradient generalized to noisy and nonunitary circuits,” Physical Review A, vol. 106, p. 062416, 12 2022.
  * [48] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, “Barren plateaus in quantum neural network training landscapes,” Nature Communications, vol. 9, no. 1, pp. 1–6, 2018.
  * [49] E. Grant, L. Wossnig, M. Ostaszewski, and M. Benedetti, “An initialization strategy for addressing barren plateaus in parametrized quantum circuits,” Quantum, vol. 3, p. 214, 12 2019.
  * [50] H.-Y. Liu, T.-P. Sun, Y.-C. Wu, Y.-J. Han, and G.-P. Guo, “Mitigating barren plateaus with transfer-learning-inspired parameter initializations,” New Journal of Physics, vol. 25, p. 013039, 1 2023.
  * [51] L. Friedrich and J. Maziero, “Avoiding barren plateaus with classical deep neural networks,” Physical Review A, vol. 106, p. 042433, 10 2022.
  * [52] C.-Y. Park and N. Killoran, “Hamiltonian variational ansatz without barren plateaus,” Quantum, vol. 8, p. 1239, 2 2024.
  * [53] J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, “Generalized Unitary Coupled Cluster Wave functions for Quantum Computation,” Journal of Chemical Theory and Computation, vol. 15, pp. 311–324, 1 2019.
  * [54] M. Benedetti, E. Lloyd, S. Sack, and M. Fiorentini, “Parameterized quantum circuits as machine learning models,” Quantum Science and Technology, vol. 4, p. 043001, 11 2019.
  * [55] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature, vol. 549, pp. 242–246, 9 2017.
  * [56] G. Dahlquist and A. Bjork, “Equidistant interpolation and the Runge phenomenon,” in Numerical Methods, pp. 101–103, Prentice Hall, 1974.

##  Appendix A Hardware Efficient Ansatz Simplification

The Hardware Efficient Ansatz [55] was initially described as a quantum
circuit comprising layers of a sequence of rotational gates
$\hat{R}_{z}-\hat{R}_{x}-\hat{R}_{z}$ and CNOT gates, with independent angle
parameters $\theta$. However, this general structure can be simplified to
layers consisting of a single rotational gate around the x-axis $\hat{R}_{x}$
and CNOT gates, in the case of a projective measurement onto the Z-axis.

The magnetisation or projective measurement onto the Z-axis of a given qubit
$\ket{\varphi}$ can be defined as the difference in the probabilities of
measuring the qubit in the computational basis states:

$\displaystyle\braket{\hat{Z}}$
$\displaystyle=\bra{\varphi}\hat{Z}\ket{\varphi}$ (23)
$\displaystyle=P(\ket{0})-P(\ket{1})$
$\displaystyle=|\braket{0}{\varphi}|^{2}-|\braket{1}{\varphi}|^{2}.$

Since the $\hat{R}_{z}$ gate only affects the phase factor of the qubit state,
the probabilities of measuring the qubit in the computational basis states
remain unchanged.

$\ket{\varphi}=\begin{pmatrix}\alpha\\\
\beta\end{pmatrix};R_{z}(\phi)=\begin{pmatrix}e^{-i\frac{\phi}{2}}&0\\\
0&e^{i\frac{\phi}{2}}\end{pmatrix},$ (24)

$R_{z}(\phi)\ket{\varphi}=\begin{pmatrix}\alpha\ e^{-i\frac{\phi}{2}}\\\
\beta\ e^{i\frac{\phi}{2}}\end{pmatrix},$ (25)

$\displaystyle\bra{\varphi}R_{z}(\phi)^{\dagger}\ \hat{Z}\
R_{z}(\phi)\ket{\varphi}$
$\displaystyle=|\braket{0}{R_{z}(\phi)\varphi}|^{2}-|\braket{1}{R_{z}(\phi)\varphi}|^{2}$
(26) $\displaystyle=|\alpha|^{2}-|\beta|^{2}$
$\displaystyle=|\braket{0}{\varphi}|^{2}-|\braket{1}{\varphi}|^{2}$
$\displaystyle=\bra{\varphi}\hat{Z}\ket{\varphi}.$

##  Appendix B Estimation of the gate complexity

In this appendix, the methodology used to determine the number of circuits and
basic quantum gates required for one iteration is detailed. The basic gates
considered are Pauli and Clifford gates.

###  B.1 Kyriienko-Inspired Algorithm

Regarding the quantum solver inspired by the work of Kyriienko et al. [30],
one iteration involves - as detailed in the Background and Methods Sections -
circuits for the functions $f_{\theta}$, its derivatives
$\displaystyle\frac{df_{\theta}}{dx}\text{ ; }\frac{d^{2}f_{\theta}}{dx^{2}}$,
as well as their respective gradients
$\displaystyle\nabla_{\theta}(f_{\theta})\text{ ;
}\nabla_{\theta}(\frac{df_{\theta}}{dx})\text{ ;
}\nabla_{\theta}(\frac{d^{2}f_{\theta}}{dx^{2}})$ for all controlled points
and boundary conditions. In this algorithm, the gradients are obtained using
the parameter shift rule method, which requires two circuits for a partial
derivative. As a result, the number of circuits per iteration is defined as
follows:

$N_{\text{Circuits/iteration}}=N_{\text{Points}}(1+2N_{\text{Parameters}})\left(N_{\text{Circuits}}(f_{\theta})+N_{\text{Circuits}}(\frac{df_{\theta}}{dx})+N_{\text{Circuits}}(\frac{d^{2}f_{\theta}}{dx^{2}})\right),$
(27)

where $N_{\text{Circuit}}$ is zero if the respective component does not appear
in the differential equation. In the case of the damped-mass-spring-system
problem solved in Results section (4.1), the differential equation is composed
by all the components mentioned above.

$\begin{cases}N_{\text{Points}}=N_{\text{Controlled points}}+1,\\\
N_{\text{Parameters}}=N_{\text{qubits }}N_{\text{layers }},\\\
N_{\text{Circuits}}(f_{\theta})=1,\\\
N_{\text{Circuits}}(\frac{df_{\theta}}{dx})=2N_{\text{qubits }},\\\
N_{\text{Circuits}}(\frac{d^{2}f_{\theta}}{dx^{2}})=4N_{\text{qubits
}}^{2}.\end{cases}$ (28)

The circuits are constructed across 5 qubits with a variational ansatz depth
of 2. Moreover, the algorithm employs a quantum circuit structure based on
basic quantum gates such as CNOT, RY and RX quantum gates. The count of gates
per circuit can be expressed as follows:

$N_{\text{Basic gates/circuit}}=N_{\text{qubits }}+2N_{\text{Parameters}}.$
(29)

###  B.2 Hadamard-Lagrange Algorithm

As outlined in the Methods section and depicted in Figure 5, each iteration of
the Hadamard-Lagrange algorithm draws inspiration from the methodology
proposed by Kyriienko et al. [30]. Consequently, the number of circuits per
iteration can be defined by Equation 27. However, due to variations in
encoding and structure, the derivatives circuits with respect of $x$ are
obtained via the Hadamard test differentiation method, necessitating two times
less circuits per derivative.

$\begin{cases}N_{\text{Circuits}}(\frac{df_{\theta}}{dx})=N_{\text{Interpolation
points }},\\\
N_{\text{Circuits}}(\frac{d^{2}f_{\theta}}{dx^{2}})=N_{\text{Interpolation
points }}^{2}.\end{cases}$ (30)

In both the extended and reduced structures, the derivative circuit
necessitates an additional gate. Notably, the reduced structure exhibits a
slightly higher count of gates overall, despite being composed of half the
number of qubits.

For the extended structure:

$N_{\text{Basic gates/circuit}}(f_{\theta})=5N_{\text{Interpolation points
}}+2N_{\text{Parameters}}$ (31)

For the reduced structure:

$N_{\text{Basic gates/circuit}}(f_{\theta})=5N_{\text{Interpolation points
}}+2N_{\text{Parameters}}+\lfloor N_{\text{Interpolation points }}/2\rfloor.$
(32)

And in both cases,

$\begin{cases}N_{\text{Basic
gates/circuit}}(\frac{df_{\theta}}{dx})=N_{\text{Basic
gates/circuit}}(f_{\theta})+1\\\ N_{\text{Basic
gates/circuit}}(\frac{d^{2}f_{\theta}}{dx^{2}})=N_{\text{Basic
gates/circuit}}(f_{\theta})+2.\\\ \end{cases}$ (33)

Additionally, in the application of this algorithm, the number of qubits can
evolve throughout the algorithm depending on the number of interpolation
points. Consequently, the estimation of the number of circuits as well as the
number of gates should be considered step-by-step.

###  B.3 Sato et al. Algorithm

The Sato et al. algorithm [32] was designed to solve the Poisson equation
example described in the Results section. This algorithm employs a
fundamentally different approach and circuit structure. For each iteration,
evaluating the loss function—defined as the total potential energy of the
system requires 3, 4, or 5 circuits for periodic, Dirichlet, or Neumann
boundary conditions, respectively. Additionally, the derivative of the loss
function for each variational parameter must be evaluated. This is achieved
using the Hadamard test differentiation method, which requires only one
modified circuit per partial derivative. Thus, the number of circuits needed
for one iteration can be expressed as follows:

$N_{\text{Circuits/iteration}}=N_{\text{Observables
}}(1+N_{\text{Parameters}}),$ (34)

where $N_{\text{Observables }}=3,4\text{ or }5$ depending on the boundary
condition, and $N_{\text{Parameters}}=N_{\text{Layers}}\times
N_{\text{Parameters/layer}}+N_{\text{Encoding qubits}}$. For the Poisson
equation solved in the Results section, 5 qubits are used for encoding, the
number of layers of the ansatz is set to 5, and the number of parameters per
layer is 8, resulting in a total of 45 parameters. To account for the gate
complexity of the algorithm, the number of basic quantum gates (Pauli and
Clifford gates) per iteration is determined as follows:

$N_{\text{Basic gates/circuit}}=3+N_{\text{Basic
gates}}(\hat{\mathcal{U}}_{b})+N_{\text{Basic
gates}}(\hat{\mathcal{U}}_{\theta})+N_{\text{Shift Circuits}}N_{\text{Basic
gates}}(\hat{\mathcal{U}}_{\text{Shift}}),$ (35)

$\text{with, }\begin{cases}N_{\text{Basic
gates}}(\hat{\mathcal{U}}_{b})=N_{\text{Encoding qubits }}+1,\\\
N_{\text{Basic gates }}(\hat{\mathcal{U}}_{\text{Shift}})=1+N_{\text{Encoding
qubits }}N_{\text{Basic gates}}(MCP),\\\ N_{\text{Basic gates
}}(\hat{\mathcal{U}}_{\theta})=N_{\text{Parameters }}N_{\text{Basic gates
}}(CRY)+N_{\text{Layers }}N_{\text{Basic gates }}(MCP).\end{cases}$ (36)

The CRY quantum gate refers to the controlled rotational gate around the Y
axis, which can be decomposed into four basic quantum gates. The MCP quantum
gate, or the multi-controlled phase gate, is equivalent to a Toffoli quantum
gate in the case of two controlled qubits, which can be decomposed into 18
basic quantum gates.

##  Appendix C Analytical Solutions

###  C.1 Damped Mass Spring System

In this study, we investigate the dynamics of a damped mass-spring system
described by the second-order differential equation, referenced as Equation
18. The corresponding characteristic equation is formulated as:

$m\lambda^{2}+b\lambda+k=0.$ (37)

For the scenario explored in this paper, the discriminant $\Delta=b^{2}-4mk$
is negative, indicating an underdamped case. Consequently, the solution takes
the form:

$f(x)=e^{\alpha x}(C_{1}\cos(\beta x)+C_{2}\sin(\beta x)),$ (38)

where $\alpha$ and $\beta$ are the real and imaginary parts of the complex
conjugate roots $\lambda_{1,2}=\alpha\pm i\beta$. The coefficients $C_{1}$ and
$C_{2}$ are determined by the initial conditions $f(x_{0})=u_{0}$ and
$f^{\prime}(x_{0})=\dot{u}_{0}$.

###  C.2 Poisson Equation

The Poisson equation is defined in equation 22, with a given source term
defined as follow:

$s(x)=\begin{cases}+(1/2)^{n/2}&\text{if }x<\frac{a+b}{2},\\\
-(1/2)^{n/2}&\text{if }x>\frac{a+b}{2},\end{cases}$ (39)

with $n=5$, over the interval $[a,b]$. In the work of Sato et al. [32], this
Poisson equation is solved over 32 nodes, defining the interval as
$[\\![0,31]\\!]$. As mentionned in the Results section, the analytical
solution can be determined by splitting the interval in half. The solution
over each half is a second degree polynomial as expressed in the following
equation:

$f(x)=\begin{cases}Ax^{2}+B_{1}x+C_{1}&\text{if }x<\frac{a+b}{2},\\\
Ax^{2}+B_{2}x+C_{2}&\text{if }x>\frac{a+b}{2},\end{cases}$ (40)

with $A=-\frac{1}{2}^{\frac{n}{2}+1}$. $B_{i},\ C_{i}$ are real coefficients
which depend on the boundary conditions. In every cases, the continuity of the
middle point $x_{m}=\frac{a+b}{2}$ for the solution function and its first
derivative gives the following definitions:

$\lim_{x\to x_{m}^{+}}f(x)=\lim_{x\to x_{m}^{-}}f(x)=f(x_{m})$ (41)
$\lim_{x\to x_{m}^{+}}\frac{df}{dx}(x)=\lim_{x\to
x_{m}^{-}}\frac{df}{dx}(x)=\frac{df}{dx}(x_{m})$ (42)

Other conditions depend on the boundary. For the periodic boundary condition,
the solution follows $f(x+T)=f(x)$, with the period $T$ defined in Sato et al.
[32] work as $T=b-a+1=32$. The Dirichlet boundary condition consists in
defining the values at the extremities of the interval. In Sato et al. [32]
work, it is defined as $f(a-1)=f(b+1)=0$. Finally, the Neumann boundary
condition indicates the values of the first order derivative at the
extremities, which is defined here as $f^{\prime}(a)=f^{\prime}(b)=0$. All
considered, the analytical solution for each boundary equation are:

  * •

Periodic :

$f_{P}(x)=\begin{cases}A\left(x^{2}-\frac{1}{2}(3a+b-1)x+\frac{1}{2}(a-\frac{1}{2})(a+b+1)\right)&\text{if
}x<x_{m},\\\ 0&\text{if }x=x_{m},\\\
-A\left(x^{2}-\frac{1}{2}(3b+a+1)x+\frac{1}{2}(b+\frac{1}{2})(a+b+1)\right)&\text{if
}x<x_{m},\\\ \end{cases}$ (43)

  * •

Dirichlet :

$f_{D}(x)=\begin{cases}A\left(x^{2}-\frac{1}{2}(3a+b-2)x+\frac{1}{2}(a-1)(a+b)\right)&\text{if
}x<x_{m},\\\ 0&\text{if }x=x_{m},\\\
-A\left(x^{2}-\frac{1}{2}(a+3b+2)x+\frac{1}{2}(b+1)(a+b)\right)&\text{if
}x<x_{m},\\\ \end{cases}$ (44)

  * •

Neumann :

$f_{D}(x)=\begin{cases}A\left(x^{2}-2ax+\frac{1}{4}(3a^{2}+2ab-b^{2})\right)&\text{if
}x<x_{m},\\\ 0&\text{if }x=x_{m},\\\
-A\left(x^{2}-2bx+\frac{1}{4}(3b^{2}+2ab-a^{2})\right)&\text{if }x<x_{m}.\\\
\end{cases}$ (45)