Practical Quantum Computing: solving the wave equation using a quantum approach
Adrien Suau Gabriel Staffelbach Henri Calandra
quant-ph




adrien.suau@cerfacs.fr


  CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, France


  LIRMM, University of Montpellier, 161 rue Ada, 34095 Montpellier, France



  CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, France



  TOTAL SA, 2 Avenue de Vignancour, 64000 Pau, France




  In the last years, several quantum algorithms that try to address the problem of partial differential equation solving have been devised. On one side, “direct” quantum algorithms that aim at encoding the solution of the PDE by executing one large quantum circuit. On the other side, variational algorithms that approximate the solution of the PDE by executing several small quantum circuits and making profit of classical optimisers.
  In this work we propose an experimental study of the costs (in terms of gate number and execution time on a idealised hardware created from realistic gate data) associated with one of the “direct” quantum algorithm: the wave equation solver devised in [PCS. Costa, S. Jordan, A. Ostrander, Phys. Rev. A 99, 012323, 2019]. We show that our implementation of the quantum wave equation solver agrees with the theoretical big-O complexity of the algorithm. We also explain in great details the implementation steps and discuss some possibilities of improvements.
  Finally, our implementation proves experimentally that some PDE can be solved on a quantum computer, even if the direct quantum algorithm chosen will require error-corrected quantum chips, which are not believed to be available in the short-term.  



Practical Quantum Computing: solving the wave equation using a quantum approach
    Henri Calandra
    December 30, 2023
===============================================================================





§ INTRODUCTION



Quantum computing has drawn a lot of attention in the last few years, following the
successive announcements from several world-wide companies about the implementation of
quantum hardware with an increasing number of qubits or reduced error rates <cit.>.

Along with the hardware improvement, new quantum algorithms were discovered, yielding
potential quantum speed-up and applications in various fields such as quantum chemistry <cit.>, linear algebra <cit.> or optimisation <cit.>. Recent works even show that differential equations may be solved
by using a quantum computer <cit.>.
But despite the large number of algorithms available, it is hard to find an actual implementation
of a quantum differential equation solver, Hamiltonian simulation being the unique exception
by solving the time-dependant Schrödinger equation.


In this work, we present and analyse a quantum wave equation solver we implemented from scratch
according to the algorithm depicted in <cit.>. During the solver implementation, we had
to look for a Hamiltonian Simulation procedure. The implementations we found being too
restricted, we decided to implement our own Hamiltonian Simulation procedure, which will also be analysed.

To the best of our knowledge, this work is the first to analyse experimentally the characteristics
of a quantum PDE solver. Such a study has already been performed on the HHL algorithm, in <cit.>. 
We checked that the practical implementation agrees with the theoretical asymptotic complexities
on several quantities of interest such as the total gate count with respect to the number of discretisation points used or the precision, the number of qubits required versus the number of discretisation points used to approximate the solution or precision of the solution when compared to a classical finite-difference solver. Finally, we verified that the execution time of the generated quantum circuit on today's accessible quantum hardware was still following the theoretical asymptotic complexities devised for the total gate count. Quantum hardware data were extracted from IBM Q chips.

We show experimentally that it is possible to solve the 1-dimensional wave equation on a quantum computer
with a time-complexity that grows as N_d^3/2log(N_d)^2 where N_d is the number
of discretisation points used to approximate the solution. But even if the asymptotic scaling is better
than classical algorithms, we found out that the constants hidden in the big-O notation were huge
enough to make the solver less efficient than classical solvers for reasonable discretisation sizes.




































§ PROBLEM CONSIDERED


We consider a simplified version of the
wave equation on the 1-dimensional line [0, 1] where the propagation speed c
is constant and equal to 1. This equation can be written as


    ∂^2/∂ t^2ϕ(x,t) = ∂^2/∂ x^2ϕ(x,t).


Moreover, we only consider solving <ref> with the Dirichlet boundary conditions

    ∂/∂ xϕ(0, t) = ∂/∂ xϕ(1, t) = 0.


No assumption is made on initial speed ϕ(x, 0) and initial velocity
∂ϕ/∂ t(x, 0).

The resolution of this simplified wave equation on a quantum computer is an
appealing problem for the first implementation of a PDE solver for several
reasons. First, the wave equation is a well-known and intensively studied
problem for which a lot of theoretical results have been verified. Secondly,
even-though it is a relatively simple PDE, the wave equation can be used to
solve some interesting problems such as seismic imaging <cit.>. Finally, the theoretical
implementation of a quantum wave equation solver has already
been studied in <cit.>.

In this paper, we present the complete implementation of a
1-dimensional wave equation solver using quantum technologies based on  library. To the best of
our knowledge, this work is the first to consider the implementation of an
entire PDE solver that can run on a quantum computer. Specifically, we explain
all the implementation details of the solver from the mathematical theory to the
actual quantum circuit used. The characteristics of the solver are then discussed and
analysed, such as the estimated gate count and estimated execution time on
real quantum hardware. We show that the implementation follows the theoretical
asymptotic behaviours devised in <cit.>. Moreover, the wave equation
solver algorithm relies critically on an efficient implementation of a Hamiltonian
simulation algorithm, which we have also implemented and analysed thoroughly.





§ IMPLEMENTATION


The algorithm used to solve the wave equation is explained in <cit.>
and uses a Hamiltonian simulation procedure.  chose
the Hamiltonian simulation algorithm described in <cit.> for its
nearly optimal theoretical asymptotic behaviour. We privileged instead
the Hamiltonian simulation procedure explained in <cit.>
for its good experimental results based on <cit.> and its simpler
implementation (detailed in <ref>).

The code has been written using , a Python library
shipped with the Quantum Learning Machine (QLM), a package developed and
maintained by Atos. It has not been extensively optimized yet, which means that
there is still a large room for possible improvements.

All the circuits used in this paper have been generated with a subset of
's gate set:

    { H, X, R_y( θ), P_h( θ), CP_h( θ), CNOT, CCNOT }

and have then been translated to the gate set

    { U_1( λ), U_2( λ, ϕ), U_3( λ, ϕ, θ), CNOT }

for U_1, U_2 and U_3 defined in Equation (7) of <cit.> as follow:

    U(λ, ϕ, θ) =
      [           cos( θ/2)      -e^iλsin( θ/2);       e^iϕsin( θ/2) e^i(λ + ϕ)cos( θ/2) ]


    U_3(λ, ϕ, θ) = U(λ, ϕ, θ)


    U_2(λ, ϕ) = U(π/2, λ, ϕ)


    U_1(λ) = U(0, 0, λ)

 
  The target gate set presented in <ref>
  does not correspond to the physical gate set implemented by IBM hardware
  (see Equation (8) of <cit.>). This choice is justified by the
  fact that IBM only provides hardware characteristics such as gate times
  for the gate set of <ref> and not for the
  real hardware gate set. 






This implementation aims at validating in practice the theoretical asymptotic
complexities of Hamiltonian simulation algorithms and providing a
proof-of-concept showing that it is possible to solve a partial differential
equation on a quantum computer.



 §.§ Sparse Hamiltonian simulation algorithm








s-sparse matrix:
  A s-sparse matrix with s∈ℕ^* is a matrix that has at most s non-zero entries per row and per column


sparse matrix:
  A sparse matrix is a s-sparse matrix with s ∈log(N), N being the size of the matrix.


In the past years, a lot of algorithms have been devised to simulate the
effect of a Hamiltonian on a quantum state <cit.>.
Among all these algorithms, only few have already been implemented for
specific cases <cit.> but to the best of our
knowledge no implementation is currently capable of simulating a generic sparse
Hamiltonian. 

The domain of application of the already existing methods being too
narrow, we decided to implement our own generic sparse Hamiltonian simulation
procedure. We based our work on the product-formula approach described in <cit.>.
One advantage of this approach is that product-formula based algorithms have
already been thoroughly analysed both theoretically <cit.>
and practically <cit.>, and several
implementations are publicly available, though restricted to Hamiltonians that
can be decomposed as a sum of tensor products of Pauli matrices.
Moreover, <cit.>
provides a lot of implementation details that allowed us to go straight to the
development step.


Our implementation is capable of simulating an arbitrary sparse Hamiltonian
provided that it has already been decomposed into a sum of 1-sparse
Hermitian matrices with either only real or only complex entries, each described
by an oracle.
The implementation has been validated with several automated tests and a more
complex case involving the simulation of a 2-sparse Hamiltonian and described
in <ref>. Furthermore, it agrees perfectly
with the theoretical complexities devised in <cit.>
as studied and verified in <ref>.






 §.§ Quantum wave equation solver


Using the Hamiltonian simulation algorithm implementation, we successfully
implemented a 1-dimensional wave equation solver using the algorithm described
in <cit.> and explained in
<ref> and <ref>.

For the specific case considered (<ref> and
<ref>), 
solving the wave equation for a time T on a quantum computer boils down to
simulating a 2-sparse Hamiltonian for a time f(T), the function f being
thoroughly described in <cit.> and <ref>.
The constructed quantum circuit can then be
applied to a quantum state representing the initial position ψ(x, 0) and
velocity ∂ϕ/∂ t (x, 0), and will evolve this state
towards a quantum state representing the final position ϕ(x, T)
and velocity ∂ϕ/∂ t(x, T).

As for the Hamiltonian simulation procedure, the practical results we obtain from the implementation of the quantum wave equation solver seems to match the theoretical asymptotic complexities. See <ref> for an analysis of the theoretical asymptotic complexities.




§ RESULTS


Using a simulator instead of a real quantum computer has several advantages. In
terms of development process, a simulator allows the developer to perform several
actions that are not possible as-is on a quantum processor such as describing a
quantum gate with a unitary matrix instead of a sequence of hardware operations.
Another useful operation that is possible on a quantum simulator and
not currently achievable on a quantum processor is efficient generic state
preparation.

Our implementation uses only standard quantum gates and does not leverage any of
the simulator-only features such as quantum gates implemented from a unitary
matrix. In other words, both the Hamiltonian simulation procedure and the
quantum wave equation solver are “fully quantum” and are readily executable on
a quantum processor, provided that it has enough qubits. As a proof, and in
order to benchmark our implementation, we translated the generated quantum circuits
to IBM Q Melbourne gate-set (see <ref>). IBM Q Melbourne <cit.>
is a quantum chip with 14 usable qubits made available by IBM the 23th
of September, 2018.



  We chose IBM Q Melbourne mainly because, at the time of writing, it was the publicly accessible quantum chip with the larger number of qubits and so was deemed to be the closest to future quantum hardware. It is important to note that even if IBM Q Melbourne has 14 qubits, the quantum circuits constructed in this paper are not runnable because they require more qubits. Consequently, because of this hardware limitation, hardware topology has also been left apart of the study.


This allowed us to have an estimation of the number of hardware gates needed to either
solve the wave equation or simulate a specific Hamiltonian on this specific hardware.
Combining these numbers and the hardware gate execution time published in <cit.>, we were able to compute a rough approximation of
the time needed to solve the considered problem presented in
<ref> and <ref> on this specific hardware.



 §.§ Hamiltonian simulation


As explained in <ref>, the Hamiltonian simulation
algorithm implemented has been first devised in <cit.>.
A quick review of the algorithm along with implementation details can be
found in <ref>.
This Hamiltonian simulation procedure requires that the Hamiltonian matrix H
to simulate can be decomposed as 

    H = ∑_j=1^m H_j

where each H_j is an efficiently simulable Hermitian matrix.

In our benchmark, we simulated the Hamiltonian described in <ref>.
According to <cit.>,
real 1-sparse Hermitian matrices with only 1 or 0 entries can be
simulated with n gates and 2 calls to the oracle, n being the
number of qubits the Hamiltonian H acts on. The exact gate count can be found in <ref> in the row .

Let O_i be the gate complexity of
the oracle implementing the ith Hermitian matrix H_i of the decomposition in
<ref>, we end up with an asymptotic complexity
of n + O_i to simulate H_i. Once again, the exact gate count is decomposed in <ref>. 

Applying the Trotter-Suzuki product-formula of order k (see
Definition <ref> in <ref>
for the definition of the Trotter-Suzuki product-formula) on the quantum circuit
simulating the Hermitian matrices produces a circuit of size


    5^k ∑_i = 1^m (n + O_i).


This circuit should finally be repeated r times in order to achieve an error
of at most ϵ, with

    r ∈5^k m τ(m τ/ϵ)^1/2k,

and τ = t max_i || H_i ||, t being the time
for which we want to simulate the given Hamiltonian and
||·|| being the spectral norm <cit.>.

Merging <ref> and
<ref> gives us the complexity


    5^2k m τ(m τ/ϵ)^1/2k∑_i = 1^m (n + O_i).


This generic expression of the asymptotic complexity can be specialized to our
benchmark case. The number of gates needed to implement the
oracles is n^2 and the chosen decomposition contains m = 2
Hermitian matrices, each with a spectral norm of 1. Replacing the symbols in
<ref> and <ref>
results in the asymptotic gate complexity of

    5^k n^2

for the circuit
simulating e^-iHt/r and a number

    r ∈5^k t ( t/ϵ)^1/2k

of repetitions, which lead to a total gate complexity of

    5^2k n^2 t ( t/ϵ)^1/2k
      .


In order to check that our implementation follows this theoretical asymptotic
behaviour, we chose to let k = 1 and plotted the number of gates generated
versus the three parameters that have an impact on the number of gates: the
number of discretisation points N_d (<ref>),
the time of simulation t (<ref>) and the
precision ϵ (<ref>). The corresponding
asymptotic complexity should be

    n^2 t^3/2/√(ϵ)
       = 
        log_2( N_d )^2 t^3/2/√(ϵ)
      .


A small discrepency can be observed in <ref>: the theoretical asymptotic number of gates is log_2(N)^2 but the experimental values seem better fitted with an asymptotic behaviour of log_2(N)^7/4. This may be caused by the asymptotic regime not being reached yet.


















 §.§ Wave equation solver





The first characteristic of the wave equation solver that needs to be checked is its validity: is
the quantum wave equation solver capable of solving accurately the wave equation as described in
<ref> and <ref>?

To check the validity of the solver, we used  simulators and Atos QLM to simulate the
quantum program generated to solve the wave equation with different values for the number of discretisation
points N_d, for the physical time t and for the precision ϵ.
<ref> shows the classical solution
versus the quantum solution and the absolute error between the two solutions
for N_d = 32, t = 0.4 and ϵ = 10^-3.
The solution obtained by the quantum solver is nearly exactly the same as the classical solution
obtained with finite differences. The error between the two solutions is of the order of 10^-7,
which is 4 orders of magnitudes smaller than the error we asked for. 








Once the validity of our solver has been checked on multiple test cases, the next interesting
property we would like to verify is the asymptotic cost: does the implemented simulator seem
to agree with the theoretical asymptotic complexities derived from <cit.>
and <cit.>?

In our specific case, the Hamiltonian H to simulate can be decomposed in two
1-sparse Hermitian matrices, both of them having a spectral norm of 1. The exact decomposition can be found in <ref>. We chose to
let the product-formula order be equal to k=1 and reuse the asymptotic complexity
found in <ref> by changing the time
of simulation t by the time f(t):


    5^2k n^2 f(t) ( f(t)/ϵ)^1/2k
      .


Following the study performed in <cit.>,

    f(t) = t/δ x = t (N_d - 1)

where δ x is the distance between two discretisation points.
Moreover, it is possible to prove (see <ref>)
that 

    n = ⌈log_2(2N_d - 1) ⌉


Replacing f(t) and n in <ref> and
<ref> gives us a gate complexity of

    5^k log_2( N_d )^2

to construct a circuit simulating e^-iHt/r and a number of repetitions

    r ∈5^k t N_d ( t N_d/ϵ)^1/2k.




Merging the two expression results in a gate complexity of


    5^2k t N_d log_2(N_d)^2 (t N_d/ϵ)^1/2k
      .


Choosing the Totter-Suzuki formula order k = 1 gives us a final complexity of

    N_d^3/2log_2(N_d)^2 t^3/2/√(ϵ)

to solve the wave equation presented in <ref>.
This theoretical result is verified experimentally in
<ref>.







































§ DISCUSSION


In this work, we focus on the practical cost of implementing a 1-dimensional quantum wave equation solver on a quantum computer. We show that a quantum computer is able to solve partial differential equations by constructing and simulating the quantum circuits described. We also study the scaling of the solver with respect to several parameters of interest and show that the theoretical asymptotic bounds are mostly verified.

In future works, one can study the possibilities of circuit optimisation. It would also be interesting to implement Neumann boundary conditions instead of Dirichlet ones. A practical implementation including a non-constant propagation speed c has also been realised during the writing of this paper. The results were encouraging but were not judged mature enought to include them in the paper. Finally, future works might want to extend the wave equation solver to 2 dimensions or more. 













  The authors would like to thank Reims University, the ROMEO HPC center, Total, the CCRT and Atos for their support by giving us access to Atos quantum simulator.




§ SUPPLEMENTARY MATERIAL



The implementation of the quantum wave equation solver is available at <https://gitlab.com/cerfacs/qaths>.
The  tool is available at <https://gitlab.com/qcomputing/qprof/qprof>.




ACM-Reference-Format






§ PRODUCT-FORMULA IMPLEMENTATION DETAILS




 §.§ Hamiltonian simulation


Hamiltonian simulation is the problem of constructing a quantum circuit that will
evolve a quantum state according to a Hamiltonian matrix, following the
Schrödinger equation. In other words, Hamiltonian simulation algorithms generate a
quantum circuit performing the unitary transformation U such that || U - e^-iHt || < ϵ,
H being a given Hamiltonian matrix, t a time of evolution and ϵ a
precision with respect to ||·||, the spectral norm.

Several quantum algorithms have been developed in the last few years to solve the problem
of s-sparse Hamiltonian simulation <cit.>.
Among these algorithms we decided to implement the product-formula approach <cit.>,
for the reasons presented in <ref>.


The product formula algorithm has three main steps: decompose, simulate,
recompose. It works by first decomposing the s-sparse Hamiltonian matrix H that
should be simulated as a sum of Hermitian matrices H_j that are considered
easy to simulate

    H = ∑_j = 0^m-1 H_j.

The second step is then to simulate each H_j separately, i.e. to create
quantum circuits implementing e^-iH_jt for all the H_j in the decomposition
in <ref>. The last step uses the simulations computed in
step two to approximate e^-iHt.

The very first questions that should be answered before starting any
implementation of the product-formula algorithm are “What is an easy to
simulate matrix?” and “What kind of Hermitian matrices are easy to
simulate?”.



 §.§ Easy to simulate matrices



One of the most desirable properties for an “easy to simulate” matrix is the
possibility to simulate it exactly, i.e. to construct a quantum circuit that
will perfectly implement e^-iHt. This property becomes a requirement when one
wants rigorous bounds on the error of the final simulation. Another enviable
property of these matrices is that they can be simulated with
a low gate number and only a few calls to the matrix oracle. 


  
  A Hermitian matrix H can be qualified as “easy to simulate” if there exist
  an algorithm that takes as input a time t and the matrix H and outputs a quantum
  circuit C(H)_t such that
  
  
  * The quantum circuit C(H)_t implements exactly the unitary transformation
    e^-iHt, i.e.
    
    | | e^-iHt - C(H)_t | | = 0.

  
  * The algorithm only needs 1 calls to the oracle of H and
    log N additional gates, N being the dimension of the matrix H.
  





With this definition of an “easy to simulate” matrix, we can now search for
matrices or group of matrices that satisfy this definition.



  §.§.§ Multiples of the identity


The first and easiest matrices that fulfil the easy to simulate matrix requirements
are the multiples of the identity matrix {α I,  α∈ℝ}
with I the identity matrix. The quantum circuit to simulate this class of matrices can be
found in <cit.>.



  §.§.§ 1-sparse Hermitian matrices


A larger class of matrices that can be efficiently and exactly simulated are the
1-sparse, integer weighted, Hermitian matrices. Quantum circuits simulating exactly 1-sparse
matrices with integer weights can be found in <cit.>.


  Procedures simulating 1-sparse matrices with real (non-integers) weights are
  also described in the paper, but these matrices do not fall in the “easy to simulate”
  category because the procedures explained are exact only if all the
  matrix weights can be represented exactly with a fixed-point representation,
  which is not always verified.



  Multiples of identity matrices presented in 
  <ref> are a special case of 1-sparse
  matrices. The two classes have been separated because more efficient
  quantum circuits exists for α I matrices.




 §.§ Decomposition of H


Once the set of “easy to simulate” matrices has been established, the next
step of the algorithm is to decompose the s-sparse matrix H as a sum of
matrices in this set.

There are two possible ways of performing this decomposition, each one with its
advantages and drawbacks: applying a procedure computing the decomposition
automatically, or decompose the matrix H beforehand and provide
the decomposition to the algorithm.

The first solution, which is to automatically construct the oracles of
the H_j matrices from the oracle of the H matrix has been studied in <cit.>
and <cit.>. Thanks to this automatic decomposition procedure, we only need to implement one oracle. This simplicity comes at the cost of a higher gate count: each call to the automatically constructed oracles of the matrices H_j will require several calls to the oracle of H along with additional gates.

On the other hand, the second solution offers more control at the cost of less abstraction and more work. The decomposition of H is not automatically computed and should be performed beforehand. Once the matrix H has been decomposed as in <ref>, the oracles for the matrices H_j should be implemented.
This means that we should now implement m oracles instead of only 1 for the first solution. The main advantage of this method over the one using automatic-decomposition is that it gives us more control, a control that can be used to optimize even more the decomposition of <ref> (less H_j in the decomposition, H_j matrices that can be simulated more efficiently, …).










All the advantages and drawbacks weighted, we chose to implement the second option for
several reasons. First, the implementation of the automatic decomposition
procedure adds a non-negligible implementation complexity to the whole
Hamiltonian simulation procedure. Moreover, the automatic decomposition
procedure can be implemented afterwards and plugged effortlessly to the
non-automatic implementation. Finally, our use-case only required to simulate
a 2-sparse Hamiltonian that can be decomposed as the sum of two 1-sparse, easy to simulate,
Hermitian matrices, which makes the manual decomposition step manageable. 



 §.§ Simulation of the H_j


Once the matrix H has been decomposed following <ref>
with each H_j being an “easy to simulate” matrix, the simulation
of H_j becomes a straightforward application of the procedures described in
<ref>.

After this step, we have access to quantum circuits implementing e^-iH_jt for
j ∈[ 0, m-1 ] and t∈ℝ.



 §.§ Re-composition of the e^-iH_jt


The ultimate step of the algorithm is to approximate the desired evolution e^-iHt with the
evolutions e^-iH_jt. In the special case of mutually commuting H_j, this step is trivial
as it boils down to use the properties of the exponential function on matrices and write
e^iHt = e^i∑_jH_jt = ∏_j e^iH_jt. But in the more realistic case where the matrices
H_j do not commute, a more sophisticated method should be used to approximate the evolution
e^-iHt. To this end, we used the first-order Lie-Trotter-Suzuki product formula defined in
Definition <ref>.


  The Lie-Trotter-Suzuki product formula approximates
  
    exp( λ∑_j=0^m-1α_j H_j  )

  with 
  
    S_2(λ) = ∏_j=0^m-1 e^α_j H_j λ / 2∏_j=m-1^0 e^α_j H_j λ / 2

  and can be generalized recursively to higher-orders
  
    S_2k( λ) =  [ S_2k-2( p_k λ) ]^2 × S_2k-1( ( 1 - 4p_k ) λ) [ S_2k-2( p_k λ) ]^2

  with p_k = ( 4 - 4^1/(2k-1))^-1 for k> 1.
  Using this formula, we have the approximation
  
    e^λ H = [ S_2k( λ/n) ]^n + |λ|^2k+1/n^2k.

  


We used the Lie-Trotter-Suzuki product formula with λ = -it to approximate
the operator e^-iHt up to an error of ϵ∈t^2k+1/n^2k.



§ HERMITIAN MATRIX CONSTRUCTION AND DECOMPOSITION


One of the main challenge in implementing a quantum wave equation solver lies in
the construction and implementation of the needed oracles. This appendix
describes the first step of the implementation process: the construction and
decomposition of the Hamiltonian matrix that will be simulated using the
Hamiltonian simulation procedure introduced in <ref>.

This appendix follows the analysis performed in <cit.> and adds
details and observations that will be refereed to in
<ref> when dealing with the actual oracle implementation.




 §.§ Hamiltonian matrix description


In order to devise the Hamiltonian matrix that should be simulated to solve the
wave equation, the first step is to discretise <ref>
with respect to space. Such a discretisation can be seen as a graph G_δ
  x whose vertices are the discretisation points and with edges between
nearest neighbour vertices. The graph G_δ x is depicted in 
<ref>.




The graph Laplacian of G_δ x, defined as

    L(G_δ x)_i,j :=
      (v_i)    if  i = j 
    
        -1    if  (i ≠ j)  (v_i  v_j) 
    
        0    otherwise

can then be used to approximate the differential operator
∂^2/∂ x^2. By using the discretisation approximation

    ∂^2 ϕ/∂ x^2(i δ x, t) ≈ϕ_i-1,t - 2 ϕ_i,t + ϕ_i+1,t/δ x^2

with ϕ_i,t = ϕ(iδ x, t), 
and approximating ϕ(x, t) with a vector ϕ = [ ϕ_i,t]_0 ⩽ i < N_d, the matrix

    A = - 1/δ x^2 L(G_δ x)

approximates the second derivative of ϕ when δ x → 0 as

    [ Aϕ]_i = ϕ((i-1)δ x, t) - 2 ϕ(iδ x, t) + ϕ((i+1)δ x, t)/δ x^2.


The approximation in <ref> is then used in
<ref> to approximate the spatial derivative operator: 

    ∂^2/∂ t^2ϕ = -1/δ x^2 L(G_δ x) ϕ.


Based on this formula, <cit.> shows that simulating

    H = 
      [   0   B; B^†   0;     ]

with

    BB^† = L(G_δ x)

constructs a quantum circuit that will evolve a part of
the quantum state it is applied on according to the discretised wave equation in
<ref>.

A matrix B satisfying <ref> can be obtained directly from the
graph G_δ x representing the discretisation. The algorithm to construct the
matrix B can be decomposed in three steps. First, the vertices (discretisation
points) should be arbitrarily ordered by assigning them a unique index in [0, N_d-1].
Then, each edge of the graph is arbitrarily oriented and indexed with indices in [0,
N_d-2]. Finally, B is computed with the following definition 

    B_ij =
      {[                                         1 if edge  j  is a self-loop of vertex  i ,;                                         1     if edge  j  has vertex  i  as source,;                                        -1       if edge  j  has vertex  i  as sink,;                                         0                                 otherwise;                                           ]..


Note that edges' orientation and vertices/edges ordering is completely arbitrary.
Changing either the edges orientation on one of the orderings will change the
matrix B but will not affect BB^† which should be equal to L(G_δ x). This
freedom in the ordering and orientation choices takes a crucial importance in
the oracle implementation as it allows us to pick the ordering/orientation that
will produce an easy-to-implement matrix B.



 §.§ Dirichlet boundary conditions


Fixing boundary conditions is a requirement for most of the partial differential
equations to admit a unique well-defined solution. There exist several boundary
conditions such as Neumann, Dirichlet, Robin or Cauchy ones. For simplicity, we
restricted ourselves to the study of <ref> with
Dirichlet boundary condition of <ref>.

In the case of Dirichlet boundary conditions on the 1-dimensional line [0,
1], the two boundary nodes at x = 0 and x = 1 can be ignored as their value
is always equal to 0. Moreover, <cit.> shows that the graph
G_δ x^D representing the discretisation with Dirichlet boundary
conditions of <ref> is simply G_δ x with
self-loops on the two outer nodes (i.e. the ones indexed 1 and N_d-2 as 0
and N_d-1 are ignored). G_δ x^D is depicted in 
<ref>. The algorithm to construct the matrix B remain
the same as explained in <ref>.





 §.§ Matrices construction


All the pieces are now in place to start building the matrix
B_d∈ℝ^ (N_c-1) × N_c. Using the definition of the
matrix B written in <ref>
and the graph G_δ x^D depicted in 
<ref> we end up with 


    B_d =
      [  1  1  0  ⋯  0;  0 -1  1  ⋱  ⋮;  ⋮  ⋱  ⋱  ⋱  0;  0  ⋯  0 -1  1;    ]


We can easily check that B_dB_d^† is equal to the
well-known discretisation matrix

    B_dB_d^† =
      [  2 -1  0  ⋯  0; -1  2  ⋱  ⋱  ⋮;  0  ⋱  ⋱  ⋱  0;  ⋮  ⋱  ⋱  2 -1;  0  ⋯  0 -1  2;    ],

which validate the method of construction of B_d.

Computing H_d, the Hamiltonian matrix that should be simulated to
evolve the quantum state according to the wave equation in
<ref> with Dirichlet boundary conditions, is now
straightforward. Using <ref>, we directly obtain


    H_d = 1/δ x[  0  ⋯  ⋯  0  1  1  0  ⋯  0;  ⋮        ⋮  0 -1  1  ⋱  ⋮;  ⋮        ⋮  ⋮  ⋱  ⋱  ⋱  0;  0  ⋯  ⋯  0  0  ⋯  0 -1  1;  1  0  ⋯  0  0  ⋯  ⋯  ⋯  0;  1 -1  ⋱  ⋮  ⋮           ⋮;  0  1  ⋱  0  ⋮           ⋮;  ⋮  ⋱  ⋱ -1  ⋮           ⋮;  0  ⋯  0  1  0  ⋯  ⋯  ⋯  0;    ]





As explained in <ref>, the Hamiltonian simulation algorithm
implemented requires that the Hamiltonian to simulate is split as a sum
of 1-sparse Hermitian matrices. There are a lot of valid decompositions for
the matrix H_d and we are free to choose the decomposition that will
simplify the most the oracle implementation or reduce the gate complexity.

We made the choice to decompose B_d as two 1-sparse matrices and
then reflect this decomposition on H_d. Let B_1 and B_-1
defined as

    B_1 =
      [ 0 1 0 ⋯ 0; ⋮ ⋱ 1 ⋱ ⋮; ⋮   ⋱ ⋱ 0; 0 ⋯ ⋯ 0 1;   ]


    B_-1 =
      [  1  0  ⋯  ⋯  0;  0 -1  ⋱     ⋮;  ⋮  ⋱  ⋱  ⋱  ⋮;  0  ⋯  0 -1  0;    ]

we have B_d = B_1 + B_-1. Let also

    H_1 = 1/δ x[     0   B_1; B_1^†     0;       ]
      ,  H_-1 = 1/δ x[      0   B_-1; B_-1^†      0;        ],

it is easy to see that H_d = H_1 + H_-1
and that both H_1 and H_-1 are 1-sparse Hermitian
matrices.

For convenience, we also define

    H_1 = 
      [     0   B_1; B_1^†     0;       ]
      ,  
      H_-1 =
      [      0   B_-1; B_-1^†      0;        ],

and H_d = H_1 + H_-1, the H_1, H_-1 and H_d matrices rescaled to contain only integer weights. These matrices have the interesting property that simulating  H_d (resp. H_1, H_-1) for a time t is equivalent to simulating H_d (resp. H_1, H_-1) for a time t/δ x. This property will be used in the following sections as it offers us the opportunity to simulate the integer-weighted matrices  H_d, H_1 and H_-1 instead of the real-weighted ones H_d, H_1 and H_-1.

Note also that a lower bound of the number of qubits needed to solve the wave equation for N_d
discretisation points can be computed from the dimensions of H_d. As
the non-empty upper-left block of matrix H_d is of dimension ( 2N_d - 1
) ×( 2N_d - 1 ), we need at least 

    ⌈log_2(2 N_d - 1) ⌉

qubits to simulate it. This estimation does not take into account ancilla qubits that may be needed to implement the oracles.




§ ORACLE CONSTRUCTION


Oracles can be seen as the interface between a quantum procedure and real-world data.
Their purpose is to encode classical data such that a quantum algorithm can process
it efficiently.



 §.§ Oracle interface


In order to work as a bridge between the classical and the quantum worlds and to be
used by the quantum algorithm, a clear interface for the oracle should be
established. 

We chose to use the interface described in <cit.>
with slight modifications improving the arity of the oracle for our specific case
of 1-sparse matrices.

More precisely, our oracles O implement the following interface

    O |x_0⟩_x |0⟩_m |0⟩_v |0⟩_s = |x_0⟩_x |m(x_0)⟩_m |v(x_0)⟩_v |s(x_0)⟩_s

with |x_0⟩_x encoding a row index as a unsigned integer, m(x) the function that
returns the column index of the only non-zero element in row x,
v(x) = | w(x) | the absolute value of the weight w(x) of the first (and only)
non-zero element in row x and

    s(x) = {
          0     if  w(x) ⩾ 0 
    
          1     else .

the sign of the first non-zero entry in row x.

Note that the quantum register are labeled with their respective usage: x for the index of of the row considered, m for the index of the column considered, v for the value of the element at (row, index) and s for the sign of the element at (row, index). A fifth label “a” is used along the paper to label a register used as an ancilla. 

The interface of the oracle O can also be obtained with 3 separate oracles
that will each take care of computing one output:

    M |x_0⟩_x |0⟩_m = |x_0⟩_x |m(x)⟩_m


    V |x_0⟩_x |0⟩_v = |x_0⟩_x |v(x)⟩_v


    S |x_0⟩_x |0⟩_s = |x_0⟩_x |s(x)⟩_s




 §.§ Optimisation of M and S



  The simulation algorithms provided by <cit.>
  have the interesting property that if the oracle V encodes a weight of
  zero for some inputs (i.e. v(x) = 0 for some x) then the outputs of
  oracles M and S are ignored for those inputs.



  
  The circuit simulating a 1-sparse m-bit-integer weighted Hamiltonian H
  depicted in <ref> is taken from <cit.>.


  
  
  In our special case of 1-bit weights (i.e. m = 1), the third quantum gate
  e^-i Z ⊗ F_m t can be written as
  
    e^-i Z ⊗ F_m t   = e^-i Z ⊗ F_1 t
       = exp[ -i 
            [  F_1    0;    0 -F_1;      ]
            t ] 
       = [ e^-iF_1t        0;        0  e^iF_1t;          ]
       =
          [     1     0     0     0;     0 e^-it     0     0;     0     0     1     0;     0     0     0 e^-it;       ].


  where 

  
    F_m =
        [     0     0     ⋯     ⋯     ⋯     ⋯     0;     0     1     ⋱                       ⋮;     ⋮     ⋱     2     ⋱                 ⋮;     ⋮           ⋱     3     ⋱           ⋮;     ⋮                 ⋱     4     ⋱     ⋮;     ⋮                       ⋱     ⋱     0;     0     ⋯     ⋯     ⋯     ⋯     0 2^m-1;       ].

  
  It follows from the matrix notation that if the second qubit e^-iZ⊗ F_1 t is
  applied on is in the state |0⟩, the gate e^-iZ⊗ F_1 t is the identity
  transformation, i.e. the unitary operation e^-iZ⊗ F_1 t sends |00⟩ (resp. |10⟩)
  to |00⟩ (resp. |10⟩).
  This means that if the oracle O does not set the last qubit to |1⟩ (i.e. the oracle encodes a weight of
  0 for the xth row of H), the quantum circuit
  depicted in <ref> can be simplified up to an identity
  transformation as the effects of O (resp. A) are reverted by O^† (resp. A^†).

  Rephrasing, if the xth row of matrix H has no non-zero entries, the effects of the oracle O is ignored, which implies that the effects of the oracles M and S that compose O are also ignored.


Using the result of Claim <ref>, we are free to
implement any transformation that best suits us for the set of inputs |x⟩ such that
the xth row of the considered Hermitian matrix (H_1 or H_-1) has no non-zero
elements as long as the oracle V implements the right transformation.

To illustrate clearly the implemented transformations we chose to encode with
M and S, the next sections will re-write the matrices H_1 and H_-1 according
to <ref> but with one 0 or -0 in each empty
row. A 0 entry at position (i,j) in the matrix means that the row i was
empty, the oracle M will map |i⟩_x to |j⟩_m and the oracle S will encode a
positive sign, i.e. |0⟩_s. The same reasoning applies for -0 entries, except
that the encoded sign is now negative, i.e. |1⟩_s. 

The following sections will explain step by step the construction of each
of the three oracles M, V and S, both for the matrix H_1 (M_1, V_1 and S_1)
and the matrix H_-1 (M_-1, V_-1 and S_-1).



 §.§ About arithmetic and logic quantum gates


Implementing the oracles M, V and S for the matrices H_1 and H_-1 requires several
arithmetic and logic quantum gates such as ,  or .

All these gates have been implemented prior to the oracle implementation and the implementation
steps are detailed in this section.



  §.§.§ The  gate


The  gate is easily implemented using only  and  (or Toffoli) gates.
The implementation used is depicted in <ref> and uses
the famous Boole algebra formula linking ,  and :
x ∨ y = ( x ∧ y).






  §.§.§ The  and  gates


Most of the research papers presenting an implementation of the  or 
gates only consider the case where the two numbers to add or subtract are stored in
quantum registers.

In our case, the oracles implementation requires an adder and subtractor that can add or
subtract to a quantum register a quantity known when the quantum circuit is generated, i.e. not necessarily encoded on a quantum state.


  Implementing a subtractor is trivial once an adder procedure is available.


  A subtractor can be implemented from a generic adder by using the identity
  
    a - b = (a' + b)'

  where ' denotes the bitwise complementation.

  The circuit resulting of the application of this identity is depicted in
  <ref> and only requires one call
  to the adder and 2n additional gates, n being the number of qubits used to
  represent one of the operands.

  




  Following Claim <ref> we will restrict the study to
  implementing an adder. Implementing a subtractor is trivial and cheap in term
  of additional quantum gates used once an adder is available.


Generation-time value
  A generation-time value is a value that is known by the programmer when generating the
  quantum circuit. Knowing a value at generation-time may allow to optimise even further
  the generated quantum circuit. The closest analog in classical programming would be
  C-like macros or recent C++ constexpr expressions.


The easiest solution to overcome the problem caused by the non-compatible input formats between
our problem (with a generation-time value) and the existing adders (with two values encoded on quantum registers)
is to encode the quantity known at generation-time into ancillary qubits and then use the
regular adder algorithms to add to a quantum register the value encoded in a second quantum register.
Even if this solution is trivial to implement, it has the huge downside of requiring log_2 b
additional ancillary qubits to temporarily store the generation-time value b.

Another answer to the problem would be to adapt a quantum adder originally devised
to add two quantum registers to a quantum adder capable of adding a constant value
to a quantum register. Several adders <cit.>
have been studied to check if they can be modified to allow a generation-time input,
i.e. if it possible to remove completely the quantum register storing the right-hand-side
(or left-hand-side) of the addition.








The task of removing the quantum register storing one of the operands
appears to be challenging for adders based on classical arithmetic like <cit.> but trivial for Draper's quantum adder
introduced in <cit.>.





  Draper's quantum adder can be adapted into an efficient adder that takes as right-hand side input
  a unsigned “generation-time” integer value and add this value to a sufficiently large quantum
  register encoding another unsigned integer.



  The original Drapper's adder as introduced in <cit.> is illustrated in
  <ref>.
  
  

  The only quantum gates using the quantum register |b⟩ are the controlled-phase gates.
  Moreover, they only use the qubits of the right-hand-side register |b⟩ as controls.

  In the case of a constant value of b known at generation time, we can replace each
  controlled-phase gate by either a phase gate if the corresponding bit of b is 1 or
  by an identity gate (or a “no-op” gate) if the bit of b is 0. Once this transformation
  has been performed, the quantum register |b⟩ is no longer used and can be safely removed
  from the circuit. 


The final quantum  gate implementation is depicted in <ref>,
requires n^2 gates and has a depth of n. Following <cit.>,
the asymptotic gate count can be improved to nlog(n) by removing the rotation with an angle below
a given threshold that depend on hardware noise.








  §.§.§ The  gate


For the same reasons exposed in the adder implementation in <ref>,
the  gate cannot be implemented using the arithmetic comparator presented
in <cit.> because removing the right-hand side qubits seems to be
a challenging task.

Instead, we use the idea from <cit.> that explain how to
implement a comparator only by using a quantum adder. The comparison algorithm works by
computing the high-bit of the expression a - b. If this high-bit is in the state |1⟩
then a < b.

In order to compute the high-bit of a - b, several options are open. The two most promising options are described in the following paragraphs.

The first option is to use a subtractor acting on n+1 qubits and behaving nicely on underflows (i.e. underflows result in cycling to the highest-value), as illustrated in
<ref>. 
This approach requires 2 calls to the subtractor and 1 additional 2-qubit quantum
gate.



Another solution would be to use <ref> to change the subtraction into
an addition and then use a specialised procedure to compute the high-bit of the addition of
two numbers a and b (a being encoded on a quantum register and b a constant).
Computing the high-bit of an addition between a quantum register and a constant can be performed
with the  gate introduced in <cit.>. This approach requires n Toffoli,  and  gates.

















Each of the described methods has its advantages and drawbacks.

For example, the first method crucially relies on a quantum subtractor, and will have the same properties as the subtractor used. In our specific case, we use the subtractor implemented with Drapper's adder <cit.> as explained in <ref>, which in turn uses the quantum Fourier transform. The main disadvantage of using the QFT when looking at practical implementation on quantum hardware is that the QFT involves phase gates with exponentially small angles. These gates may be implemented correctly up to a given threshold, but very small rotation angles will inevitably not be as precise as normal rotation angles due to the hardware limitations in precision. This problem can be circumvented by using an approximate QFT algorithm <cit.> that will cut all the rotation gates that have a rotation angle smaller than a given threshold from the generated circuit but the algorithm will not be exact anymore (small probability of incorrect result).

On the other hand, the  gate involves only ,
controlled- and Toffoli gates. This restriction makes this implementation
more robust than the first one to hardware approximations. Another difference is the
connectivity needed by the approaches: the first method relies on a adder implemented with the
quantum Fourier transform, which use an all-to-all connectity whereas the  gate, once
the qubits correctly ordered, only contains gates on adjacent qubits. As a side note, the exclusive
use of logical gates , controlled- and Toffoli may allow us to
simulate efficiently the  gate on classical hardware as it only involves classical
arithmetic.

As a last word, in the future, the QFT may be implemented directly into the
hardware chips to make it more efficient because it is one of the most used
quantum procedure (and so one of the best candidate for optimisation). Taking
this possibility into account seems a little premature right now but may have a high
impact on the efficiency and precision of the first solution presented.

After summarising all the drawbacks and advantages, we decided to use the
arithmetic comparator for its linear number of gates, because it is based on
arithmetic which does not involve exponentially small rotation angles and because
the need to have n-1 dirty qubits to lend to the procedure is not an issue in
our implementation.



  §.§.§ The  gate


The last gate the oracle implementation will need is an  gate, testing the
equality between an integer stored in a quantum register and a generation-time constant
integer.

This gate has been implemented with a multi-controlled Tofolli gate and a few  gates
before and after the control qubits of the Toffoli gates that should be equal to |0⟩.
The  gates are necessary because a raw Toffoli gate set its target qubit only when all
its control are in the state |1⟩, but we want each control qubit to be equal to a specific
bit of the generation-time constant integer, which can be either |0⟩ or |1⟩.

An implementation example is available in <ref>.



Implementing a  gate controlled by n qubits can be done with only one ancilla qubit
or n-2 garbage qubits and requires n , controlled- or Toffoli gates <cit.>.



 §.§ Oracles for H_1


As noted in <ref>, the oracles M_1 and S_1 can be
optimized by using the fact that they can encode anything for |x⟩_x when the
xth row of H_1 is empty.

We decided to use this optimization opportunity to add regularity to the description
of the H_1 matrix. The implemented matrix H_1, denoted as H_1^impl, is described
in <ref>.


  All indices start at 0. The first row of a matrix has the index 0, the second
  row the index 1 and so on.
  This convention is used to match Python's indexing that starts at 0.









    H_1^impl =
        [                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             [            ^N_c          ^N_c+1 ^2^q-(2N_c + 1);                 ]; [            [origin=c]90N_c{.;        [origin=c]90N_c + 1{.; [origin=c]902^q - (2N_c+1){.;                              ][ 0 ⋯ ⋯ 0 0 1 0 ⋯ 0 0 ⋯ ⋯ ⋯ ⋯ 0; ⋮     ⋮ ⋮ ⋱ ⋱ ⋱ ⋮ ⋮         ⋮; ⋮     ⋮ ⋮   ⋱ ⋱ 0 ⋮         ⋮; 0 ⋯ ⋯ 0 0 ⋯ ⋯ 0 1 0         ⋮; 0 ⋯ ⋯ 0 0 ⋯ ⋯ ⋯ 0 0         ⋮; 1 ⋱   ⋮ ⋮       ⋮ 0         ⋮; 0 ⋱ ⋱ ⋮ ⋮       ⋮ ⋮         ⋮; ⋮ ⋱ ⋱ 0 ⋮       ⋮ ⋮         ⋮; 0 ⋯ 0 1 0 ⋯ ⋯ ⋯ 0 0 ⋯ ⋯ ⋯ ⋯ 0; 0 ⋯ ⋯ 0 0 0 ⋯ ⋯ 0 0 ⋯ ⋯ ⋯ ⋯ 0; ⋮       ⋱ 0 ⋱   ⋮ ⋮         ⋮; ⋮         ⋱   ⋱ ⋮ ⋮         ⋮; ⋮           ⋱ 0 0 0         ⋮; ⋮             ⋱ 0 0 ⋱       ⋮; 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 0 0 0 0 ⋯ ⋯ 0;   ] ]



According to the shape of the matrix in <ref>, the oracle M_1
should implement the transformation

    M_1| x ⟩_x | 0 ⟩_m ↦| x ⟩_x ⊗| x + (N_c + 1) ⟩_m    if  x < (N_c + 1) 
    | x ⟩_x ⊗| x - (N_c + 1) ⟩_m    else
    .


M_1 can be easily implemented with the quantum circuit depicted in 
<ref>.



The oracle V cannot be simplified using the results from
Claim <ref>. It should implement
the transformation written in <ref>.


    V_1| x ⟩_x | 0 ⟩_v ↦| x ⟩_x | 1 ⟩_v    if  (x < 2N_c + 1)  (x ≠ N_c) 
    | x ⟩_x | 0 ⟩_v    else
    .


The implementation of the oracle V_1 is depicted in <ref>. The first part, Set, sets the weight qubit to 1 for all |x⟩_x such that x < 2N_c+1. As this does not correspond to the correct expression of V, the second part Correct is here to set the weight register back to |0⟩_v when x == N_c. 




The last oracle left to implement in order to be able to simulate H_1 is S_1, the
oracle encoding the signs of the non-zero entries of H_1.
The convention used to encode the sign of an entry has been taken from <cit.>
and is: a positive sign is encoded as |0⟩_s, a negative sign is encoded as |1⟩_s. 
As shown in <ref>, H_1 only contains positive non-zero
entries so the sign oracle S_1 should implement the simple transformation of
<ref>: the identity.


    S_1| x ⟩_x | 0 ⟩_s ↦| x ⟩_x | 0 ⟩_s




 §.§ Oracles for H_-1


The matrix H_-1 has less regularity than H_1, which will lead to a more complex
implementation. The implemented matrix H_-1, denoted as H_-1^impl,
is described in <ref>.






    H_-1^impl =
        [                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               [            ^N_c          ^N_c+1 ^2^q-(2N_c + 1);                 ]; [            [origin=c]90N_c{.;        [origin=c]90N_c + 1{.; [origin=c]902^q - (2N_c+1){.;                              ][  0  ⋯  ⋯  0  1  0  ⋯  ⋯  0  0  ⋯  ⋯  ⋯  ⋯  0;  ⋮        ⋮  0 -1  ⋱     ⋮  ⋮              ⋮;  ⋮        ⋮  ⋮  ⋱  ⋱  ⋱  ⋮  ⋮              ⋮;  0  ⋯  ⋯  0  0  ⋯  0 -1  0  ⋮              ⋮;  1  0  ⋯  0  0  ⋯  ⋯  ⋯  0  ⋮              ⋮;  0 -1  ⋱  ⋮  ⋮           ⋮  ⋮              ⋮;  ⋮  ⋱  ⋱  0  ⋮           ⋮  ⋮              ⋮;  ⋮     ⋱ -1  0           ⋮  ⋮              ⋮;  0  ⋯  ⋯  0 -0  ⋱        ⋮  0  ⋯  ⋯  ⋯  ⋯  0;  0  ⋯  ⋯  ⋯  0 -0  ⋱     ⋮  0  ⋯  ⋯  ⋯  ⋯  0;  ⋮              ⋱  ⋱  ⋱  ⋮  ⋮              ⋮;  ⋮                 ⋱  ⋱  0  ⋮              ⋮;  ⋮                    ⋱ -0  0              ⋮;  ⋮                       ⋱ -0  ⋱           ⋮;  0  ⋯  ⋯  ⋯  ⋯  ⋯  ⋯  ⋯  ⋯  0 -0  0  ⋯  ⋯  0;    ] ]



Following the placement of the non-zero and the 0 or -0 entries in the matrix
H_-1^impl of <ref>, the oracle M_-1 should implement
the transformation

    M_-1| x ⟩_x | 0 ⟩_m ↦| x ⟩_x | x + N_c ⟩_m    if  x < N_c 
    | x ⟩_x | x - N_c ⟩_m    else
    .


This transformation is quite similar to the one implemented by the oracle M_1 in
<ref>: N_c + 1 from the transformation M_1 has
been replaced by N_c in the transformation M_-1. Thanks to this similarity,
the implementation of M_-1 will be a nearly-exact copy of the implementation of
M_1. The full implementation of the M_-1 oracle is depicted in 
<ref>.



The weight oracle V_-1 is the simplest to implement for the matrix H_-1, even if
it cannot take advantage of the optimisation discussed in
Claim <ref>.
The transformation that should be implemented by the oracle V_-1 is shown in
<ref>.


    V_-1| x ⟩_x | 0 ⟩_v ↦| x ⟩_x | 1 ⟩_v    if  x < 2N_c 
    | x ⟩_x | 0 ⟩_v    else
    .


The implementation of the weight oracle V_-1 is illustrated in 
<ref>.



The last oracle left to implement is S_-1, the sign oracle. Due to the sign irregularity
in the matrix H_-1^impl, the implementation of S_-1 is more involved and requires several ancillary
qubits.
According to the shape of the matrix H_-1^impl, the sign oracle S_-1 should implement the
transformation defined in <ref>.


    S_-1| x ⟩_x | 0 ⟩_s ↦| x ⟩_x | 0 ⟩_s    if  (x = 0)  (x = N_c)
    | x ⟩_x | 1 ⟩_s    else
    .


An implementation of the oracle S_-1 is illustrated in <ref>. 




§ IMPLEMENTATION OF HIGHER-ORDER LAPLACIANS


The results shown in this paper have all been generated using the second-order discretisation formula shown in <ref>. Higher-order formulas are studied in <cit.>.


  As shown in <cit.>, higher-order discretisations with Neumann boundary conditions are not implementable using the algorithm described in <ref>.


In this appendix we replace the second-order formula given in <ref> and used all along the paper by the fourth-order formula given in <cit.> and re-written below


    ∂^2 ϕ/∂ x^2(i δ x, t) ≈1/δ x^2( 5/2ϕ_i,t - 4/3(ϕ_i-1,t + ϕ_i+1,t) + 1/12( ϕ_i-2,t + ϕ_i+2,t) ).


We are left to devise the matrix B_d^4 that satisfy B_d^4 B_d^4^† = Δ_4 where Δ_4 is the discretisation matrix arising from the fourth-order finite-differences approximation in <ref>.

<cit.> devised an analytic formula for B̂_d^4, the B_d^4 matrix with periodic boundary conditions, using the matrix Ŝ representing the cyclic permutation {1, 2, …, N} with entries Ŝ_i,j = δ_i, (j+1  N) as shown in <ref>.


    Ŝ =
      [ 0 0 ⋯ ⋯ ⋯ 0 1; 1 0         0; 0 ⋱ ⋱       ⋮; ⋮ ⋱ ⋱ ⋱     ⋮; ⋮   ⋱ ⋱ ⋱   ⋮; ⋮     ⋱ ⋱ ⋱ ⋮; 0 ⋯ ⋯ ⋯ 0 1 0;   ]


With this definition of Ŝ, the analytic formula for B̂_d^4 is given in <ref>, with b and c being solution of <cit.>. The exact values for b and c are:

    b = ±1/2√(3)√(7 ± 4 √(3))
    
      c = ±√(7 ± 4√(3)/12)

with the ± signs that can be chosen freely. Note that, b = ±1/2√(3)√(7 ± 4√(3)) and c = 1/12b are irrational because √(3)√(7 ± 4√(3)) = √(3)√(2 + √(3))^2 = √(3)(2 + √(3)) = 2√(3) + 3 is irrational.

<ref> shows the matrix shape with its entries.

    B̂_d^4    = c Ŝ - (b + c) * 𝕀 + b Ŝ^†
       ≈[ b + c     b     0     ⋯     ⋯     ⋯     0     c;     c b + c     b     ⋱                       0;     0     c     ⋱     ⋱     ⋱                 ⋮;     ⋮     ⋱     ⋱     ⋱     ⋱     ⋱           ⋮;     ⋮           ⋱     ⋱     ⋱     ⋱     ⋱     ⋮;     ⋮                 ⋱     ⋱     ⋱     b     0;     0                       ⋱     c b + c     b;     b     0     ⋯     ⋯     ⋯     0     c b + c;       ]


Because periodicity has not been studied in the main use-case of this paper, we would like to also remove the need of periodic boundary conditions in this higher-order laplacian discretisation. This can be achieved by removing the upper-right entry of Ŝ by changing it from 1 to 0. The resulting matrix S is shown in <ref>.


    S =
      [ 0 0 ⋯ ⋯ ⋯ ⋯ 0; 1 0         ⋮; 0 ⋱ ⋱       ⋮; ⋮ ⋱ ⋱ ⋱     ⋮; ⋮   ⋱ ⋱ ⋱   ⋮; ⋮     ⋱ ⋱ ⋱ ⋮; 0 ⋯ ⋯ ⋯ 0 1 0;   ]


Using the exact same formula we can devise B_d^4:

    B_d^4    = c S - (b + c) * 𝕀 + b S^†
       ≈[ b + c     b     0     ⋯     ⋯     ⋯     ⋯     0;     c b + c     b     ⋱                       ⋮;     0     c     ⋱     ⋱     ⋱                 ⋮;     ⋮     ⋱     ⋱     ⋱     ⋱     ⋱           ⋮;     ⋮           ⋱     ⋱     ⋱     ⋱     ⋱     ⋮;     ⋮                 ⋱     ⋱     ⋱     b     0;     ⋮                       ⋱     c b + c     b;     0     ⋯     ⋯     ⋯     ⋯     0     c b + c;       ]


Replacing B_d^4 in <ref> we obtain

    H_d^4 = 1/δ x[     0     ⋯     ⋯     ⋯     0 b + c     b     0     ⋯     0;     ⋮                       ⋮     c b + c     b     ⋱     ⋮;     ⋮                       ⋮     0     c     ⋱     ⋱     0;     ⋮                       ⋮     ⋮     ⋱     ⋱ b + c     b;     0     ⋯     ⋯     ⋯     0     0     ⋯     0     c b + c; b + c     c     0     ⋯     0     0     ⋯     ⋯     ⋯     0;     b b + c     ⋱     ⋱     ⋮     ⋮                       ⋮;     0     ⋱     ⋱     ⋱     0     ⋮                       ⋮;     ⋮     ⋱     ⋱ b + c     c     ⋮                       ⋮;     0     ⋯     0     b b + c     0     ⋯     ⋯     ⋯     0;       ].


One of the main difference with the second-order approximation used all along this paper is that, with the fourth-order approximation, the entries of the matrix H_d^4 are no longer multiples of a common number α∈ℝ. This means that we cannot write down H_d^4 as an integer weighted matrix multiplied by a real number, and so the trick used in <ref> to simulate the integer weighted matrix H_d for a time α t is no longer applicable.

Consequently, and independently of the decomposition we use for Ĥ_̂d̂^̂4̂, at least one of the matrices in the decomposition of Ĥ_̂d̂^̂4̂ will not be “easy to simulate” as defined in Definition <ref>.

Ultimately, the main consequence of this observation is that we will have to use a real-weighted hamiltonian simulation procedure. Such a procedure can be found in <cit.> but requires to approximate the real-weighted entries with a fixed-point representation that has at least 2 evident caveats:



  * It is impossible to encode the irrational numbers b and c exactly with a fixed-point representation. This means that we add another layer of approximation, even before the approximation caused by the use of a product-formula.

  * The hamiltonian simulation procedures used for real numbers requires more qubits. More precisely, the number of additional qubits required depends on the desired precision ϵ and grows as log_2(1/ϵ). 



  Even if the H_d^4 matrix seems quite hard to simulate, it is still a 3-sparse matrix. This means that it is still managable to hand-write the oracles. Moreover, having a small number of matrices in the decomposition helps in reducing the error introduced by product-formulas.





§ OPTIMISATION OF THE IMPLEMENTATION


Once the correctness of the implementation validated, one of the most important remaining work is to try to optimise the implementation. The optimisation of a software is often performed as an iterative task.

The first step is to define a figure of merit, a quantity we want to minimise during the optimisation process. Among the most obvious figures of merit are the total number of gates, the number of CNOT gates or the total execution time of the quantum program. More complex quantities can also be considered, such as the execution time using error correction codes or the final state fidelity. In this paper we decided to take into account an estimation of the total execution time of the quantum program on an imaginary device that shares today's chips characteristics.

The second step of the optimisation process consists in isolating the subroutines that contribute the most to the figure of merit. As an example, if the quantum program spend 90% of the total execution time in one subroutine, this subroutine should be the first place to look for optimisations.

After the isolation of one or two subroutines, the actual optimisation can take place. The goal of this third step is to decrease the impact of the subroutines considered on the overall figure of merit without changing the final result of the implementation.

Finally, once the optimisation is performed, the optimisation process can be repeated by re-starting at the second step, until the program is considered sufficiently optimised.

One of the main difficulty we encountered when applying this optimisation process was to correctly isolate the most time-consumming subroutines. In classical computing, this step is usually performed with specialised tools such as  or a more advanced profiler, but no such tool exist for quantum programs. In order to fill this gap we developped , a tool that analyses a quantum program and generates a report similar to the one generated by . Using  and some of the various tools compatible with , we plotted the call-graph shown in <ref>.

From this call-graph, it is clear that the adder is the most costly subroutine and that it should be optimised. The adder internally uses the Quantum Fourier Transform (QFT), which takes more than 50% of the total execution time. The issue is that the QFT implementation is already very concise and we do not expect to be able to optimise it enough to cut significantly its overall cost. This leads us to the conclusion that a new algorithm that do not require the QFT should be used to implement an adder. Such an algorithm can be found in <cit.>.

Changing the implementation of the adder from Draper's adder to the arithmetic-based adder from <cit.> improves drastically the total execution time of the quantum program and produce the  call-graph in <ref>.





§ GATE COUNT ANALYSIS





 §.§ Precise subroutines gate counts



<ref> summarise the gate count and ancilla qubit requirements for all the major subroutines used in the wave equation solver implementation. Using the entries of this table, it is
possible to compute an estimation of the number of gates required to solve the wave equation. 

As explained in <ref>, we need to simulate each of the 1-sparse Hamiltonians in the decomposition. Aggregating the estimates in <ref> we obtain the costs in <ref> for the Hamiltonian simulation part. Note that the cost of the adder has been voluntarily omited from the computations in order to be able to compare the cost with different adder implementations. Let a be the gate cost of the adder implementation chosen, the cost of simulating the Hamiltonian needed to solve the 1-dimensional wave equation is: 82n - 35 + 12a_𝚃𝚘𝚏𝚏𝚘𝚕𝚒 Toffoli gates, 84n + 21 + 9[ 0, n-1 ] + 12a_𝙲𝙽𝙾𝚃  gates and n + 12a_1-𝚚𝚞𝚋𝚒𝚝 1-qubit gates.

Choosing an adder implementation and simplifying the gate counts by omitting negligible terms we obtain the gate counts summarised in <ref>. It is interesting to note that even if the arithmetic-based adder adds huge constants in the gate count, it does not change the asymptotic complexity whereas Draper's adder changes the number of  gates required from n to n^2.










 §.§ Impact of the precision requirements



The gate counts presented in <ref>, <ref> and <ref> are only valid when the precision of the solver is not accounted for. When the solver precision matters, an additional step that consists is splitting the Hamiltonian Simulation into r steps needs to be performed as noted in <cit.>.

Several bounds exist to determine a r∈ℕ^* that will analytically ensure that the
maximum allowable error ϵ is not exceeded. The definition of such bounds can be found
in <cit.> and <cit.>.

The first bound has been devised by analytically bounding the error of simulation due to the
Trotter-Suzuki formula approximation by ϵ_0

    ||exp( -it ∑_j=0^m-1 H_j ) - [ S_2k( -it/r) ]^r||⩽ϵ_0

and then let ϵ_0 ⩽ϵ for a given desired precision ϵ.
If we let Λ = max_j || H_j|| and

    τ = 2m5^k-1Λ| t |

then 

    r^ana_2k = ⌈max{τ, √(eτ^2k+1/3ϵ)}⌉.

This bound is called the analytic bound.

A better bound called the minimised bound can be devised by searching
for the smallest possible r that satisfies the conditions detailed in
<cit.>. This bound is rewritten in Equation
(<ref>).


    r^min_2k = min{ r ∈ℕ^* : τ^2k+1/3r^2kexp( τ/r) < ϵ}


Another bound involving nested commutators of the H_i is described in <cit.> and gives

    r_2k^comm∈α_comm^1/2k t^1 + 1/2k/ϵ^1/2k

where k is the order of the product-formula used, t the time of simulation, ϵ the error and

    α_comm = ∑_i_0, i_1, …, i_p = 0^m-1||[ H_i_p, …[ H_i_1, H_i_0] ] ||.


Once the value of r has been computed, the quantum circuit simulating the matrix H for a time t/r should be repeated r times. This adds a factor of r in front of all the gate counts computed in <ref>, <ref> and <ref>.



 §.§ Impact of error-correction



When error-correction is studied, two gates are particularly important: T and Toffoli gates. The T gate has a prohibitive cost when compared to the Clifford quantum gates and implementing a Toffoli gate requires 7 of such T gates as noted in <cit.> and <cit.>.

<ref> summarise the cost of the non Clifford quantum gates used in the implementation of the 1-dimensional wave equation solver. The rotation gates need to be approximated. One solution to approximate the R_n and P_h gates is given in <cit.>. In order to obtain practical results as opposed to theoretical ones, we chose to use the number computed in <cit.>.

The final T-count is summarised in <ref>. From <ref> it is clear that the  implementation is more efficient than the  one. 





