# Quantum algorithms to simulate quadratic classical Hamiltonians and optimal
control

Hari Krovi

###### Abstract

Simulation of realistic classical mechanical systems is of great importance to
many areas of engineering such as robotics, dynamics of rotating machinery and
control theory. In this work, we develop quantum algorithms to estimate
quantities of interest such as the kinetic energy in a given classical
mechanical system in the presence of friction or damping as well as forcing or
source terms, which makes the algorithm of practical interest. We show that
for such systems, the quantum algorithm scales polynomially with the logarithm
of the dimension of the system. We cast this problem in terms of Hamilton’s
equations of motion (equivalent to the first variation of the Lagrangian) and
solve them using quantum algorithms for differential equations. We then
consider the hardness of estimating the kinetic energy of a damped coupled
oscillator system. We show that estimating the kinetic energy at a given time
of this system to within additive precision is ${\mathsf{BQP}}$ hard when the
strength of the damping term is bounded by an inverse polynomial in the number
of qubits. In other words, when combined with prior work on coupled
oscillators, this shows that we can get a quantum advantage even in the
presence of damping.

We then consider the problem of designing optimal control of classical
systems, which can be cast as the second variation of the Lagrangian. In this
direction, we first consider the Riccati equation, which is a nonlinear
differential equation ubiquitous in control theory. We give an efficient
quantum algorithm to solve the Riccati differential equation well into the
nonlinear regime. To our knowledge, this is the first example of any nonlinear
differential equation that can be solved when the strength of the nonlinearity
is asymptotically greater than the amount of dissipation. We then show how to
use this algorithm to solve the linear quadratic regulator problem, which is
an example of the Hamilton-Jacobi-Bellman equation. Our results show that
quantum algorithms for differential equations, especially for ordinary linear
differential equations, can go quite far in solving problems of practical
relevance.

##  1 Introduction

Quantum algorithms have been designed for several important problems such as
cryptography, quantum Hamiltonian simulation, condensed matter physics, semi-
definite programming and optimization. Important primitives have been
developed such as quantum phase estimation, quantum walks, Fourier transforms,
Hamiltonian simulation and quantum linear systems algorithms (QLSA). The
latter two have a close relationship and common techniques can be used in
implementing them. In fact, it was shown in [1, 2] that one can include
quantum walks, phase estimation and Hamiltonian simulation to show that they
can all be cast as the problem of implementing functions of Hamiltonians. This
has recently been generalized to a large set of non-Hermitian matrices in [3].
Hamiltonian simulation essentially solves the Schrödinger equation, which is a
linear ordinary differential equation (ODE). Efficient techniques for this
have been developed in a series of papers [4, 5, 6, 7, 8] with current
techniques having an optimal dependence on the parameters. Hamiltonian
simulation algorithms for time-dependent Hamiltonians have been developed in
[9, 10, 11].

The HHL algorithm to solve linear systems of equations was described in the
seminal paper [12]. This algorithm has been generalized in subsequent works
(see [13] for a review), where the matrix representing the linear system can
be non-Hermitian, rectangular and even singular. The dependence on various
parameters such as the condition number and error has been improved in [14,
15, 16]. This algorithm has been applied to problems in topological data
analysis [17, 18], machine learning [19] and more recently in solving
differential equations (described in detail below). An important subroutine
developed in the process is a protocol to load classical data called block-
encoding. This method takes classical data such as matrix entries of sparse
and structured matrices and creates a unitary that represents the matrix in
the quantum algorithm. Block-encoding has been used in earlier quantum
algorithms such as [15] and recently many explicit methods have been developed
[1, 8, 20, 21, 22]. We use these techniques in this paper as well.

An approach to solve differential equations using the linear systems
algorithms was given in [23]. This was made into a high-precision algorithm in
[24] using truncated Taylor series to approximate $\exp(At)$. These algorithms
have been generalized to allow for arbitrary time-independent matrices in
[25]. An important quantity used in this work is the norm of the exponential
of the matrix i.e., $\|\exp(At)\|$. It was shown that as long as this norm of
the exponential of the matrix $A$ over the time domain of interest (denoted
$C(A)$ in [25]) is polynomially bounded, one can bound the condition number of
the resulting linear system. Using a stronger condition that the log-norm of
$A$ is negative, this has been extended to time-dependent matrices in [26,
27]. More recently, an approach to reduce time-dependent differential
equations to Hamiltonian simulation is given in [28, 29] where optimal
dependence on parameters is achieved assuming $A$ has a negative log-norm. It
is interesting to note that for one of the applications considered here (of
damped coupled oscillators), the matrix involved in the differential equation
does not have a negative log-norm and needs the weaker condition of bounded
$C(A)$. A general theory of quantum algorithms for linear ODEs is given in
[30] where a fast-forwarding theorem similar to the one in Hamiltonian
simulation is proved. Recently, there has been work on using Hamiltonian
simulation algorithms to solve general ODEs. In [31], the authors reduce a
general class of partial differential equations to Hamiltonian simulation.

Many of the above algorithms essentially solve a _quantum_ problem i.e., given
an initial state as a quantum state, and sparse access oracles, the goal is to
create a quantum state that encodes the normalized solution state. However, in
practical applications, one usually has access to classical data that must be
encoded into the quantum state and at the end of the simulation we would also
need to extract classical data out of the quantum state. Many quantum
algorithms are designed within a specific mathematical framework and applying
them to problems of practical interest is not always straightforward. The
comprehensive review [32] covers the above mentioned algorithms and many more
from the point of view of end-to-end applications. In considering such end-to-
end applications, one needs proof or strong evidence of three different
aspects namely, quantum efficiency, classical hardness and practical utility.

The recent work on estimating the kinetic energy of coupled oscillators [33]
is an example of such an end-to-end application. In that paper, the authors
consider a problem of practical interest and show that quantum algorithms for
Hamiltonian simulation can be used to provide evidence of a quantum advantage.
The specific problem considered is that of estimating kinetic energy of a
network of coupled oscillators. The problem of estimating the kinetic energy
(or more precisely, the decision version where one has to decide whether the
kinetic energy is above some value or below another with the promise that the
gap is a constant independent of the number of qubits) has been shown to be
${\mathsf{BQP}}$ complete and therefore, one does not expect an efficient
classical algorithm unless ${\mathsf{BQP}}={\mathsf{BPP}}$. While this problem
has a quantum advantage, an important caveat is that the oscillators are
assumed to be ideal and there is no damping in the system. Going beyond this
idealized setting i.e., simulating oscillators in the presence of damping
would increase the practical utility of the quantum algorithm. In this work,
we show that even in the presence of damping in coupled oscillators (or
dissipation in other systems that can be modeled by the similar equations),
there exist efficient quantum algorithms to approximate the kinetic energy of
the system.

One might wonder if the presence of damping makes the problem of approximating
the kinetic energy of an oscillator network amenable to classical algorithms
i.e., it may be that the problem of kinetic energy estimation of damped
coupled oscillators is in ${\mathsf{BPP}}$. We show that this is unlikely to
be the case. Specifically, we show that even in the presence of damping,
approximating the kinetic energy of an arbitrary sparse oscillator network is
${\mathsf{BQP}}$ hard (provided that the strength of damping is inverse
polynomial in the number of qubits). This means that even under the realistic
conditions, this problem has a quantum advantage. In order to construct a
quantum algorithm for this problem, we cast this problem as ordinary
differential equation (ODE) and use the recently developed quantum algorithm
for ODEs [25].

In order to apply quantum algorithms for linear ODEs of the form
$\dot{x}=Ax+b$, one needs to overcome certain obstacles. As pointed out in
[33], if the quantum algorithm for a linear ODE requires that $A$ has a
negative log-norm, then one cannot get an efficient algorithm. The algorithm
in [25] requires that the norm of the exponential of $A$ is bounded and works
even when the log-norm of $A$ is positive. However, this encoding means that
the norm of the exponential of $A$ is proportional to the condition number of
the matrix $V$, which encodes the spring constants. The algorithm of [33] did
not depend on the condition number of $V$, but rather only the norm of $V$.
This is the price one has to pay to be able to simulate damping (which is not
amenable to Hamiltonian simulation). Fortunately, the work of [33] also showed
the the hardest cases of the coupled oscillator problem (i.e., the instances
which are ${\mathsf{BQP}}$ hard) have a constant condition number for $V$.
This also carries over for the case with damping. We also show that in a
different basis, we can design a quantum algorithm to estimate the kinetic
energy with a slightly worse dependence on error, but no dependence on the
condition number of $\sqrt{V}$. Therefore, depending on the condition number
of $\sqrt{V}$, one can choose one of these bases. These algorithms are
applicable to systems other than just damped forced coupled oscillators. For
the specific case of damped and forced coupled oscillators, we show that one
can remove the dependence on the condition number of $V$ and the worse
dependence on error by using the specific structure of the problem (by
leveraging the results for block-encoding in [33]).

Quantum algorithms for nonlinear differential equations have been constructed
in [34, 25, 35, 36] and quantum algorithms for partial differential equations
(PDEs) in [37, 38, 39]. For nonlinear ODEs, quantum algorithms such as [34,
25] can be applied to very general types of nonlinear ODEs (with a polynomial
nonlinearity). However, these algorithms can only handle a small amount of
nonlinearity i.e., they require $R<1$, where $R$ is (roughly) the ratio of the
strength of the nonlinearity to the strength of dissipation. In [34], it was
also shown that no quantum algorithm can efficiently solve nonlinear ODEs for
$R>\sqrt{2}$. This was improved to $R\geq 1$ in [40] matching the algorithm.
While this is true for arbitrary nonlinear differential equations, one might
hope that for specific nonlinear ODEs, higher nonlinearity may be efficiently
solvable using quantum algorithms. However, there are currently no quantum
algorithms that can solve nonlinear ODEs for $R>1$. In this paper, we remedy
this situation and show that a specific nonlinear differential equation called
the Riccati equation used in optimal control can be solved for
$R=\operatorname{polylog}(N)$, where $N$ is the dimension of the matrices
involved in the equation. The solution to the Riccati equation, which is
ubiquitous in classical optimal control, has practical applications such as
the stability of magnetohydrodynamics [41].

The solution to the Riccati equation can also be used to solve another problem
in optimal control known as the linear quadratic regulator problem (LQR). This
problem is about finding a controller that can regulate a process such as in
an airplane by minimizing a cost function. It turns out that this can be cast
as a Riccati equation by viewing it as a special case of the Hamilton-Jacobi-
Bellman (HJB) equation. The HJB equation is a non-linear partial differential
equation that generalizes the Hamilton-Jacobi equation from the theory of
calculus of variations. It gives the necessary and sufficient conditions for
optimality of a control problem. In its generality, this problem can be highly
nonlinear. However, in the case of the LQR problem, it becomes tractable using
quantum algorithms.

The differential equations used in classical mechanical systems described
above and those in optimal control are linked via a theory known as Calculus
of Variations [42, 43] as mentioned above. This theory has two equivalent
approaches - a Lagrangian formalism and a Hamiltonian formalism. We use the
Hamiltonian formulation of this theory since it gives us linear ODEs. We can
then apply quantum algorithms for linear ODEs to solve the above problems of
simulating classical mechanical systems under realistic conditions and of
solving the Riccati and Hamilton-Jacobi-Bellman equation (for the LQG).
Interestingly, in the above applications, we need the linear ODE solvers that
are based on linear systems rather then Hamiltonian simulation. This is
because we need to extract classical information from the unnormalized final
state (which can be done via the so-called history state that keeps track of
the solution at various times). Recent approaches via Hamitonian simulation
directly produce the normalized final state and it would be interesting to see
if they can be modified for these applications.

The rest of the paper is organized as follows. In Section 2, we give a basic
background in calculus of variations, the Hamiltonian formalism, optimal
control, Riccati equation and the HJB equation. Then in Section 3, we present
some background and definitions needed on oracle access and block encoding as
well as results on quantum linear systems algorithms that we will need later
on. In Section 4, we give our quantum algorithm to estimate the kinetic energy
of damped coupled oscillators and other realistic mechanical systems. We also
give a proof of the ${\mathsf{BQP}}$ hardness of additive approximation of the
kinetic energy. Then in Section 5, we give quantum algorithms to solve the
vector and matrix versions of the Riccati equation and its application to the
LQR problem. Finally, in Section 6, we present our conclusions and open
questions.

##  2 Preliminaries

###  2.1 Calculus of variations

The starting point of the theory of calculus of variations is the Lagrangian
which is a functional that depends on the position $q$, the velocity $\dot{q}$
and time $t$. The quantities $q$ and $\dot{q}$ may also depend on time. The
goal is to find the minimum of the following integral of the Lagrangian.

$J(q)=\int_{t_{i}}^{t_{f}}L(q,\dot{q},t)dt\,,$ (2.1)

subject to some initial and final conditions on $q(t_{i})$ and $q(t_{f})$. In
classical mechanics, this is called the _principle of least action_. The
choice of $q$ and $\dot{q}$ is not unique. Using canonical transformations,
one can move to other coordinate systems such as polar coordinates. In this
paper, we use the above “Cartesian” coordinates so that the solution
corresponds to physical positions and velocities.

For the above minimization problem, one can derive a necessary condition
similar to the condition that the vanishing at a point of the first derivative
of a function gives its extremum. For the functional above, this is equivalent
to the vanishing of the first variation, which leads to the famous Euler-
Lagrange equations of motion. These equations can be written as

$\frac{\partial L}{\partial q_{i}}=\frac{d}{dt}\frac{\partial
L}{\partial\dot{q}_{i}}\,,$ (2.2)

for all coordinate indices $i$. For example, when the mechanical system is a
set of masses connected by springs, the Lagrangian describes a system of
coupled oscillators. The equations of motion for this system were recently
explored in [33], where quantum algorithms for simulation and the complexity
of estimating energies were given. The equations of motion in this case become

$M\ddot{q}(t)=Fq(t)\,,$ (2.3)

where $M$ is a diagonal matrix of masses and $F$ is a matrix that captures the
force due to the spring constants.

There is a different set of equations, called Hamilton’s canonical equations,
that are first order unlike the Euler-Lagrange equations (which are second
order equations). Hamilton’s equations are equivalent to the Euler-Lagrange
equations of motion and are derived from the classical Hamiltonian. To define
the classical Hamiltonian, we first need to define the momentum $p$ which is
conjugate to the position coordinate vector $q$ as follows.

$p_{i}=\frac{dL}{d\dot{q}_{i}}\,,$ (2.4)

where $p_{i}$ are the components of $p$. Then one takes the Legendre transform
to obtain the Hamiltonian as follows.

$\mathcal{H}(q,p,t)=p\dot{q}-L(q,\dot{q},t)\,,$ (2.5)

where $p\dot{q}$ is the dot product of $p$ and $\dot{q}$. The set of $2N$
coordinates $q_{i}$ and $p_{i}$ constitute the classical phase space of the
mechanical system. For instance, the classical Hamiltonian of a system of
coupled oscillators can be written as

$\mathcal{H}=\frac{1}{2}p^{T}M^{-1}p+q^{T}Kq\,,$ (2.6)

where $K$ is the potential term of the spring constants. The equations of
motion coming from the classical Hamiltonian, also known as Hamilton’s
canonical equations, can be written as

$\dot{q}=\nabla_{p}\mathcal{H}(q,p,t)\,,\hskip
7.22743pt\dot{p}=-\nabla_{q}\mathcal{H}(q,p,t)\,,$ (2.7)

where

$\nabla_{q}=\begin{pmatrix}\partial_{q_{1}}\\\ \vdots\\\
\partial_{q_{N}}\end{pmatrix}\,,\hskip
7.22743pt\nabla_{p}=\begin{pmatrix}\partial_{p_{1}}\\\ \vdots\\\
\partial_{p_{N}}\end{pmatrix}\,.$ (2.8)

This set of equations is equivalent to the Euler-Lagrange equations but are
first order rather than second order equations. In the following, we present
some known results [44] on quadratic Hamiltonians that are relevant to this
paper.

###  2.2 Quadratic Hamiltonians

Suppose we have a quadratic Hamiltonian of the form

$\mathcal{H}(q,p,t)=\frac{1}{2}(q^{T}Q_{1}q+q^{T}Q_{2}p+p^{T}Q_{3}q+p^{T}Q_{4}p)\,.$
(2.9)

We can write this compactly as

$\mathcal{H}(z,t)=\frac{1}{2}z^{T}Qz\,,$ (2.10)

where

$Q=\begin{pmatrix}Q_{1}&Q_{2}\\\ Q_{3}&Q_{4}\end{pmatrix}\,,\hskip
7.22743ptz=\begin{pmatrix}q\\\ p\end{pmatrix}\,.$ (2.11)

We can then exploit the rich symplectic structure of the phase space of such
systems. Below, we give some basic information about this structure that we
will use in our algorithms later.

For any integer $n$, a symplectic matrix $S$ is a $2n\times 2n$ matrix (over
the reals) such that

$S^{T}JS=J\,,$ (2.12)

where $S^{T}$ is the transpose of $S$ and $J$ is defined as the following
block matrix

$J=\begin{pmatrix}0&I\\\ -I&0\end{pmatrix}\,.$ (2.13)

The set of all symplectic matrices forms a group under multiplication denoted
Sp$(2n,\mathbb{R})$. This continuous group also has the structure of a Lie
group whose Lie algebra is given by

$\mathfrak{sp}(2n,\mathbb{R})=\\{X\in M_{2n\times 2n}:JX+X^{T}J=0\\}\,.$
(2.14)

In terms of block matrices, the symplectic algebra can be written as

$X=\begin{pmatrix}A&B\\\ C&D\end{pmatrix}\,,$ (2.15)

such that

$\displaystyle A=-D^{T}$ (2.16) $\displaystyle C=C^{T}$ (2.17) $\displaystyle
B=B^{T}\,.$ (2.18)

The exponential map takes the lie algebra to the Lie group of symplectic
matrices i.e., $\exp(X)=S$, where $X\in\mathfrak{sp}(2n,{\mathbb{R}})$ and
$S\in$Sp$(2n,{\mathbb{R}})$. Without any additional constraints on $X$, the
exponential map is not a diffeomorphism, which means that one cannot
exponentiate to get to an arbitrary symplectic matrix. In terms of these
quantities, we can write Hamilton’s equation (2.7) as follows

$\frac{dz}{dt}=JQz\,,$ (2.19)

where $Q$ defines the quadratic Hamiltonian in (2.10). The solution to this

$z(t)=\exp(JQt)z(0)\,,$ (2.20)

where $z(0)$ is the initial condition at $t=0$.

When $Q$ is a positive definite matrix, Williamson’s theorem [44] implies that
$JQ$ is diagonalizable (not necessarily unitarily diagonalizable) and its
eigenvalues are all purely imaginary. Loosely speaking, $J$ plays the role of
the complex unit $i=\sqrt{-1}$. When $Q$ is a positive definite matrix, we can
convert the evolution in (2.20) to a quantum Hamiltonian evolution as follows.
First, define $H=i\sqrt{Q}J\sqrt{Q}$. It is easy to see that when $Q$ is
positive, $H$ is Hermitian (using the fact that $J^{\dagger}=-J$). Now define

$y(t)=\sqrt{Q}z(t)\,.$ (2.21)

Then (2.20) becomes

$y(t)=\sqrt{Q}\exp(JQt)(\sqrt{Q})^{-1}y(0)=\exp(\sqrt{Q}J\sqrt{Q}t)y(0)=\exp(iHt)y(0)\,.$
(2.22)

This shows that $y(t)$ evolves according to the quantum Hamiltonian $H$.
However, $y(t)$ may not always be easy to deal with from the point of view of
extracting useful information in a quantum algorithm. At times, we might be
interested in obtaining the solution

$x=\begin{pmatrix}q\\\ \dot{q}\end{pmatrix}\,.$ (2.23)

It can be seen that

$z=\begin{pmatrix}I&0\\\ 0&M\end{pmatrix}x\,,$ (2.24)

where recall that

$z=\begin{pmatrix}q\\\ p\end{pmatrix}\,.$ (2.25)

This means that from (2.19)

$\begin{pmatrix}I&0\\\ 0&M\end{pmatrix}\frac{dx}{dt}=JQ\begin{pmatrix}I&0\\\
0&M\end{pmatrix}x\,,$ (2.26)

This gives us

$\frac{dx}{dt}=\begin{pmatrix}I&0\\\
0&M^{-1}\end{pmatrix}JQ\begin{pmatrix}I&0\\\ 0&M\end{pmatrix}x\,.$ (2.27)

When $Q$ has a special block-diagonal form below (such as in the case of
coupled oscillators)

$Q=\begin{pmatrix}V&0\\\ 0&M^{-1}\end{pmatrix}\,,$ (2.28)

we have

$\frac{dx}{dt}=\begin{pmatrix}0&I\\\ -M^{-1}V&0\end{pmatrix}x\,.$ (2.29)

###  2.3 Quadratic Hamiltonians with dissipation

Non-conservative systems more accurately model real-world mechanical systems.
In this subsection, we consider dissipation, which could arise due to drag,
damping or other resistive forces on conservative systems. In such systems,
one defines the Lagrangian $L$ in the usual way without including dissipation
and then adds terms such as the Rayleigh dissipation function [45]. However,
the Rayleigh dissipation function does not account for forcing or source
terms, which are added separately. All of this modifies the first variation as
follows.

$\frac{d}{dt}\frac{\partial L}{\partial\dot{q}_{i}}-\frac{\partial L}{\partial
q^{i}}=r_{i}\,,$ (2.30)

where $r_{i}$ is models the role of dissipation and sources. The equation of
motion can be written as

$M\ddot{q}+R\dot{q}+Vq+f=0\,.$ (2.31)

In the above equation, $M$ comes from the kinetic energy term and corresponds
the mass matrix, $R$ is a matrix capturing damping or dissipation, $V$ is the
potential energy matrix and $f$ corresponds to sources or forcing terms.
Several interesting systems such as damped coupled oscillators [45], dynamics
of rotating machinery [46] and gyroscopic systems [47] have equations of
motion of this type. To convert this system of second order differential
equations to Hamilton’s first order equations, consider a Lagrangian as
follows.

$L=\frac{1}{2}\dot{q}^{T}M\dot{q}-\frac{1}{2}q^{T}Vq\,.$ (2.32)

We can define the Hamiltonian in the same way as before (via a Legendre
transformation) as

$H=p\cdot\dot{q}-L\,,$ (2.33)

where

$p_{i}=\frac{\partial L}{\partial\dot{q}_{i}}=M_{i}\dot{q}_{i}\,.$ (2.34)

Now one can define the following modification of Hamilton’s canonical
equations [45].

$\displaystyle\dot{q}_{j}=\frac{\partial H}{\partial p_{j}}\,,$ (2.35)
$\displaystyle\dot{p}_{j}=-\frac{\partial H}{\partial q_{j}}+r_{j}\,,$ (2.36)

where $r_{j}$ captures the dissipation and source terms and comes from the
matrix $R$. In our case, these are the terms involving $R$ and $f$. This gives
us

$r_{j}=-\sum_{i}R_{ji}\dot{q}_{i}-f_{j}\,.$ (2.37)

This becomes

$\frac{dz}{dt}=\begin{pmatrix}0&M^{-1}\\\
-V&-RM^{-1}\end{pmatrix}z-\tilde{f}\,,$ (2.38)

where

$z=\begin{pmatrix}q\\\ p\end{pmatrix}\,,\hskip 7.22743pt\text{and}\hskip
7.22743pt\tilde{f}=\begin{pmatrix}0\\\ f\end{pmatrix}\,.$ (2.39)

Next, we consider $y=\sqrt{Q}z$. In this basis, we get

$\frac{dy}{dt}=\sqrt{Q}\begin{pmatrix}0&M^{-1}\\\
-V&-RM^{-1}\end{pmatrix}\sqrt{Q^{-1}}y-\sqrt{Q}\tilde{f}\,,$ (2.40)

Let $b=-\sqrt{Q}\tilde{f}$. Using

$Q=\begin{pmatrix}V&0\\\ 0&M^{-1}\end{pmatrix}\,,$ (2.41)

we get

$\frac{dy}{dt}=\begin{pmatrix}0&\sqrt{VM^{-1}}\\\
-\sqrt{M^{-1}V}&-\sqrt{M^{-1}}R\sqrt{M^{-1}}\end{pmatrix}y+b\,.$ (2.42)

Finally, switching to $x=(q,\dot{q})^{T}$, we get

$\frac{dx}{dt}=\begin{pmatrix}0&I\\\ -M^{-1}V&-M^{-1}R\end{pmatrix}x+b\,,$
(2.43)

where (it turns out due to the form of $\tilde{f}$)

$b=-\begin{pmatrix}I&0\\\ 0&M^{-1}\end{pmatrix}\tilde{f}\,.$ (2.44)

This can be written as

$\frac{dx}{dt}=Ax+b\,,$ (2.45)

where

$A=\begin{pmatrix}0&I\\\ -M^{-1}V&-M^{-1}R\end{pmatrix}\,.$ (2.46)

and $b$ is the transformed source term.

###  2.4 Optimal control

In this subsection, we continue with calculus of variations and discuss
applications to optimal control. In this setting, a _control_ parameter is
included in the Lagrangian which controls the position variable through a
differential equation. The goal is to find the value of the control parameter
as a function of time such that a cost functional is minimized. Such problems
have many practical applications in operations research [43].

In mathematical terms, the problem of optimal control can be described as
follows. Let $q(t)$ be the position variable as before and $u(t)$ be the
control parameter. Suppose that the two are related by an ODE as follows.

$\dot{q}=F(q,u,t)\,,$ (2.47)

where $F$ is some function (which may be nonlinear). Later, we will relax that
condition and in fact, assume a linear relationship. For the general problem
definition however, we will assume $F$ to be arbitrary. The problem is to find
the optimal $u(t)$ such that the total cost $J(u)$ is minimized, where $J(u)$
is defined as

$J(u)=\int_{t_{i}}^{t_{f}}L(q,\dot{q},t)dt+K(t_{f},q_{f})\,,$ (2.48)

Here $L$ is the Lagrangian (also the running cost) and $K(t_{f},q_{f})$ is the
terminal cost function that depends on the terminal time $t_{f}$ and the
terminal position $q_{f}$. Here the final state $q_{f}$ and the final time
$t_{f}$ are either fixed or free. This cost functional is sometimes called the
Bolza functional [43]. One can also convert this to a problem of maximization
by taking the negative of the Lagrangian.

While it is not easy to solve the above problem for arbitrary functions and
cost functionals, there are useful necessary and sufficient conditions that
can be used to check if a given curve is optimal. Often, these conditions
narrow down the search space making it possible to find a local minimum. A
local minimum is defined as the set of all $u^{\ast}$ such that
$J(u^{\ast})\leq J(u)$ for all $|u^{\ast}-u|\leq\epsilon$.

There are three main necessary conditions for any curve $u(t)$ to be a
minimum. The first of these conditions is that the first variation of the
Lagrangian must equal zero leading to the Euler-Lagrange equations described
above. The second is called the Legendre condition and it states that the
second derivative of the Lagrangian with respect to $\dot{q}$ is positive.
More precisely, define the matrix $L_{\dot{q}\dot{q}}$ as

$[L_{\dot{q}\dot{q}}]_{i,j}=\frac{\partial^{2}\mathcal{H}}{\partial\dot{q}_{i}\partial\dot{q}_{j}}\,,$
(2.49)

then the Legendre condition states that

$L_{\dot{q}\dot{q}}\geq 0\,.$ (2.50)

Finally, the third necessary condition is called Weierstrass condition and it
states that the Lagrangian must be locally convex. To state it in mathematical
form, define the so called Weierstrass excess function as follows.

$E(q,v,w,t)=L(q,v,t)-L(q,w,t)-(w-v)L_{v}(q,v,t)\,,$ (2.51)

where $v=\dot{q}$. Then the necessary condition states that $E(q,v,w,t)\geq 0$
for all $t\in[t_{i},t_{f}]$ and for all $w$.

Sufficient conditions are similar to the above necessary conditions but
stronger. The following set of conditions is sufficient for a curve $q(t)$ to
be the (strong local) minimum.

  1. 1.

$q(t)$ and $\dot{q}(t)$ satisfy the Euler-Lagrange equations (or equivalently,
$q$ and $p$ satisfy Hamilton’s canonical equations),

  2. 2.

$L_{\dot{q}\dot{q}}>0$,

  3. 3.

$E(x,v,w,t)\geq 0$, for all $w$ and where $|x-q|\leq\epsilon$ and
$|v-\dot{q}|\leq\epsilon$.

  4. 4.

There are no conjugate points in the interval $[t_{i},t_{f}]$.

The first condition is just the first variation and the second comes from
taking the second variation of the Lagrangian. As pointed out above, these two
conditions are similar to the vanishing of the first derivative and the
positivity of the second derivative as sufficient conditions for the existence
of a minimum of a smooth function. The third condition is same as before. Here
the last condition is sometimes called the Jacobi condition, which we explain
next since this is relevant to our algorithm for the Riccati equation.

Conjugate points are points on a manifold which need not minimize geodesics.
The Jacobi equation can be used to determine if there exist conjugate points
in any interval. Applied to the above Lagrangian, the definition states that
$t_{0}$ and $t_{1}$ are conjugate if there exists a function $h(t)$ satisfying

$\frac{d}{dt}(L_{\dot{q}\dot{q}}\dot{h}+L_{\dot{q}q}h)-L_{q\dot{q}}\dot{h}-L_{qq}h=0\,,$
(2.52)

such that $h(t_{0})=0$ and $h(t_{1})=0$. It turns out that the non-existence
of conjugate points is related to the existence of a solution of the Riccati
equation (which we explain in Section 2.7).

The two main approaches to solving problems in optimal control are the
Pontryagin principle approach (also called the maximum principle) and the
Hamilton-Jacobi-Bellman approach. The approach using the Pontryagin principle
uses Hamilton’s canonical equations along with inequality conditions, where as
the Hamilton-Jacobi-Bellman equation is a partial differential equation. These
are described in detail in [43] for instance. Here, we focus on a special case
called the linear quadratic regulator (LQR) problem, where the cost function
or Lagrangian is quadratic and the ODE that describes the relationship between
the position $q(t)$ and the control $u(t)$ is linear. We give a quantum
algorithm for this problem by first solving a closely related differential
equation called the Riccati equation. Since our quantum algorithm to solve the
LQR problem uses the Riccati equation, we describe this equation next. The
Riccati equation is intimately related to the optimal control problem as
explained in [48].

###  2.5 Riccati equation and flows on the Grassmannian

The Grassmannian (denoted Gr$(N,M)$) is the set of $M$ dimensional vector
spaces of an $N$ dimensional vector space (over the real or complex numbers).
While the field for these vector spaces can be arbitrary such as the
quaternions, we restrict attention to the complex numbers in this paper. The
Grassmannian has the structure of a smooth manifold since it can be viewed as
the homogeneous space GL$(N)/($GL$(N-M)\times$GL$(M))$. The action of the
general linear group on the Grassmannian is transitive and gives rise to a
vector bundle where the base manifold is the Grassmannian and the fibers are
isomorphic to GL$(M)$. This vector bundle is called the Stiefel manifold
$\mathbb{V}(N,M)$. It can be viewed as the set of $M$ frames i.e., the set of
tuples of $M$ linearly independent vectors in ${\mathbb{C}}^{N}$. This set of
points constitute the non-compact Stiefel manifold. The compact version of
this consists of the set of _orthonormal_ $M$ frames.

Any point in $\mathbb{V}(N,M)$ can be represented by an $N\times M$ matrix $y$
consisting of linearly independent columns. Two such points have the same span
if they differ by an element of the general linear group. This means that the
projection map $\pi:\mathbb{V}(N,M)\rightarrow$Gr$(N,M)$ takes $y$ to the
vector space consisting of the span of the columns of $y$. One can think of
the matrix $y$ as representing a point on the Grassmannian up to
multiplication by an element of GL$(M)$ i.e., $y$ and $\tilde{y}$ represent
the same point on Gr$(N,M)$ if $\tilde{y}=yg$, where $g\in$GL$(M)$. Strictly
speaking, the point $y$ lies in the Stiefel manifold and points in GL$(M)$ is
on the fiber above the point $\pi(y)$ on the Grassmannian.

In terms of their action on vectors $z$ of complex numbers in
${\mathbb{C}}^{N}$, the coordinates $y$ can be defined as

$z_{1}=y\cdot z_{2}\,,$ (2.53)

where $z_{1}$ is the set of the first $N-M$ components of $z$ and $z_{2}$ is
the rest of the $M$ components. Choosing these coordinates $y$ is equivalent
to choosing a set of $M$ independent vectors in $\mathbb{C}^{N}$. To view the
Riccati equation as a flow on $Gr(M,N)$, first consider a vector field on the
Steifel manifold given by

$z^{\prime}=Az\,,$ (2.54)

where $A$ is some $N\times N$ matrix and $z^{\prime}$ is the derivative of $z$
(with respect to some parameter $t$). Let $A$ have the block form

$A=\begin{pmatrix}F_{1}&F_{0}\\\ F_{2}&F_{3}\end{pmatrix}\,,$ (2.55)

where $F_{0}$ is $N-M\times M$, Then we have

$\displaystyle z_{1}^{\prime}=F_{1}z_{1}+F_{0}z_{2}\,,$ (2.56) $\displaystyle
z_{2}^{\prime}=F_{2}z_{1}+F_{3}z_{2}\,.$ (2.57)

Taking the derivative of (2.53) and using (2.56), we get

$F_{1}z_{1}+F_{0}z_{2}=z_{1}^{\prime}=y^{\prime}z_{2}+yz_{2}^{\prime}=y^{\prime}z_{2}+y(F_{2}z_{1}+F_{3}z_{2})\,.$
(2.58)

Using (2.53) again, we get

$y^{\prime}z_{2}=(F_{1}y+F_{0}-yF_{2}y-yF_{3})z_{2}\,.$ (2.59)

Since $z_{2}$ is arbitrary, we obtain the Riccati equation as a flow defined
on the Grassmannian. With this description of the Riccati equation, it becomes
easy to see a natural linearization method which we describe next.

###  2.6 Riccati equation and the Möbius transformation

Here, we introduce a well-known transformation called the Möbius
transformation which linearizes the Riccati equation. This is a slight
generalization of the one derived in [49] and [50]. This version is more
convenient for our quantum algorithm since it allows for a constant term.

######  Proposition 1.

Define the following initial value problem for a linear differential equation.

$\dot{x}=Ax+b\,,$ (2.60)

where

$\displaystyle x(t)=\begin{pmatrix}u(t)\\\ v(t)\end{pmatrix}\,,\hskip
21.68121ptA=\begin{pmatrix}F_{1}&F_{0}\\\ F_{2}&F_{3}\end{pmatrix}\,,\hskip
21.68121ptb=\begin{pmatrix}e\\\ f\end{pmatrix}\,,$ (2.61)

and where $u$ and $e$ are $N\times M$ matrices, $v$ and $f$ are $M\times M$
matrices. The initial conditions are $u(0)=y_{0}$ and $v(0)=I$. Then whenever
$v(t)^{-1}$ is defined, $y(t)=(u(t)+w)v(t)^{-1}$ is a solution to the matrix
Riccati equation (3.15), where $w$ can be picked to be any time-independent
vector and $e$ and $f$ are defined as

$e=F_{1}w\,\text{ and }\,f=F_{2}w\,.$ (2.62)

###### Proof.

From the definition of $y$ as

$y=(u+w)v^{-1}\,,$ (2.63)

we see that taking the derivative of $y$ gives us

$\dot{y}=\dot{u}v^{-1}-(u+w)v^{-1}\dot{v}v^{-1}\,.$ (2.64)

From (2.60) and (2.61), we have

$\displaystyle\dot{u}=F_{1}u+F_{0}v+e$ (2.65)
$\displaystyle\dot{v}=F_{2}u+F_{3}v+f\,.$ (2.66)

Replacing $\dot{u}$ and $\dot{v}$ in (2.64) above, we get

$\displaystyle\dot{y}$
$\displaystyle=(F_{1}u+F_{0}v+e)v^{-1}+(u+w)v^{-1}(F_{2}u+F_{3}v+f)v^{-1}$
(2.67) $\displaystyle=F_{0}+F_{1}y+yF_{2}y+yF_{3}\,,$ (2.68)

where we used the definitions of $e$ and $f$ from (2.62). The last equation is
the Riccati equation with the initial condition $y(0)=u(0)v(0)^{-1}=y_{0}$. ∎

When using the Möbius transformation to solve the Riccati equation, one can
run into singularities i.e., places where $v(t)$ is not invertible. This does
not mean that the solution to the Riccati equation does not exist since $y(t)$
is always defined. In [49], it was shown that this is due to the solution
crossing coordinate patches in the Grassmannian manifold. If one considers the
solution on the single patch, then there is no singularity in the solution. It
turns out that in our applications, the quantity $v(t)$ is invertible and in
fact, the condition number can be assumed to be well-behaved. As we explain
next, for the problem of finding the solution of an optimal control problem,
the fact that the Jacobi condition is satisfied is equivalent to $v$ in (2.61)
being non-singular.

###  2.7 Conjugate points, Jacobi’s condition and the Riccati equation

In this subsection, we describe the relationship between conjugate points and
the Riccati equation. Given a Lagrangian $L$, conjugate points are defined as
follows.

######  Definition 1.

A point $t=c$ is called a conjugate point in the interval $(a,b)$ if there
exists a non-trivial solution to the following boundary value problem in
$(a,c)$.

$\displaystyle\frac{d}{dt}[L_{\dot{q}\dot{q}}\dot{h}+L_{\dot{q}q}h]-L_{q\dot{q}}\dot{h}-L_{qq}h=0$
(2.69) $\displaystyle h(a)=h(c)=0\,,$ (2.70)

where $L(q,\dot{q},t)$ is the Lagrangian and terms such as $L_{\dot{q}q}$ are
its partial derivatives.

The Jacobi condition for an interval $(a,b)$ states that there are no
conjugate points in that interval. It turns out that there exists a Riccati
equation corresponding to the Jacobi equation [51]. It can be defined as
follows (via its linearization)

$\displaystyle\dot{U}(t)=AU(t)+BV(t)$ (2.71)
$\displaystyle\dot{V}(t)=CU(t)-A^{\dagger}V(t)\,,$ (2.72)

where $A,B$ and $C$ are defined as

$\displaystyle A=-L_{\dot{q}\dot{q}}^{-1}L_{\dot{q}q}$ (2.73) $\displaystyle
B=L_{\dot{q}\dot{q}}^{-1}$ (2.74) $\displaystyle
C=L_{qq}-L_{q\dot{q}}L_{\dot{q}\dot{q}}^{-1}L_{\dot{q}q}\,.$ (2.75)

The above is the matrix Riccati equation. Define the corresponding vector
Riccati equation

$\displaystyle\dot{u}(t)=Au(t)+Bv(t)$ (2.76)
$\displaystyle\dot{v}(t)=Cu(t)-A^{\dagger}v(t)\,,$ (2.77)

where $u$ and $v$ are vectors. By defining

$\displaystyle u=h$ (2.78) $\displaystyle
v=L_{\dot{q}\dot{q}}\dot{h}+L_{\dot{q}q}h\,,$ (2.79)

the Jacobi equation (2.69) is equivalent to the set of equations (2.71). To
see this, notice that the Jacobi equation is

$\dot{v}-L_{q\dot{q}}\dot{u}-L_{qq}u=0\,.$ (2.80)

Now Riccati equation can be written as

$\displaystyle v=B^{-1}(\dot{u}-Au)$ (2.81)
$\displaystyle\dot{v}=Cu-A^{\dagger}v\,,$ (2.82)

which leads to

$\dot{v}=(C-A^{\dagger}B^{-1}A)u-A^{\dagger}B^{-1}\dot{u}\,.$ (2.83)

Now using the definitions of $A,B$ and $C$, we have

$\displaystyle C-A^{\dagger}B^{-1}A=L_{qq}$ (2.84) $\displaystyle
A^{\dagger}B^{-1}=L_{q\dot{q}}\,,$ (2.85)

which shows that the Jacobi equation is equivalent to the Riccati equation.
The following lemma was proved in [48].

######  Lemma 1.

If the Jacobi condition is satisfied i.e., there are no conjugate points, then
$V(t)$ in the Riccati equation (2.71) is invertible.

In this paper, we assume a quantitative version of the Jacobi condition where
$V(t)$ is invertible and has a condition number bounded by $\kappa_{V}$.

###  2.8 From the Riccati equation to the Hamilton-Jacobi-Bellman equation

The Hamilton-Jacobi-Bellman equation is a nonlinear partial differential
equation of the form below.

$\frac{\partial V^{\ast}(x(t),t)}{\partial
t}=-\min_{u(t)}\Big{(}C(x(t),u(t),t)+\Big{[}\frac{\partial
V^{\ast}(x(t),t)}{\partial t}\Big{]}^{\dagger}f(x(t),u(t),t)\Big{)}\,.$ (2.86)

Here $C(x(t),u(t),t)$ is a cost function, which we assume is a quadratic
function of the following form.

$C(x(t),u(t),t)=u^{\dagger}(t)R(t)u(t)+x^{\dagger}(t)Q(t)x(t)\,,$ (2.87)

and that $V^{\ast}$ is of the form

$V^{\ast}(x(t),t)=x^{\dagger}(t)P(t)x(t)\,,$ (2.88)

and $f(x,u,t)$ has a linear form i.e.,

$f(x(t),u(t),t)=F(t)x(t)+G(t)u(t)\,.$ (2.89)

Using this and re-writing the HJB equation, we get

$x^{\dagger}\dot{P}x(t)=-\min_{u(t)}[u^{\dagger}Ru+x^{\dagger}Qx+2x^{\dagger}PFx+2x^{\dagger}PGu]\,.$
(2.90)

The next result can be found in [52].

######  Lemma 2.

Using the quadratic form for the cost function above and the linear form for
$f$, the solution to the HJB equation (2.86) can be obtained from a matrix
Riccati equation.

###### Proof.

Using the identity

$\displaystyle u^{\dagger}Ru+x^{\dagger}Qx+2x^{\dagger}PFx+2x^{\dagger}PGu=$
$\displaystyle(u+R^{-1}G^{\dagger}Px)^{\dagger}R(u+R^{-1}G^{\dagger}Px)+$
$\displaystyle x^{\dagger}(Q-PGR^{-1}G^{\dagger}P+PF+F^{\dagger}P)x\,.$ (2.91)

Since $R$ is positive definite and the second term above does not depend on
$u$, the minimum of (2.90) is obtained at

$\bar{u}(t)=-R^{-1}G^{\dagger}Px\,.$ (2.92)

Using this in (2.90), we get

$x^{\dagger}\dot{P}x(t)=-x^{\dagger}(Q-PGR^{-1}G^{\dagger}P+PF+F^{\dagger}P)x\,.$
(2.93)

Since this holds for all $x(t)$, we have

$\dot{P}=PGR^{-1}G^{\dagger}P-PF-F^{\dagger}P-Q\,.$ (2.94)

This is a matrix Riccati equation. Using its solution $P(t)$, we can get the
solution of the HJB equation by setting

$\bar{u}(t)=-R^{-1}G^{\dagger}Px(t)\,.$ (2.95)

∎

In applications to control problems such as linear quadratic control, one
needs to set a boundary condition rather than an initial condition i.e.,
$P(t_{f})=P_{f}$ and obtain a solution $P(t)$ for all time $t\leq t_{f}$.

##  3 Computational preliminaries and problem statements

###  3.1 Block encoding

We first define a useful concept called block-encoding (see [1] for more
details) and cite some results from the literature that we use later in the
paper.

######  Definition 2.

Suppose $A$ is an $a$-qubit operator, $\alpha,\epsilon\in\mathbb{R}_{+}$ and
$b\in\mathbb{N}$, then an $a+b$ qubit unitary $U$ is said to be an
$(\alpha,b,\epsilon)$ block encoding of $A$ if

$\|A-\alpha(\bra{0}^{b}\otimes I)U(\ket{0}^{b}\otimes I)\|\leq\epsilon\,.$
(3.1)

######  Definition 3.

We denote by $\bar{A}$ the Hermitian complement of $A$ defined below

$\bar{A}=\begin{pmatrix}0&A\\\ A^{\dagger}&0\end{pmatrix}\,.$ (3.2)

Next, we state the relevant results such as block-encoding a sparse matrix,
adding, multiplying and inverting block-encoded matrices.

######  Lemma 3 ([1], Lemma 48).

Suppose $A$ is an $s_{r}$-row-sparse and $s_{c}$-column-sparse $a$-qubit
matrix where each element of $A$ has absolute value at most 1 and where we
have oracle access to the positions and values of the nonzero entries, then we
have an implementation of $(\sqrt{s_{r}s_{c}},a+3,\epsilon)$ block encoding of
$A$ with a single use of $O_{r}$ and $O_{c}$ and two uses of $O_{A}$ and using
$O(a+\log^{2.5}(s_{r}s_{c}/\epsilon))$ elementary gates and
$O(b+\log^{2.5}(s_{r}s_{c}/\epsilon))$ ancillas.

The next result is from [8] and gives the complexity of inverting a matrix.

######  Lemma 4 ([8], Lemma 9).

Let $A$ be a matrix with condition number $\kappa\geq 2$. Let $H$ be the
Hermitian complement of $A$ and suppose that $I/\kappa\leq H\leq I$. Let

$\delta=o(\epsilon/\kappa^{2}\log^{3}(\kappa^{2}/\epsilon))\,.$ (3.3)

If $U$ is an $(\alpha,a,\delta)$ block encoding of $H$ that has gate
complexity $T_{U}$, then we can implement a

$(2\kappa,a+O(\log(\kappa^{2}\log(1/\epsilon))),\epsilon)$ (3.4)

block encoding of $H^{-1}$ with gate complexity

$O(\alpha\kappa(a+T_{U})\log^{2}(\kappa^{2}/\epsilon))\,.$ (3.5)

######  Remark 1.

In the above lemma, when $H$ is a Hermitian complement of a matrix $A$, then
after inversion, the resulting block encoded matrix satisfies

$\|A^{-1}-2\kappa(\bra{0}^{b}\otimes\bra{1}\otimes
I)U(\ket{0}^{b}\otimes\ket{0}\otimes I)\|\leq\epsilon\,,$ (3.6)

where the middle qubit corresponds to the qubit needed for the Hermitian
complement. By appending an $X$ gate to the block encoding for that qubit,
this can be brought into the above form of block encoding of $A^{-1}$ (as in
Definition 2).

The next result gives the complexity of performing matrix arithmetic on block-
encoded matrices.

######  Lemma 5 ([8], Lemma 6, [1], Lemmas 52, 54).

If $A$ has an $(\alpha,a,\epsilon)$ block encoding with gate complexity
$T_{A}$ and $B$ has a $(\beta,b,\delta)$ block encoding with gate complexity
$T_{B}$, then

  1. 1.

$\bar{A}$ has an $(\alpha,a+1,\epsilon)$ block encoding that can be
implemented with gate complexity $O(T_{A})$.

  2. 2.

$A+B$ has an $(\alpha+\beta,a+b,\beta\epsilon+\alpha\delta)$ block encoding
that can be implemented with gate complexity $O(T_{A}+T_{B})$.

  3. 3.

$AB$ has an $(\alpha\beta,a+b,\alpha\delta+\beta\epsilon)$ block encoding with
gate complexity $O(T_{A}+T_{B})$.

Next, we need the following theorem to implement the QLSA.

######  Theorem 1 ([16], Theorem 19).

Let $A$ be such that $\|A\|=1$ and $\|A^{-1}\|=\kappa$. Given a oracle block
encoding of $A$ and an oracle for implementing $\ket{b}$, there exists a
quantum algorithm which produces the normalized state $A^{-1}\ket{b}$ to
within an error $\epsilon$ using

$O(\kappa\log 1/\epsilon)$ (3.7)

calls to the oracles.

######  Remark 2.

In using Lemma 4 and Theorem 1, we will apply them to matrices with norms that
are $O(1)$ rather than $1$. This does not change the asymptotic scaling of the
gate complexities.

The next theorem from [8] gives a quantum algorithm to estimate the norm of
the solution state.

######  Theorem 2.

Let $Ax=b$ be an $N\times N$ linear system with an invertible matrix $A$ with
sparsity $s$ and condition number $\kappa$. Given an oracle that computes s
normalized version of $b$, there is a quantum algorithm that outputs a
quantity $\tilde{x}$ such that

$\Big{|}\tilde{x}-\|A^{-1}b\|\Big{|}\leq\gamma\|A^{-1}b\|\,,$ (3.8)

with constant probability and with gate complexity

$O(\frac{\kappa}{\gamma}(sT_{U}\log^{2}\frac{\kappa}{\gamma}+T_{b})\log^{3}\kappa\log\log\kappa)\,,$
(3.9)

where

$T_{U}=O(\log
N+\log^{2.5}(\frac{s\kappa}{\gamma}\log\frac{\kappa}{\gamma}))\,.$ (3.10)

Next we define the oracle access framework that is commonly assumed in quantum
algorithms in the literature. Suppose the matrix $A$ in linear ODEs such as
(2.45) is a block $2\times 2$ matrix composed of blocks $F_{0}$ through
$F_{3}$, we assume that the entries of these matrices are accessed through the
following oracles.

######  Definition 4.

We assume that each $F_{m}$ can be accessed through the following oracles.

$\displaystyle O_{r}\ket{i,k}=\ket{i,r_{ik}}\,,$ (3.11) $\displaystyle
O_{c}\ket{i,k}=\ket{i,c_{ik}}\,,$ (3.12) $\displaystyle
O_{a}\ket{i,j}\ket{0}^{\otimes b}=\ket{i,j,F_{m}(i,j)}\,,$ (3.13)

where $r_{ij}$ (respectively $c_{ji}$) is the $j^{th}$ nonzero entry in the
$i^{th}$ row (resp. column). If there are fewer than $j$ entries, then
$r_{ij}$ (resp. $c_{ji}$) is set to $j+2^{n}$. Here $F_{m}(i,j)$ is a binary
representation of the $(i,j)$ matrix entry of $F_{m}$. Here $F_{m}$ is one of
$F_{0}$, $F_{1}$, $F_{2}$ or $F_{3}$.

Finally, we define an oracle $O_{init}$ to prepare a normalized version of the
initial state $z_{0}$. We assume that the norm of $z_{0}$ is known.
Specifically, let $O_{init}$ be any unitary that maps $\ket{1}\ket{\phi}$ to
$\ket{1}\ket{\phi}$ for any state $\ket{\phi}$ and $\ket{0}\ket{0}$ to
$\ket{0}\ket{\bar{z}_{0}}$, where $\bar{z}_{0}=z_{0}/\|z_{0}\|$.

######  Remark 3.

One interesting and useful aspect of linearization of the Riccati equation
using a Möbius transformation is that we do not need to be able to prepare a
quantum state proportional to $F_{0}$ as was the case with Carleman
linearization [34]. We only need to access its entries via an oracle, which is
easier.

###  3.2 Problem statements

In this subsection, we define precisely the problems for which we have
efficient algorithms. The first problem is to simulate the evolution of a
classical physical system with a quadratic Hamiltonian with dissipation and
forcing terms.

######  Problem 1.

Given oracle access to matrices in the following equation of motion

$M\ddot{q}+R\dot{q}+Vq+f=0\,,$ (3.14)

and an oracle to prepare the normalized versions of the initial state $q(0)$
and $f$. We are also given estimates of their norms and promised that the
damping matrix $R$ is positive. The problem is to

  1. 1.

produce the normalized quantum state proportional to $(q(T),p(T))^{T}$ at some
time $T$, and

  2. 2.

output an estimate of the kinetic energy at time $T$ to within an additive
error $\epsilon$.

In optimal control, we are first interested in the following initial value
problem of the Riccati equation.

######  Problem 2.

Consider the following initial value problem for the vector Riccati equation

$\dot{y}=F_{0}+F_{1}y+yF_{2}y+yF_{3}\,,\hskip 28.45274pty(0)=y_{0}\,,$ (3.15)

where $y$ is an $N$ dimensional vector with initial value $y_{0}$ and the
coefficients $F_{0}$, $F_{1}$, $F_{2}$ and $F_{3}$ are time-independent. Given
oracle access to the entries of $F_{i}$, the problem is to output a quantum
state proportional to $y(T)$ for some $T$ to within a given $\ell_{2}$-norm
error $\epsilon$.

The next problem we consider is to find the solution of the matrix Riccati
equation as a block-encoded unitary.

######  Problem 3.

Given the matrix Riccati equation

$\dot{y}=F_{0}+F_{1}y+yF_{2}y+yF_{3}\,,\hskip 28.45274pty(0)=y_{0}\,,$ (3.16)

where $y$ is an $N\times M$ dimensional matrix with initial value matrix
$y_{0}$ and the coefficients $F_{0}$, $F_{1}$, $F_{2}$ and $F_{3}$ are time-
independent. Given oracle access to the entries of $F_{i}$ and $y_{0}$, the
problem is to output a block encoded unitary proportional to $y(T)$ for a
given $T$ to within a given $\ell_{2}$-norm error $\epsilon$.

###  3.3 Norm of the exponential

When implementing quantum algorithms to solve linear ODEs of the form

$\dot{z}=Az+b\,,$ (3.17)

It was shown in [25] that the gate complexity depends on the norm of the
exponential of the matrix $A$ i.e., it depends on

$C(A)=\sup_{t\in[0,T]}\|\exp(At)\|\,,$ (3.18)

where $T$ is the time to which we need to simulate the dynamics. Here, we
prove some bounds on this norm when $A$ is a block $2\times 2$ matrix as is
the case in all our applications here.

We first consider the case when $F_{0}=0$.

######  Lemma 6.

Suppose

$A=\begin{pmatrix}F_{1}&0\\\ F_{2}&F_{3}\end{pmatrix}\,.$ (3.19)

Then we have

$C(A)\leq C_{d}(1+C_{d}\|F_{2}\|T)\,,$ (3.20)

where

$C_{d}=\max\\{\sup_{t\in[0,T]}\|\exp(F_{1}t)\|,\sup_{t\in[0,T]}\|\exp(F_{3}t)\|\\}\,.$
(3.21)

###### Proof.

Split the matrix $A$ as follows

$A=\begin{pmatrix}F_{1}&0\\\ 0&F_{3}\end{pmatrix}+\begin{pmatrix}0&0\\\
F_{2}&0\end{pmatrix}=A_{1}+A_{2}\,,$ (3.22)

where $A_{1}$ contains the dissipative part of the evolution and $A_{2}$
contains the nonlinear term. From the definition in (3.21) above, we see that
$C(A_{1})=C_{d}$. Now using Dyson’s formula ([53] or [54]), we have

$\exp(At)=\exp(A_{1}t)+\int_{0}^{t}\exp(A_{1}(t-s))A_{2}\exp(As)ds\,.$ (3.23)

We can define a recursive version of this by using the formula to define
$\exp(As)$. We get (see [54])

$\exp(At)=\sum_{n=0}^{\infty}U_{n}(t)\,,$ (3.24)

where

$U_{n+1}(t)=\int_{0}^{t}\exp(A_{1}(t-s))A_{2}U_{n}(s)ds\,,$ (3.25)

and $U_{0}(t)=\exp(A_{1}t)$. It turns out that in our case, $U_{n}=0$ for
$n\geq 2$. To see this, note that

$\exp(A_{1}t)A_{2}\exp(A_{1}t^{\prime})A_{2}=0\,,$ (3.26)

since $A_{1}$ is a block diagonal matrix and $A_{2}$ is a strictly lower
triangular block matrix. Therefore, we have

$\exp(At)=\exp(A_{1}t)+\int_{0}^{t}\exp(A_{1}(t-s))A_{2}\exp(A_{1}s)ds\,.$
(3.27)

Taking norms, we get

$C(A)\leq C_{d}(1+C_{d}T\|F_{2}\|)\,,$ (3.28)

where we have used (3.21) and the fact that $\|A_{2}\|\leq\|F_{2}\|$. ∎

Note that we get a similar result when $F_{2}=0$ instead of $F_{0}$. We now
show a weaker result when both $F_{2}$ and $F_{0}$ are non-zero.

######  Lemma 7.

Assume that $F_{0}\neq 0$ and $F_{2}\neq 0$. Also assume that

$\|\exp(At)\|\leq\exp(\mu t)\,.$ (3.29)

Then we have

$\|\exp(At)\|\leq\exp((\mu+\|F_{2}\|+\|F_{0}\|)t)\,.$ (3.30)

###### Proof.

As before, split $A$ as follows

$A=\begin{pmatrix}F_{1}&0\\\ 0&F_{3}\end{pmatrix}+\begin{pmatrix}0&F_{0}\\\
F_{2}&0\end{pmatrix}=A_{1}+A_{2}\,,$ (3.31)

we have that $A_{1}$ contains the dissipative part of the evolution and
$A_{2}$ contains the nonlinear and constant terms. From [54], for any two
matrices $B$ and $E$ if $\|\exp(Bs)\|\leq\exp(\beta s)$, then

$\|\exp((B+E)s)\|\leq\exp((\beta+\|E\|)s)\,.$ (3.32)

We now use the fact that $\|\exp(A_{1}t)\|\leq\exp(\mu t)$, where $\mu$ is the
log norm of $A_{1}$. Using this, we see that

$\|\exp((A_{1}+A_{2})t)\|\leq\exp((\mu+\|F_{0}\|+\|F_{2}\|)t)\,.$ (3.33)

∎

######  Remark 4.

In order to have $C(A)$ be polynomially bounded i.e.,
$C(A)=\operatorname{poly}(n)$, we need

$\exp((\mu+\|F_{0}\|+\|F_{2}\|)t)=O(n^{k})\,,$ (3.34)

for some $k$ for all $t\in[0,T]$. If $\mu+\|F_{0}\|+\|F_{2}\|>0$, then we need

$(\mu+\|F_{0}\|+\|F_{2}\|)T=O(k\log n)\,.$ (3.35)

If $\mu+\|F_{0}\|+\|F_{2}\|<0$, then the condition in (3.34) is satisfied
already.

######  Remark 5.

We compare the above quantity to $R$, which was defined in [34] (see also [25]
for a generalized version) as follows.

$R=\frac{1}{|\mu|}\Big{(}\|F_{2}\|\|u_{in}\|+\frac{\|F_{0}\|}{\|u_{in}\|}\Big{)}\,,$
(3.36)

where $\mu$ is the log norm of $F_{1}$, which is assumed to be negative and
$\|u_{in}\|$ is the norm of the initial state. This means we need to re-define
$F_{0}$ and $F_{2}$ as $\bar{F_{0}}=F_{0}/\|u_{in}\|$ and
$\bar{F_{2}}=\|u_{in}\|F_{2}$. The requirement that $R<1$ (by Carleman
linearization) is equivalent to

$\mu+\|\bar{F_{2}}\|+\|\bar{F_{0}}\|\leq 0\,.$ (3.37)

In contrast, here we need (from (3.35))

$\mu+\|F_{2}\|+\|F_{0}\|=\frac{1}{T}O(\log n)\,,$ (3.38)

which grows logarithmically in $n$, which is the number of bits needed to
encode the problem i.e., $n=\log N$. The condition (3.38) is more general than
(and includes) (3.37). In many cases, it can lead to more nonlinearity and
even positive log-norm $\mu$.

######  Remark 6.

When $\mu(A_{1})<0$ and $F_{0}=0$, we have an even weaker constraint on the
nonlinearity. From (3.28), we have

$1+\|F_{2}\|T\leq\operatorname{poly}(n)\,,$ (3.39)

since $C_{d}\leq 1$. This means that the strength of the nonlinear term can be
as large as $\operatorname{poly}(n)/T$.

##  4 Quantum algorithms for classical Hamiltonians

In this section, we present our quantum algorithms to solve the problem of
preparing the final state after classical evolution and also to estimate the
kinetic energy as described in Problem 1. In this approach, we essentially
solve the differential equation coming from Hamilton’s canonical equations in
(2.45), namely an equation of the type

$\frac{dx}{dt}=Ax+b\,,$ (4.1)

using the quantum ODE solver recently developed in [25]. Specifically, we have
the following result.

######  Theorem 3.

There exists a quantum algorithm that solves Problem 1 of preparing a quantum
state given oracle access to $M$, $V$, $R$ and an oracle to prepare a
normalized version of $f$. Let $d$ be the dimension of these matrices and $s$
be the sparsity. Let

$k=O(\log T\|A\|(1+\frac{T\|f\|}{\|x(T)\|})\,,\hskip 7.22743pt\text{and}\hskip
7.22743pt\kappa_{L}=T\|A\|C(A)\,.$ (4.2)

The query complexity of this algorithm is given by

$T_{Q}=O(g\kappa_{L}\log\kappa_{L}\operatorname{poly}(s,k,\log d,\log
1/\epsilon))\,,$ (4.3)

and gate complexity $T_{G}$ is greater by a factor of at most

$O(\operatorname{poly}(k,\log(\frac{1}{\epsilon}),\log(\|A\|T))\,,$ (4.4)

where

$C(A)\leq\max\\{\|\sqrt{V}\|,\|\sqrt{M}\|\\}\cdot\max\\{\|\sqrt{V^{-1}}\|,\|\sqrt{M^{-1}}\|\\}\,,$
(4.5) $g=\frac{\max_{t\in[0,T]}\|x(t)\|}{\|x(T)\|}\,,$ (4.6)

where $x(t)$ is the solution at time $t$, $T$ is the final time and

$\|A\|\leq 2\max\\{\|M^{-1}V\|,\|M^{-1}R\|,1\\}\,.$ (4.7)

The query complexity of estimating the kinetic energy is

$T_{K}=O(\frac{1}{\epsilon}\kappa_{L}\operatorname{poly}(\log\kappa_{L},s,k,\log
d))\,,$ (4.8)

with the gate complexity greater by a fraction of at most
$\operatorname{poly}(\log 1/\epsilon,\log T\|A\|,k)$.

###### Proof.

Using the result from [25], we can obtain a quantum state $\ket{\psi}$ that is
close to the solution of the equation

$\frac{dx}{dt}=Ax+b\,,$ (4.9)

such that $\|\ket{\psi}-x(T)/\|x(T)\|\|\leq\epsilon$. To do so, we need a
block encoding of $A$ and an oracle to prepare the normalized version of $b$
and the normalized initial state. A block-encoding of $A$ can implemented by
first block-encoding the sparse matrices

$\begin{pmatrix}0&I\\\ -M^{-1}V&0\end{pmatrix}\,\hskip
7.22743pt\text{and}\hskip 7.22743pt\begin{pmatrix}0&0\\\
0&-M^{-1}R\end{pmatrix}\,,$ (4.10)

and then adding the two block-encoded unitaries. The initial state $x(0)$ in
this basis can be prepared easily using oracles to prepare $q(0)$ and
$\dot{q}(0)$ as normalized states similar to the description in [25]. The
oracle to prepare $b$ can be done using the method described in [33]. Recall
that

$b=\begin{pmatrix}0\\\ M^{-1}f\end{pmatrix}\,.$ (4.11)

The method in [33] which uses inequality testing from [55] can construct an
oracle for $b$ with error $\epsilon$ with
$O(\sqrt{\|M^{-1}\|}\operatorname{polylog}(d,1/\epsilon))$ queries to the
oracle for $f$.

Using this, the final state $\ket{\psi}$ can be obtained with query complexity

$O(gT\|A\|C(A)\log(\|A\|C(A)T)\operatorname{poly}(s,k,\log
d,\log(1/\epsilon)))\,,$ (4.12)

where

$C(A)=\sup_{t\in[0,T]}\exp(At)\,,$ (4.13)

and

$g=\frac{\max_{t\in[0,T]}\|x(t)\|}{\|x(T)\|}\,.$ (4.14)

Next, we estimate $C(A)$ in terms of the norms of $M$, $R$, $V$ and $f$.

##### Dissipation $R=0$.

First, as a special case, we consider the case when $R=0$. In this case, for
the classical mechanical system, we have from (2.29) that

$A=\begin{pmatrix}0&I\\\ -M^{-1}V&0\end{pmatrix}\,.$ (4.15)

The complexity of the algorithm depends on $C(A)$, which we now estimate. We
have

$A=I_{M}^{-1}JQI_{M}\,,$ (4.16)

where

$\displaystyle I_{M}=\begin{pmatrix}I&0\\\ 0&M\end{pmatrix}\,,$ $\displaystyle
Q=\begin{pmatrix}V&0\\\ 0&M^{-1}\end{pmatrix}\,.$ (4.17)

Therefore, we have

$A=I_{M}^{-1}Q^{-1/2}(Q^{1/2}JQ^{1/2})Q^{1/2}I_{M}\,.$ (4.18)

This means

$\exp(At)=(\sqrt{Q}I_{M})^{-1}\exp(Q^{1/2}JQ^{1/2}t)(\sqrt{Q}I_{M})\,.$ (4.19)

Taking norms and using the fact that $\sqrt{Q}J\sqrt{Q}$ is an anti-Hermitian
operator (as discussed in Section 2.2), we get

$\|\exp(At)\|\leq\|(\sqrt{Q}I_{M})^{-1}\|\|\exp(Q^{1/2}JQ^{1/2}t)\|\|\sqrt{Q}I_{M}\|=\kappa(\sqrt{Q}I_{M})\,,$
(4.20)

where in the above equation, we used the fact that the exponential of
$\sqrt{Q}J\sqrt{Q}$ has unit norm and where $\kappa(\sqrt{Q}I_{M})$ is the
condition number of $\sqrt{Q}I_{M}$. Finally, we have

$\sqrt{Q}I_{M}=\begin{pmatrix}\sqrt{V}&0\\\ 0&\sqrt{M}\end{pmatrix}\,.$ (4.21)

##### Dissipation $R\neq 0$.

In this case, we have that

$A=\begin{pmatrix}0&I\\\ -M^{-1}V&-M^{-1}R\end{pmatrix}\,.$ (4.22)

To bound $C(A)$ in this case, first we split $A$ as follows

$A=\begin{pmatrix}0&I\\\
-M^{-1}V&-M^{-1}R\end{pmatrix}=I_{M}^{-1}\Bigg{[}\begin{pmatrix}0&M^{-1}\\\
-V&0\end{pmatrix}+\begin{pmatrix}0&0\\\
0&-RM^{-1}\end{pmatrix}\Bigg{]}I_{M}\,.$ (4.23)

The first term inside the square brackets is just $JQ$ and let us denote the
second term $A_{2}$. Suppose we conjugate $A_{2}$ by $Q^{1/2}$, we get

$Q^{1/2}A_{2}Q^{-1/2}=\begin{pmatrix}0&0\\\
0&-M^{-1/2}RM^{-1/2}\end{pmatrix}\,.$ (4.24)

Let

$A_{R}=\begin{pmatrix}0&0\\\ 0&-M^{-1/2}RM^{-1/2}\end{pmatrix}\,,$ (4.25)

and similarly conjugating $JQ$, we get

$A_{I}=Q^{1/2}JQ^{1/2}$ (4.26)

Using this, we can see that

$A=I_{M}^{-1}Q^{-1/2}(Q^{1/2}JQ^{1/2})Q^{1/2}I_{M}+I_{M}^{-1}Q^{-1/2}A_{R}Q^{1/2}I_{M}\,.$
(4.27)

Therefore

$\exp(At)=(\sqrt{Q}I_{M})^{-1}\exp(A_{I}+A_{R})(\sqrt{Q}I_{M})\,.$ (4.28)

Hence

$\|\exp(At)\|\leq\kappa(\sqrt{Q}I_{M})\|\exp(A_{R}+A_{I})\|\,.$ (4.29)

Now a result from [53] states that for any matrix $B$

$\|\exp(B)\|\leq\|\exp((B+B^{\dagger})/2)\|\,.$ (4.30)

Letting $B=A_{R}+A_{I}$, we get $(B+B^{\dagger})/2=A_{R}$ (because $A_{I}$ is
anti-Hermitian). Therefore, we have

$C(A)\leq\kappa(I_{M}\sqrt{Q})C(A_{R})\,.$ (4.31)

Since we assume that $R$ is a positive matrix, we have $M^{-1}RM^{-1}$ is
positive and hence $A_{R}$ is negative. This means $\|\exp(A_{R}t)\|\leq 1$.
Therefore,

$C(A)\leq\kappa(\sqrt{Q}I_{M})\,.$ (4.32)

From the block structure of $A$, we can bound $\|A\|$ as follows. We have

$\|A\|=\sup_{v:\|v\|=1}\|Av\|\,.$ (4.33)

We can assume block structure for $v$ as $v=(v_{1},v_{2})^{T}$ such that
$\|v_{1}\|^{2}+\|v_{2}\|^{2}=1$. We now have

$\displaystyle\|A\|$
$\displaystyle=\sup_{(v_{1},v_{2})^{T}:\|v_{1}\|^{2}+\|v_{2}\|^{2}=1}\|(v_{2},-M^{-1}Vv_{1}-M^{-1}Rv_{2})^{T}\|$
(4.34)
$\displaystyle\leq\sup_{(v_{1},v_{2})^{T}:\|v_{1}\|^{2}+\|v_{2}\|^{2}=1}\sqrt{\|v_{2}\|^{2}+\max\\{\|M^{-1}V\|,\|M^{-1}R\|\\}^{2}\|v_{1}+v_{2}\|^{2}}$
(4.35)
$\displaystyle\leq\sup_{(v_{1},v_{2})^{T}:\|v_{1}\|^{2}+\|v_{2}\|^{2}=1}\sqrt{2}\max\\{\|M^{-1}V\|,\|M^{-1}R\|,1\\}\|v_{1}+v_{2}\|\,.$
(4.36)

We also have $\|v_{1}+v_{2}\|^{2}\leq 2(\|v_{1}\|^{2}+\|v_{2}\|^{2})$. Using
this, we can see that

$\|A\|\leq 2\max\\{\|M^{-1}V\|,\|M^{-1}R\|,1\\}\,.$ (4.37)

##### Estimation of kinetic energy.

The kinetic energy of the system at a given time is defined as

$K(t)=\frac{1}{2}\dot{q}(t)^{\dagger}M\dot{q}(t)\,.$ (4.38)

To estimate the kinetic energy of the final state, we need to use the so
called history state which is the quantum state produced by the algorithm
(from which we can get the normalized solution state by measuring the time
register). Here we will not measure the time register, but perform a Hadamard
test on the history state. This state can be written as

$\ket{\psi}=\frac{1}{N_{\psi}}\Bigg{(}\sum_{i=0}^{m}\ket{i,x_{i}}\Bigg{)}\,,$
(4.39)

where $N_{\psi}$ is the normalization and $m=T\|A\|$ denotes the number of
time steps needed to obtain an accurate solution (to within an additive
error). Using Theorem 2, we can estimate the norm of the above state to some
multiplicative error $\gamma$ i.e., there exists a quantum algorithm that
produces $\hat{N}_{\psi}$ such that

$|\hat{N}_{\psi}-N_{\psi}|\leq\gamma N_{\psi}\,.$ (4.40)

We can now measure the state $\ket{\psi}$ using the following procedure. We
first block-encode as $U_{M}$, the matrix

$\frac{1}{2}\ket{m}\bra{m}\otimes\begin{pmatrix}0&0\\\ 0&M\end{pmatrix}$
(4.41)

Then we need to estimate the inner product

$\frac{1}{2}\dot{q}^{\dagger}M\dot{q}=\frac{1}{2}\begin{pmatrix}q&\dot{q}\end{pmatrix}\begin{pmatrix}0&0\\\
0&M\end{pmatrix}\begin{pmatrix}q\\\ \dot{q}\end{pmatrix}\,,$ (4.42)

by estimating the inner product

$\bra{\psi}U_{M}\ket{\psi}=\frac{1}{N^{2}_{\psi}\|M\|}\frac{1}{2}\hat{\dot{q}}^{\dagger}M\hat{\dot{q}}\,,$
(4.43)

where $\hat{\dot{q}}$ is the velocity obtained from the numerical solution and
$\|M\|$ is the sub-normalization in the block-encoding. In other words, the
inner product is a scaled version of the kinetic energy. We can now use the
estimate of the norm to get the kinetic energy
$\hat{K}=\|M\|\hat{N}_{\psi}^{2}\bra{\psi}U_{M}\ket{\psi}$. The total error in
this estimate is

$|\hat{K}-K|=\frac{1}{2}\Bigg{|}\frac{\hat{N}^{2}_{\psi}}{N^{2}_{\psi}}\hat{\dot{q}}^{\dagger}M\hat{\dot{q}}-\dot{q}^{\dagger}M\dot{q}\Bigg{|}\,.$
(4.44)

This can be written as

$|\hat{K}-K|=\frac{1}{2}\Bigg{|}\frac{\hat{N}^{2}_{\psi}}{N^{2}_{\psi}}(\hat{\dot{q}}^{\dagger}M\hat{\dot{q}}-\dot{q}^{\dagger}M\dot{q})+(\frac{\hat{N}^{2}_{\psi}}{N^{2}_{\psi}}-1)\dot{q}^{\dagger}M\dot{q}\Bigg{|}\,.$
(4.45)

Using (4.40), we have $\hat{N}_{\psi}\leq(1+\gamma)N_{\psi}$ and this means

$|\hat{N}_{\psi}^{2}-N_{\psi}^{2}|=|\hat{N}_{\psi}-N_{\psi}|(\hat{N}_{\psi}+N_{\psi})\leq\gamma(2+\gamma)N_{\psi}^{2}\,.$
(4.46)

Therefore,

$|\hat{K}-K|\leq\frac{1}{2}(1+\gamma(2+\gamma))|(\hat{\dot{q}}^{\dagger}M\hat{\dot{q}}-\dot{q}^{\dagger}M\dot{q})|+\frac{1}{2}\gamma(2+\gamma)|\dot{q}^{\dagger}M\dot{q}|\,.$
(4.47)

Now let us focus on the difference of inner products in the first term. This
can be written as

$|\hat{\dot{q}}^{\dagger}M\hat{\dot{q}}-\dot{q}^{\dagger}M\dot{q}|=\sum_{i}m_{i}|\hat{\dot{q}}_{i}^{2}-\dot{q}_{i}^{2}|\leq
m_{\max}\sum_{i}|\hat{\dot{q}}_{i}^{2}-\dot{q}_{i}^{2}|\,,$ (4.48)

where $m_{\max}$ is the maximum value of the diagonal matrix of masses $M$.
Using Cauchy-Schwartz inequality, we can write this as

$m_{\max}\sum_{i}|\hat{\dot{q}}_{i}^{2}-\dot{q}_{i}^{2}|\leq
m_{\max}\sqrt{\sum_{i}(\hat{\dot{q}}_{i}-\dot{q}_{i})^{2}}\sqrt{\sum_{i}(\hat{\dot{q}}_{i}+\dot{q}_{i})^{2}}=m_{\max}\|\hat{\dot{q}}-\dot{q}\|\|\hat{\dot{q}}+\dot{q}\|\,,$
(4.49)

where $\|\cdot\|$ is the $\ell_{2}$ norm. Using the fact that

$\|\hat{\dot{q}}-\dot{q}\|\leq\epsilon_{1}\|\dot{q}\|\,,$ (4.50)

for some $\epsilon_{1}$ we get

$|\hat{\dot{q}}^{\dagger}M\hat{\dot{q}}-\dot{q}^{\dagger}M\dot{q}|\leq
m_{\max}\epsilon_{1}(2+\epsilon_{1})\|\dot{q}\|^{2}\,.$ (4.51)

Putting it together and choosing $\gamma$ such that
$\gamma=\epsilon_{1}=\epsilon/8$, we get

$|\hat{K}-K|\leq\epsilon K_{\max}\,,$ (4.52)

where

$K_{\max}=m_{\max}\|\dot{q}\|_{2}^{2}\,.$ (4.53)

The gate complexity of kinetic energy estimation follows from the complexity
of preparation of the state $\ket{\psi}$, norm estimation of $N_{\psi}$, block
encoding $U_{M}$ and the Hadamard test. ∎

######  Remark 7.

The run-time for the case with zero damping i.e., when $R=0$, depends on the
condition number of the matrix $\sqrt{Q}$. The fact that for the encoding used
here, the run-time scales as the inverse of the smallest eigenvalue is already
pointed out in [33], where it was argued that it is at least linear. Here we
show that it is no more than linear (up to log factors).

It is interesting to note that the problem shown to be ${\mathsf{BQP}}$
complete in [33] also has a constant condition number for $\sqrt{Q}$ and hence
the above algorithm and encoding would be efficient for this case.

######  Remark 8.

It is also worth noting that the run-time of the quantum algorithm from [25]
for the linear differential equation (4.9) is dependent on the norm of the
exponential of $A$ (denoted $C(A)$ here) and does not require $A$ to have a
negative log-norm.

Next we give a result that improves the above kinetic energy estimation by
giving a multiplicative or relative approximation to the kinetic energy if
there is no source term i.e., $b=0$.

######  Corollary 1.

If $b=0$ in (4.9), then there exists a quantum algorithm that outputs an
estimate $\hat{K}$ for the kinetic energy such that

$|K-\hat{K}|\leq\epsilon K\,,$ (4.54)

where $K$ is the actual kinetic energy.

###### Proof.

From (4.47) in the proof of the $b\neq 0$ case, we have that the estimate
$\hat{K}$ satisfies

$|K-\hat{K}|\leq\frac{1}{2}(1+\gamma(2+\gamma))\sum_{i}m_{i}|\dot{q}_{i}^{2}-\hat{\dot{q}}_{i}^{2}|+\frac{1}{2}\gamma(2+\gamma)|\dot{q}^{\dagger}M\dot{q}|\,.$
(4.55)

We show in Lemma 13 in Appendix A that for any $i$

$|\dot{q}_{i}-\hat{\dot{q}}_{i}|\leq\epsilon_{1}\dot{q}_{i}\,,$ (4.56)

for some constant $\epsilon_{1}$. Using this, we get

$|\dot{q}_{i}^{2}-\hat{\dot{q}}_{i}^{2}|=|\dot{q}_{i}-\hat{\dot{q}}_{i}||\dot{q}_{i}+\hat{\dot{q}}_{i}|\leq\epsilon_{1}(2+\epsilon_{1})\dot{q}_{i}^{2}\,,$
(4.57)

which gives

$|K-\hat{K}|\leq(\epsilon_{1}(2+\epsilon_{1})(1+\gamma(2+\gamma))+\gamma(2+\gamma))K\,.$
(4.58)

By choosing $\gamma=\epsilon_{1}=\epsilon/16$, we get

$|K-\hat{K}|\leq\epsilon K\,.$ (4.59)

∎

###  4.1 Improved kinetic energy estimation

In this subsection, we explain how to use Quantum Phase Estimation (QPE) and
Qubitization [7] and the improvements to this implemented in [56, 57]. In this
approach, we have an algorithm whose complexity does not have a dependence on
the condition number of $\sqrt{Q}$ as with the previous approach. The downside
is that it has a worse dependence on the error. Recall that in the $y$ basis,
we have the differential equation given by (2.42) reproduced below.

$\frac{dy}{dt}=\begin{pmatrix}0&\sqrt{VM^{-1}}\\\
-\sqrt{M^{-1}V}&-\sqrt{M^{-1}}R\sqrt{M^{-1}}\end{pmatrix}y+b=\tilde{A}y+b\,.$
(4.60)

In order to solve this equation, we need to be able to block-encode
$\sqrt{M^{-1}V}$. This can be done using the following result shown in [33].

######  Theorem 4.

Given oracle access to the sparse matrix $A$, one can implement a block-
encoding of $\sqrt{A}$ with query complexity

$T_{Q}(A)=O(\|A\|_{\max}s\log(1/\epsilon)\frac{1}{\epsilon}\min(\sqrt{A^{-1}},\frac{1}{\epsilon}))\,,$
(4.61)

and gate complexity greater by a factor of at most

$O(\operatorname{polylog}(\frac{d}{\epsilon}))\,.$ (4.62)

Using this result, we prove the following theorem.

######  Theorem 5.

There exists a quantum algorithm that solves Problem 1 of estimating to within
an additive error of $\epsilon$, the kinetic energy of a classical mechanical
system described by its equations of motion as follows

$M\ddot{q}+R\dot{q}+Vq+f=0\,,$ (4.63)

We assume that we are given oracle access to $M$, $V$, $R$ and an oracle to
prepare a normalized version of $f$. Let $d$ be the dimension of these
matrices and $s$ be the sparsity. Let

$k=O(\log T\|\tilde{A}\|(1+\frac{T\|f\|}{\|x(T)\|})\,,\hskip
7.22743ptm=T\|\tilde{A}\|\,.$ (4.64)

The query complexity of this algorithm is given by

$T_{Q}=O(\frac{m}{\epsilon^{2}}\log(m)\operatorname{poly}(s,k,\log d))\,,$
(4.65)

and gate complexity $T_{G}$ is greater by a factor of at most

$O(\operatorname{poly}(k,\log(\frac{1}{\epsilon}),\log(m))\,,$ (4.66)

where

$\|\tilde{A}\|\leq 2\max\\{\|M^{-1}V\|,\|M^{-1}R\|,1\\}\,.$ (4.67)

###### Proof.

We can first block encode $\sqrt{M^{-1}V}$ as a unitary $U_{1}$ using Theorem
4. Using this block-encoding, we can block encode the following matrix

$\begin{pmatrix}0&\sqrt{VM^{-1}}\\\ -\sqrt{M^{-1}V}&0\end{pmatrix}\,.$ (4.68)

This can be implemented by the following circuit

$U_{2}=cU_{1}^{\dagger}(iY\otimes I)cU_{1}\,,$ (4.69)

where $Y$ is the Pauli matrix on a single qubit (also denoted $\sigma_{y}$)
and $cU_{1}$ is the controlled version of the unitary $U_{1}$ controlled by
the first qubit. The identity in the tensor product is of the same dimension
as $U_{1}$. Next, we block encode the sparse matrix

$\begin{pmatrix}0&0\\\ 0&-\sqrt{M^{-1}}R\sqrt{M^{-1}}\end{pmatrix}\,,$ (4.70)

as $U_{3}$. Finally, we use block-encoded matrix arithmetic to perform
$U_{2}+U_{3}$. To extract the kinetic energy, we first write it as follows.
Initial state preparation and the oracle to prepare $b$ can be done as before
using the techniques from [33]. Recall that $y$ is

$y=\begin{pmatrix}\sqrt{V}q\\\ \sqrt{M^{-1}}p\end{pmatrix}\,.$ (4.71)

Therefore, the kinetic energy can be written as

$K=\frac{1}{2}\dot{q}^{T}M\dot{q}=\frac{1}{2}y^{T}Py\,,$ (4.72)

where

$P=\begin{pmatrix}0&0\\\ 0&I\end{pmatrix}\,,$ (4.73)

which is the projector onto the second block of $y$. Similar to the previous
approach, we can implement a block-encoding (denoted $U_{P}$) of $P$ tensored
with the time register as follows (one can also implement a unitary that puts
a sign on the basis states corresponding to $P$).

$\ket{m}\bra{m}\otimes P\,.$ (4.74)

In a manner similar to the proof of Theorem 3, we prepare the history state
that has the unnormalized approximate solution $\hat{y}$ at time $T$. This
state can be written as

$\ket{\phi}=\frac{1}{N_{\phi}}\sum_{i=0}^{m}\ket{i,y_{i}}\,,$ (4.75)

and the query complexity of preparing this state is

$O(T\|\tilde{A}\|C(\tilde{A})\log(\|\tilde{A}\|C(\tilde{A})T)\operatorname{poly}(s,k,\log
d))\,,$ (4.76)

where

$C(\tilde{A})=\sup_{t\in[0,T]}\exp(\tilde{A}t)\,.$ (4.77)

We now estimate the quantity $C(\tilde{A})$ in this case. First, we can write
$\tilde{A}$ as

$\tilde{A}=\begin{pmatrix}0&\sqrt{VM^{-1}}\\\
-\sqrt{M^{-1}V}&-\sqrt{M^{-1}}R\sqrt{M^{-1}}\end{pmatrix}=\begin{pmatrix}0&\sqrt{VM^{-1}}\\\
-\sqrt{M^{-1}V}&0\end{pmatrix}+\begin{pmatrix}0&0\\\
0&-\sqrt{M^{-1}}R\sqrt{M^{-1}}\end{pmatrix}\,.$ (4.78)

We denote the terms in the above sum $\tilde{A}_{I}$ and $\tilde{A}_{R}$. Now
using the fact that

$\|\exp(\tilde{A}t)\|\leq\|\exp(\frac{\tilde{A}+\tilde{A}^{\dagger}}{2}t)\|\,,$
(4.79)

and the fact that the first term in the sum on the right hand side of (4.78)
i.e., $\tilde{A}_{I}$ is anti-Hermitian, we get

$\|\exp(\tilde{A}t)\|\leq\|\exp(\tilde{A}_{R}t)\|\,.$ (4.80)

Now since $\tilde{A}_{R}$ is negative definite, we have
$\|\exp(\tilde{A}t)\|\leq 1$, which makes $C(\tilde{A})\leq 1$.

Next, we estimate the norm $N_{\phi}$ using Theorem 2. Finally, we perform the
Hadamard test to estimate the inner product of $\ket{\phi}$ with $U_{P}$. This
produces the following estimate of the normalized kinetic energy.

$\hat{K}=\frac{N_{\phi}}{2}\bra{\phi}U_{P}\ket{\phi}\,.$ (4.81)

The error analysis is similar to the one in the proof of Theorem 3. The main
difference is that the block-encoding of $\sqrt{M^{-1}V}$ leads to a query
complexity that depends on $1/\epsilon_{PE}$, where $\epsilon_{PE}$ is the
phase estimation error. Choosing $\epsilon_{PE}=\epsilon$, we get the
quadratic dependence of the run-time on the error. The overall query
complexity is as stated in the theorem. ∎

###  4.2 Coupled oscillators with damping

In this section, we focus specifically on coupled oscillators with damping and
forcing. Using the specific structure of the matrices in this problem, in
[33], a different square-root of $Q$ for this system was shown to be
efficiently implementable as a block-encoding. There, it was shown that there
exists a matrix $B$ which is $d\times d^{2}$ such that
$B^{\dagger}B=\sqrt{M^{-1}}V\sqrt{M^{-1}}$ and that this matrix can be block-
encoded without the overhead of $1/\epsilon$ that phase-estimation requires.
Here we show that one can include damping and solve this problem using
differential equation solvers. Going back to the $x$ basis, we have the
equation

$\frac{dx}{dt}=Ax+b\,,$ (4.82)

where

$A=\begin{pmatrix}0&I\\\ -M^{-1}V&-M^{-1}R\end{pmatrix}\,.$ (4.83)

We now define the basis, which we denote $\tilde{x}$, as follows.

$\tilde{x}=\begin{pmatrix}\sqrt{M}&0\\\ 0&\sqrt{M}\end{pmatrix}x\,.$ (4.84)

In this basis, the differential equation becomes

$\frac{d\tilde{x}}{dt}=\begin{pmatrix}\sqrt{M}&0\\\
0&\sqrt{M}\end{pmatrix}\begin{pmatrix}0&I\\\
-M^{-1}V&-M^{-1}R\end{pmatrix}\begin{pmatrix}\sqrt{M^{-1}}&0\\\
0&\sqrt{M^{-1}}\end{pmatrix}\tilde{x}+\begin{pmatrix}\sqrt{M}&0\\\
0&\sqrt{M}\end{pmatrix}b\,.$ (4.85)

Since $b=(0,\sqrt{M}f)^{T}$, we have

$\frac{d\tilde{x}}{dt}=\begin{pmatrix}0&I\\\
-M^{-1}VM^{-1}&-M^{-1}RM^{-1}\end{pmatrix}\tilde{x}-\begin{pmatrix}0\\\
\sqrt{M^{-1}}f\end{pmatrix}\,.$ (4.86)

We will write this as

$\frac{d\tilde{x}}{dt}=(J\tilde{Q}+\tilde{A}_{R})\tilde{x}+\tilde{b}\,,$
(4.87)

where $J$ was defined before and

$\tilde{Q}=\begin{pmatrix}M^{-1}VM^{-1}&0\\\ 0&I\end{pmatrix}\,,$ (4.88)

and

$\tilde{A}_{R}=\begin{pmatrix}0&0\\\ 0&-M^{-1}RM^{-1}\end{pmatrix}\,.$ (4.89)

Now define the matrix $B$ as follows (see [33]).

$BB^{\dagger}=\sqrt{M^{-1}}V\sqrt{M^{-1}}\,,$ (4.90)

where $B$ is a $d\times d(d+1)/2$ dimensional matrix (recall that $V$ is a
$d\times d$ matrix). Let us denote $d_{1}=d(d+1)/2$ for brevity. Now define
the $2d\times(d_{1}+d)$ matrix $\tilde{B}$ as follows.

$\tilde{B}=\begin{pmatrix}B_{d\times d_{1}}&0_{d\times d}\\\ 0_{d\times
d_{1}}&I_{d\times d}\end{pmatrix}\,.$ (4.91)

Note that we have

$\tilde{B}\tilde{B}^{\dagger}=\begin{pmatrix}BB^{\dagger}&0_{d\times d}\\\
0_{d\times d}&I_{d\times d}\end{pmatrix}\,,$ (4.92)

which is exactly $\tilde{Q}$. Now define the basis $\tilde{y}$ as
$\tilde{B}^{\dagger}\tilde{x}$. In this basis, we have

$\frac{d\tilde{y}}{dt}=\tilde{B}^{\dagger}(J\tilde{B}\tilde{B}^{\dagger}+\tilde{A}_{R})\tilde{x}+\tilde{B}^{\dagger}\tilde{b}\,.$
(4.93)

It is easy to check that

$\tilde{B}^{\dagger}\tilde{A}_{R}=\tilde{A}_{R}=\tilde{A}_{R}\tilde{Q}=\tilde{A}_{R}\tilde{B}^{\dagger}\tilde{B}\,.$
(4.94)

Using this, we get

$\frac{d\tilde{y}}{dt}=(\tilde{B}^{\dagger}J\tilde{B}+\tilde{A}_{R})\tilde{y}+\tilde{B}^{\dagger}\tilde{b}\,,$
(4.95)

where

$\tilde{B}^{\dagger}\tilde{b}=\begin{pmatrix}0\\\
B^{\dagger}\sqrt{M^{-1}}f\end{pmatrix}\,.$ (4.96)

Now that we have the differential equation in the $\tilde{y}$ basis, we can
prove the following result. In the following, rather than make queries to a
oracle for $V$, we assume access to an oracle for spring constants $S$ as in
[33] (where it is denoted $K$) as well as an oracle for $M$.

######  Theorem 6.

When we restrict Problem 1 to the specific case of damped coupled oscillators
with forcing, there exists a quantum algorithm that produces an estimate of
the kinetic energy at time $T$ to within an additive error $\epsilon$ that has
query complexity

$T_{Q}=O(\frac{m}{\epsilon}\operatorname{poly}(k,\log\frac{1}{\epsilon},\log
m))$ (4.97)

and gate complexity greater by a factor of at most

$\operatorname{polylog}(m,\frac{1}{\epsilon},d,\frac{M_{max}}{M_{\min}})$
(4.98)

where

$m=Ts\max\bigg{(}\frac{R_{\max}}{M_{\min}},\sqrt{\frac{S_{\max}}{M_{\min}}}\bigg{)}\,,$
(4.99)

and

$k=O(\log m(1+\frac{T\|f\|}{\|\tilde{y}(T)\|})\,.$ (4.100)

###### Proof.

The main idea is to solve the differential equation for $\tilde{y}$ and obtain
the history state. From the history state, we can get the kinetic energy by
block-encoding the projector onto the velocity space tensored with the final
time i.e.,

$U_{P}=\bra{m}\ket{m}\otimes P_{V}\,,$ (4.101)

where $P_{V}$ is the projector on to the velocity component of $\tilde{y}$.
The differential equation can be solved using the algorithm of [25] provided
some conditions are satisfied. First, we need a block-encoding of the matrix

$\tilde{B}^{\dagger}J\tilde{B}+\tilde{A}_{R}=\begin{pmatrix}0&B^{\dagger}\\\
-B&-\sqrt{M^{-1}}R\sqrt{M^{-1}}\end{pmatrix}\,.$ (4.102)

This can be done by block-encoding $\tilde{B}^{\dagger}J\tilde{B}$ and
$\tilde{A}_{R}$ separately and adding them using matrix arithmetic. The block-
encoding of $\tilde{B}^{\dagger}J\tilde{B}$ is described in [33] and takes
$O(\log d+\log^{2}(M_{\max}/M_{\min}\epsilon))$ and provides a block-encoding
with parameters

$(\sqrt{2sS_{\max}/M_{\min}},2\log d+\log 1/\epsilon,\epsilon)\,.$ (4.103)

The initial state preparation and an oracle for $\sqrt{M^{-1}}f$ is also
similar to [33]. Finally, we need to compute the norm of the matrix
exponential of the matrix in (4.102). We can see that this matrix has an anti-
Hermitian and a Hermitian components. The anti-Hermitian component is
$\tilde{B}^{\dagger}J\tilde{B}$ and the Hermitian part is $\tilde{A}_{R}$. As
in the previous sections, we can bound the norm of the exponential by

$\exp((\tilde{B}^{\dagger}J\tilde{B}+\tilde{A}_{R})t)\leq\exp(\tilde{A}_{R}t)\leq
1\,,$ (4.104)

since $R$ is a positive matrix. This gives the stated complexity. The analysis
of the error obtained in the kinetic energy estimation is the same as before.
∎

###  4.3 Hardness of kinetic energy estimation

In this subsection, we prove the main hardness result that estimating the
kinetic energy of a $d$ dimensional coupled oscillator system governed by
(3.14) i.e.,

$M\ddot{q}+R\dot{q}+Vq=0\,,$ (4.105)

is ${\mathsf{BQP}}$ hard when the damping strength (or the norm of $R$) is
bounded by an inverse polynomial in $\log d$. To prove this result, we need
the following result proved in [33], which shows the following. Here we give
an equivalent formulation of the hardness result.

######  Theorem 7.

There exists a family of instances of coupled oscillators such that it is
${\mathsf{BQP}}$ hard to estimate the normalized kinetic energy $K/E$ to
within additive precision $\epsilon$, where $K$ is the kinetic energy and $E$
is the total energy.

######  Remark 9.

The above formulation follows from the fact that in [33], it was shown that
the normalized kinetic energy is equal to the output probability of a quantum
circuit and estimating the latter to within additive precision is
${\mathsf{BQP}}$ hard.

######  Remark 10.

It was shown in [33] that the instances can be chosen such that the masses are
all equal to identity and the matrix of coupling constants $V$ has a condition
number independent of the number of oscillators. The initial conditions can
also be chosen to be $q(0)=(0,\dots,0)^{T}$ and
$\dot{q}(0)=(1,-1,0,\dots,0)^{T}$.

When the masses are all unity, the kinetic energy can be written as

$K=\frac{1}{2}\dot{q}^{\dagger}\dot{q}\,,$ (4.106)

where $\dot{q}$ is given by

$\dot{q}(t)=\cos(\sqrt{V}t)\dot{q}(0)\,.$ (4.107)

Before we prove our result, we need the following lemma which derives a
formula for the exponential of a matrix of the form $JQ$, where $Q$ is block
diagonal.

######  Lemma 8.

Let

$Q=\begin{pmatrix}V&0\\\ 0&I\end{pmatrix}\,,$ (4.108)

and (as before)

$J=\begin{pmatrix}0&I\\\ -I&0\end{pmatrix}\,.$ (4.109)

Then we have

$\exp(JQt)=\begin{pmatrix}\cos(\sqrt{V}t)&\sqrt{V^{-1}}\sin(\sqrt{V}t)\\\
-\sqrt{V}\sin(\sqrt{V}t)&\cos(\sqrt{V}t)\end{pmatrix}\,.$ (4.110)

###### Proof.

Exponentiating $JQt$, we have

$\exp(JQt)=\sum_{k=0}^{\infty}\frac{(JQt)^{k}}{k!}\,.$ (4.111)

To calculate this, observe that

$(JQ)^{2}=\begin{pmatrix}0&I\\\ -V&0\end{pmatrix}\begin{pmatrix}0&I\\\
-V&0\end{pmatrix}=\begin{pmatrix}-V&0\\\ 0&-V\end{pmatrix}=I\otimes(-V)\,,$
(4.112)

where $I$ is a $2\times 2$ identity matrix. Now consider one of the terms of
the right hand side of $(\ref{eq:exp})$. When $k=2\ell$ is even, we have

$(JQt)^{k}=t^{2\ell}(-1)^{\ell}(I\otimes V^{\ell})\,.$ (4.113)

Therefore, the even terms in (4.111) are

$\sum_{\ell=0}^{\infty}\frac{t^{2\ell}(-1)^{\ell}}{(2\ell)!}(I\otimes\sqrt{V}^{2\ell})=I\otimes\cos(\sqrt{V}t)\,.$
(4.114)

Similarly, the odd terms can be calculated to be

$(JQt)^{2\ell+1}=JQ(JQ)^{2\ell}=JQt^{2\ell+1}(-1)^{\ell}(I\otimes\sqrt{V}^{2\ell})\,.$
(4.115)

This can be written as

$t^{2\ell+1}(-1)^{\ell}J\begin{pmatrix}\sqrt{V}&0\\\
0&\sqrt{V^{-1}}\end{pmatrix}(I\otimes\sqrt{V}^{2\ell+1})\,.$ (4.116)

Therefore, the odd terms can be written as

$J\begin{pmatrix}\sqrt{V}&0\\\
0&\sqrt{V^{-1}}\end{pmatrix}\sum_{\ell=0}^{\infty}\frac{t^{2\ell+1}(-1)^{\ell}}{(2\ell+1)!}(I\otimes\sqrt{V}^{2\ell+1})=J\begin{pmatrix}\sqrt{V}&0\\\
0&\sqrt{V^{-1}}\end{pmatrix}(I\otimes\sin(\sqrt{V}t))\,.$ (4.117)

Putting them together, we get

$\exp(JQt)=\begin{pmatrix}\cos(\sqrt{V}t)&\sqrt{V^{-1}}\sin(\sqrt{V}t)\\\
-\sqrt{V}\sin(\sqrt{V}t)&\cos(\sqrt{V}t)\end{pmatrix}\,.$ (4.118)

∎

From the above lemma, we see that for the initial state given by
$q(0)=(0,\dots,0)^{T}$ and $\dot{q}(0)=(1,-1,0,\dots,0)^{T}$, we have

$\begin{pmatrix}q(t)\\\
\dot{q}(t)\end{pmatrix}=\exp(JQt)\begin{pmatrix}q(0)\\\
\dot{q}(0)\end{pmatrix}\,,$ (4.119)

where

$Q=\begin{pmatrix}V&0\\\ 0&I\end{pmatrix}\,.$ (4.120)

Using Lemma 8, we can write

$\dot{q}(t)=\cos(\sqrt{V}t)\dot{q}(0)\,.$ (4.121)

We can now prove our main hardness result.

######  Theorem 8.

For the classical mechanical system defined in Problem 1, when we restrict to
damped coupled oscillators, let $K(T)$ denote the kinetic energy at time $T$
and $E$ denote the total energy at time $t=0$. There exists a family of
instances such that it is $\mathsf{BQP}$ hard to estimate the normalized
kinetic energy $K(T)/E$ to within an additive error $\epsilon$ when
$T=\operatorname{polylog}(d)$ and $\|R\|T\leq\epsilon$ (where $d$ is the
dimension of the system).

###### Proof.

When all the mass terms are unity, the kinetic energy can be written as

$K(t)=\frac{1}{2}\dot{q}^{T}(t)\dot{q}(t)\,,$ (4.122)

which can be written as

$K(t)=\frac{1}{2}x(t)^{T}Px(t)\,,$ (4.123)

where $x(t)=(q(t),\dot{q}(t))^{T}$ and $P$ is the projector on the subspace
spanned by $\dot{q}(t)$. When $R=0$, we also have that

$x(t)=\exp(JQt)x(0)\,,$ (4.124)

which gives

$K(t)=\frac{1}{2}x(0)^{T}\exp(-QJt)P\exp(JQt)x(0)\,,$ (4.125)

using the fact that $J^{T}=-J$ and $Q^{T}=Q$. Similarly, for non-zero $R$, we
have

$K_{R}(t)=\frac{1}{2}x(0)^{T}\exp(-QJt+A_{R}t)P\exp(JQt+A_{R}t)x(0)\,,$
(4.126)

where

$A_{R}=\begin{pmatrix}0&0\\\ 0&-R\end{pmatrix}\,.$ (4.127)

Now using Dyson’s formula, we have

$\exp(JQt+A_{R}t)=\exp(JQt)+\int_{0}^{t}\exp(JQ(t-s))A_{R}\exp(JQs+A_{R}s)ds\,,$
(4.128)

and similarly for $\exp(-QJt+A_{R}t)$. Plugging this into (4.126) above, we
get

$\displaystyle 2K_{R}(t)$ (4.129)
$\displaystyle=x(0)^{T}\exp(-QJt)P\exp(JQt)x(0)$ (4.130)
$\displaystyle+\int_{0}^{t}x(0)^{T}\exp(-QJt)P\exp(JQ(t-s))A_{R}\exp(JQs+A_{R}s)x(0)ds$
(4.131)
$\displaystyle+\int_{0}^{t}x(0)^{T}\exp(-QJ(t-s))A_{R}\exp(-QJs+A_{R}s)P\exp(JQt)x(0)ds$
(4.132)
$\displaystyle+\int_{0}^{t}\int_{0}^{t}x(0)^{T}\exp(-QJ(t-s_{1}))A_{R}\exp(-QJs_{1}+A_{R}s_{1})P$
$\displaystyle\exp(JQ(t-s_{2})A_{R}\exp(JQs_{2}+A_{R}s_{2})x(0)ds_{1}ds_{2}\,.$
(4.133)

The first term (4.130) can be seen to be $K(t)$. The second term (4.131) is

$\displaystyle\int_{0}^{t}x(0)^{T}\exp(-QJt)P\exp(JQ(t-s))A_{R}\exp(JQs+A_{R}s)x(0)ds$
(4.134)
$\displaystyle=\int_{0}^{t}x(0)^{T}\sqrt{Q}\exp(-\sqrt{Q}J\sqrt{Q}t)\sqrt{Q^{-1}}P\sqrt{Q^{-1}}\exp(\sqrt{Q}J\sqrt{Q}(t-s))$
$\displaystyle\sqrt{Q}A_{R}\sqrt{Q^{-1}}\exp(\sqrt{Q}J\sqrt{Q}s+\sqrt{Q}A_{R}\sqrt{Q^{-1}}s)\sqrt{Q}x(0)ds\,.$
(4.135)

Note that since

$P=\begin{pmatrix}0&0\\\ 0&I\end{pmatrix}\,,\hskip
7.22743ptQ=\begin{pmatrix}V&0\\\ 0&I\end{pmatrix}\,\hskip
7.22743pt\text{and}\hskip 7.22743ptA_{R}=\begin{pmatrix}0&0\\\
0&-R\end{pmatrix}\,,$ (4.136)

we have

$\sqrt{Q^{-1}}P\sqrt{Q^{-1}}=P\,,\hskip 7.22743pt\text{and}\hskip
7.22743pt\sqrt{Q}A_{R}\sqrt{Q^{-1}}=A_{R}\,.$ (4.137)

The above equation gives us

$\|\exp(\sqrt{Q}J\sqrt{Q}s+\sqrt{Q}A_{R}\sqrt{Q^{-1}}s)\|\leq\|\exp(A_{R}s)\|\leq
1\,,$ (4.138)

where we have used the fact that $\sqrt{Q}J\sqrt{Q}$ is anti-Hermitian and
that $A_{R}$ is negative definite. All of this gives us that the norm of the
second term is

$\|x(0)\sqrt{Q}\|^{2}\|R\|t\,.$ (4.139)

Now, we have that

$\|x(0)\sqrt{Q}\|^{2}=x(0)^{T}Qx(0)=q^{T}Vq+\dot{q}^{T}\dot{q}=2E\,,$ (4.140)

where $E$ is the total energy of the classical mechanical system. The norm of
the third term (4.132) and the fourth term (4.133) can be bounded in a similar
way and we get

$|K-K_{R}|\leq 2E\|R\|t+E\|R\|^{2}t^{2}\,.$ (4.141)

Therefore

$\frac{1}{E}|K-K_{R}|\leq 2\|R\|t(1+\frac{1}{2}\|R\|t)\,.$ (4.142)

Suppose there is a quantum algorithm that can compute an estimate $\hat{K}$
such that

$|\hat{K}-\frac{K_{R}}{E}|\leq\epsilon\,,$ (4.143)

Then $\hat{K}$ also satisfies

$|\hat{K}-\frac{K}{E}|\leq|\hat{K}-\frac{K_{R}}{E}|+\frac{1}{E}|K-K_{R}|\leq\epsilon+2\|R\|t(1+\frac{1}{2}\|R\|t)\,,$
(4.144)

Therefore, if

$2\|R\|t(1+\frac{1}{2}\|R\|t)\leq\epsilon\,,$ (4.145)

i.e., if $\|R\|t\leq\epsilon/4$, then we have

$|\hat{K}-\frac{K}{E}|\leq 2\epsilon\,.$ (4.146)

Thus the estimate $\hat{K}$ also gives an additive approximation to the
kinetic energy when $R=0$, which is a ${\mathsf{BQP}}$ hard problem by Theorem
7. ∎

##  5 Quantum algorithms for optimal control

In this section, we present quantum algorithms for the vector and matrix
Riccati equations. Like our quantum algorithm for the estimation of kinetic
energy of damped mechanical systems, our quantum algorithm for the vector
Riccati equation highlights the fact that it is sometimes useful to have the
history state that gives access to the unnormalized solution state. Note that
in this section, we reuse previous variables such as $x$, $y$, $A$ and $b$
with new definitions.

###  5.1 Vector Riccati equation

Recall the vector Riccati equation from (3.15) i.e.,

$\dot{y}=F_{0}+F_{1}y+yF_{2}y+yF_{3}\,,\hskip 28.45274pty(0)=y_{0}\,,$ (5.1)

where $y$ and $F_{0}$ are $d$ dimensional column vectors, $F_{1}$ is a
$d\times d$ matrix, $F_{2}$ is a $d$ dimensional row vector and $F_{3}$ is a
complex number. Suppose the linear differential equation obtained by the
Möbius transformation from Section 2.6 is

$\dot{x}=Ax+b\,,\hskip 72.26999ptx(0)=x_{0}\,,$ (5.2)

where

$A=\begin{pmatrix}F_{1}&F_{0}\\\ F_{2}&F_{3}\end{pmatrix}\,,$ (5.3)

and

$b=\begin{pmatrix}F_{1}w\\\ F_{2}w\end{pmatrix}\,,$ (5.4)

where $w$ is any time-independent vector. Here

$x=(u,v)^{T}\,,$ (5.5)

is such that $y=(u+w)v^{-1}$. We now show how to construct oracles to access
the entries of $A$ and prepare a state proportional to $b$.

######  Lemma 9.

Using the oracles for $F_{i}$, $i=0,1,2,3$, we can construct an oracle $O_{A}$
and $O_{b}$ with a constant overhead in gate complexity.

###### Proof.

To get $O_{A}$, we can use the $F_{i}$ and apply them conditioned on the row
and column index. Specifically, $O_{A}$ is

$O_{A}=\sum_{i,j=1}^{d-1}\ket{i,j}\bra{i,j}\otimes
F_{1}+\sum_{i=1}^{d-1}\ket{i,d}\bra{i,d}\otimes
F_{0}+\sum_{j=1}^{d-1}\ket{d,j}\bra{d,j}\otimes
F_{2}+\ket{d,d}\bra{d,d}\otimes F_{3}\,.$ (5.6)

To create $O_{b}$, we can take $\ket{w}$ to be the all-zeros basis state
$\ket{0}$ and therefore

$b=\begin{pmatrix}e\\\ f\end{pmatrix}\,,$ (5.7)

where the $d$ dimensional vector $e$ and the number $f$ are

$e(i)=F_{1}(i,1)\,,\text{ and }\,f=F_{2}(1)\,.$ (5.8)

Since $F_{i}$ are assumed to be sparse with a constant number of nonzero
entries in every row or column, $O_{b}$ can be constructed using standard
techniques with constant gate complexity. ∎

To solve the above equation, we will use the result from [25] stated below.
Note that this version is to produce the history state and not the normalized
solution state.

######  Theorem 9.

Given an instance of a linear inhomogeneous equation (5.2) where $A$ has
sparsity $s$, dimension $d$, $C(A)=\sup_{t\in[0,T]}\exp(At)$ and oracles as
given in Definition 4, there exists a quantum algorithm that produces a
quantum state proportional to

$\ket{z}\sim\sum_{i=0}^{m-1}\ket{i,x_{i}}+\sum_{i=m}^{2m-1}\ket{i,x_{m}}\,,$
(5.9)

where $m=T\|A\|$ is the number of time steps, $x_{i}$ are $\epsilon$
approximate solutions at time steps $i=0,\dots,m$. Let

$k=O(\log T\|A\|(1+\frac{T\|b\|}{\|x(T)\|})\,.$ (5.10)

This algorithm succeeds with a constant success probability and with query
complexity

$O\Big{(}mC(A)\operatorname{poly}\Big{(}\log
d,s,k,\log(\frac{1}{\epsilon}),\log(mC(A))\Big{)}\Big{)}\,,$ (5.11)

and the gate complexity is greater by a factor of at most

$O\Big{(}\operatorname{poly}\Big{(}k,\log(\frac{1}{\epsilon}),\log(m)\Big{)}\Big{)}\,.$
(5.12)

The quantum algorithm to solve the Riccati equation (3.15) consists of the
following steps. In the steps below, for any unnormalized vector $\phi$, the
normalized version is denoted by $\ket{\phi}$.

  1. 1.

We take three registers, one for the time step, one for the Taylor series
level and one for the solution of the differential equation (and other states
such as the initial state). Define $z_{in}$ as

$z_{in}=\ket{0,0,x_{0}}+h\sum_{i=1}^{m}\ket{i,1,b}+\sum_{i=m+1}^{2m}\ket{i,0,w}\,,$
(5.13)

where $\ket{w}=\ket{0}$, $x_{0}=[y_{0},1]^{T}$ (with $\ket{x_{0}}$ being the
normalized version of $x_{0}$) and where $y_{0}$ is the initial state of the
Riccati equation in (3.15) and $h=1/\|A\|$. Let us denote

$\psi_{in}=\ket{0,0,x_{0}}+h\sum_{i=1}^{m}\ket{i,1,b}\,,$ (5.14)

and

$w_{in}=\sum_{i=m+1}^{2m}\ket{i,0,w}\,.$ (5.15)

In the first step, we construct the state $\ket{z_{in}}$.

  2. 2.

Next, construct a linear system for the ODE (5.2) from [25] i.e.,

$Lz=z_{in}\,,$ (5.16)

and solve the linear system to obtain the final state $\ket{z}$. This
unnormalized state can be written as

$z=L^{-1}z_{in}=L^{-1}\psi_{in}+L^{-1}w_{in}\,.$ (5.17)

Note that $w_{in}=Lw_{in}$ since $L$ acts as the identity when the first
register is in the support of the states $m+1$ to $2m$. Therefore

$z=L^{-1}\psi_{in}+w_{in}=\sum_{i=1}^{m}\ket{i,0,x_{i}}+\sum_{i=m+1}^{2m}\ket{i,0,x_{m}}+\sum_{i=m+1}^{2m}\ket{i,0,w}\,.$
(5.18)

  3. 3.

Measure the time register and post-select outcomes $m+1,\dots,2m$ to get a
state proportional to

$\psi_{T}=x_{m}+w\,.$ (5.19)

  4. 4.

Finally, perform the two outcome measurement of the computational basis with
outcomes $P=\\{1,\dots,d\\}$ and $I-P=\\{d+1\\}$. Post-select on the outcome
$k\in\\{1,\dots,d\\}$ (recall that $x_{m}$ is $d+1$ dimensional) to get the
state $\ket{\hat{y}}$, which is a normalized version of the approximate
solution $\hat{y}$. We show in the next section that $\hat{y}$ is the
approximate version of the actual solution $y(t)$ of equation (3.15). We can
write this state as

$\ket{\hat{y}}=\frac{P\psi_{T}}{\|P\psi_{T}\|}\,,$ (5.20)

where $P$ denotes the orthogonal projector onto the first $d$ components of
the $d+1$ dimensional state $\psi_{T}$.

Next we analyze all the aspects of the quantum algorithm for the Riccati
equation such as its accuracy, success probability and the gate and query
complexity of its circuit implementation.

####  5.1.1 Solution error

In this subsection, we bound the norm of the error between the normalized
solution and the normalized approximate solution i.e,

$\|\ket{y}-\ket{\hat{y}}\|=\Bigg{\|}\frac{P(x+w)}{\|P(x+w)\|}-\frac{P(\psi_{T})}{\|P(\psi_{T})\|}\Bigg{\|}\,,$
(5.21)

where $P$ is the projector onto the components $1,\dots,d$ of the $d+1$
dimensional vectors $x_{m}$, $w$ and $x$. We first bound the error between
unnormalized versions of these states in the lemma below. First recall that
$\|w\|=1$, $P(x+w)=u$ and $(I-P)(x+w)=v$, where $u$ and $v$ are from (2.61).
We use this notation below.

######  Lemma 10.

Suppose that the relative error in estimating $x$ is $\epsilon_{1}$ i.e.,

$\|x-\hat{x}\|\leq\epsilon_{1}\|x\|\,,$ (5.22)

then the error between $y$ and $\hat{y}$ can be bounded as

$\|y-\hat{y}\|\leq\epsilon\,,$ (5.23)

where

$\epsilon=2\epsilon_{1}(1+\frac{1}{\|y\|}+\frac{1}{\|u\|})\,,$ (5.24)

###### Proof.

First consider

$\displaystyle\|P(x+w)-P(\hat{x}+w)\|$ (5.25) $\displaystyle=\|Px-P\hat{x}\|$
(5.26) $\displaystyle\leq\epsilon_{1}\|x\|$ (5.27)
$\displaystyle=\epsilon_{1}\|(x+w)-w\|$ (5.28)
$\displaystyle\leq\epsilon_{1}\Bigg{(}1+\frac{1}{\|y\|}+\frac{1}{\|u\|}\Bigg{)}\|u\|\,,$
(5.29)

where in the last line, we used the fact that $\|w\|=1$ and $y=uv^{-1}$.

Using this, we can now bound $\|\ket{y}-\ket{\hat{y}}\|$ as follows.

$\displaystyle\|\ket{y}-\ket{\hat{y}}\|$
$\displaystyle=\Bigg{\|}\frac{P(x+w)}{\|P(x+w)\|}-\frac{P(\psi_{T})}{\|P(\psi_{T})\|}\Bigg{\|}$
(5.30)
$\displaystyle=\Bigg{\|}\frac{P(x+w)}{\|P(x+w)\|}-\frac{P(\hat{x}+w)}{\|P(\hat{x}+w)\|}\Bigg{\|}$
(5.31)
$\displaystyle\leq\Bigg{\|}\frac{P(x+w)}{\|P(x+w)\|}-\frac{P(\hat{x}+w)}{\|P(x+w)\|}\Bigg{\|}+\Bigg{\|}\frac{P(\hat{x}+w)}{\|P(x+w)\|}-\frac{P(\hat{x}+w)}{\|P(\hat{x}+w)\|}\Bigg{\|}$
(5.32)
$\displaystyle\leq\frac{\epsilon}{2}+\|P(\hat{x}+w)\|\frac{\Big{|}\|P(\hat{x}+w)\|-\|P(x+w)\|\Big{|}}{\|P(\hat{x}+w)\|\|P(x+w)\|}$
(5.33) $\displaystyle=\epsilon\,.$ (5.34)

∎

####  5.1.2 Probability of success

######  Lemma 11.

The probability of obtaining the state $\hat{y}/\|y\|$, where $\hat{y}$ is the
approximate solution to the Riccati equation is bounded below as follows.

$p_{success}\geq\frac{1}{108g^{2}}\,,$ (5.35)

where

$g=\frac{\max_{t\in[0,T]}\|u(t)\|^{2}}{\|u(T)\|^{2}}\,.$ (5.36)

###### Proof.

The probability of getting $\hat{y}$ can be broken down into two parts. The
first part is the probability of getting the state $\psi_{T}$ and the second
is the projector onto the first $d$ components. The overall probability can be
written as follows.

$p_{\mathrm{success}}=\frac{m\|x_{m}+w\|^{2}}{\|z\|^{2}}\frac{\|P(x_{m}+w)\|^{2}}{\|x_{m}+w\|^{2}}=\frac{m\|P(x_{m}+w)\|^{2}}{\|z\|^{2}}\,,$
(5.37)

where we use $x_{m}$ for $\hat{x}$ from (5.18). We now upper bound the
denominator. Recall from (5.18) that

$\|z\|^{2}\leq\sum_{i=1}^{m}\|x_{i}\|^{2}+m\|x_{m}+w\|^{2}\,.$ (5.38)

Using Theorem 8 from [25], we can bound the above quantity in terms of the
maximum norm of the solution. Using this, we have

$\displaystyle\|z\|^{2}$ $\displaystyle\leq
m(1+\gamma)^{2}\max_{t\in[0,T]}\|x(t)\|^{2}+m(1+\gamma)^{2}\|x(T)+w\|^{2}$
(5.39) $\displaystyle\leq
4m(1+\gamma)^{2}\Big{[}\max_{t\in[0,T]}\|x(t)\|^{2}+1\Big{]}$ (5.40)
$\displaystyle\leq 12m(1+\gamma)^{2}u_{max}^{2}\,.$ (5.41)

The probability can now be bounded as (assuming $\epsilon,\gamma<1/2$)

$\displaystyle p_{\mathrm{success}}$
$\displaystyle\geq\frac{m\|P(x_{m}+w)\|^{2}}{12m(1+\gamma)^{2}u_{max}^{2}}$
(5.42)
$\displaystyle\geq\frac{(1-\epsilon)^{2}\|u(T)\|^{2}}{12(1+\gamma)^{2}u_{max}^{2}}$
(5.43) $\displaystyle\geq\frac{1}{108g^{2}}\,,$ (5.44)

where

$g=\frac{\max_{t\in[0,T]}\|u(t)\|^{2}}{\|u(T)\|^{2}}\,.$ (5.45)

∎

####  5.1.3 Main result

Collecting these results together, we have our main result for Riccati
equations. As pointed out in Section 3.3, the result below gives a quantum
algorithm for the Riccati nonlinear equation in the highly nonlinear regime.

######  Theorem 10.

For the problem of producing a quantum state proportional to the solution of
the Riccati differential equation described in Problem 2, the quantum
algorithm for the Riccati equation produces a quantum state proportional to
the solution to within additive error $\epsilon$ with a constant success
probability with gate complexity given by

$O\Big{(}gT\|A\|C(A)sk\operatorname{polylog}\Big{(}d,k,\frac{1}{\epsilon},T\|A\|C(A)\Big{)}\Big{)}\,,$
(5.46)

and the query complexity is

$O\Big{(}gT\|A\|C(A)sk\operatorname{polylog}\Big{(}k,\frac{1}{\epsilon},T\|A\|C(A)\Big{)}\Big{)}\,,$
(5.47)

where

$\displaystyle C(A)\leq\begin{cases}C_{d}(1+C_{d}T\|F_{2}\|),&\text{if
$F_{0}=0$ and $F_{2}\neq 0$}\\\ C_{d}(1+C_{d}T\|F_{0}\|),&\text{if $F_{2}=0$
and $F_{0}\neq 0$}\\\ \exp(\mu+\|F_{0}\|+\|F_{2}\|),&\text{if both
nonzero}\end{cases}$ (5.48)
$\displaystyle\|A\|\leq\max(\|F_{0}\|+\|F_{1}\|,\|F_{2}\|+\|F_{3}\|)\,,$
(5.49)

where $\mu$ is the log-norm of $F_{1}$ and $C_{d}$ is defined in (3.21) and

$k=O(\log(1+\frac{\|b\|Te^{2}}{\|x(T)\|}))\,.$ (5.50)

###### Proof.

We discuss the construction of the initial state first. In a manner similar to
[24, 25], we first construct the state

$\frac{\|x_{0}\|}{N_{\mathrm{init}}}\ket{0,0,0,0}+\frac{\sqrt{m}h}{N_{\mathrm{init}}}\ket{0,1,0,0}+\frac{\sqrt{m}}{N_{\mathrm{init}}}\ket{0,0,1,0}\,,$
(5.51)

where

$N_{\mathrm{init}}=\sqrt{\|x_{0}\|^{2}+mh^{2}\|b\|^{2}+m}\,.$ (5.52)

Next, controlled on the second and third registers being in $00$ and $10$
respectively, we apply $O_{x}$ and $O_{b}$ to get

$\frac{\|x_{0}\|}{N_{\mathrm{init}}}\ket{0,0,0,\bar{x}_{0}}+\frac{\sqrt{m}h}{N_{\mathrm{init}}}\ket{0,1,0,\bar{b}}+\frac{\sqrt{m}}{N_{\mathrm{init}}}\ket{0,0,1,0}\,.$
(5.53)

Now conditioned on the second and third registers being in $10$ apply the
rotation on the first register that takes

$\ket{0}\rightarrow\frac{1}{\sqrt{m}}\sum_{i=0}^{m-1}\ket{i}\,,$ (5.54)

and conditioned on the second and third registers being in $01$, apply the
rotation

$\ket{0}\rightarrow\frac{1}{\sqrt{m}}\sum_{i=m}^{2m-1}\ket{i}\,.$ (5.55)

This produces the initial state

$\ket{z_{0}}=\frac{1}{N_{\mathrm{init}}}\Big{[}\ket{0,0,x_{0}}+h\sum_{i=1}^{m-1}\ket{i,1,b}+\sum_{i=m+1}^{2m}\ket{i,0,w}\Big{]}\,.$
(5.56)

The next step of the algorithm is to implement the quantum algorithm described
in [25]. The implementation complexity of that algorithm is proved in [25].
The bound on $\|A\|$ comes from its block structure and the bound on $C(A)$
comes from Lemma 6 and Lemma 7. This gives the query and gate complexities in
the theorem. ∎

###  5.2 Matrix Riccati equation

In this subsection, we give an algorithm to implement the solution of the
matrix Riccati equation described in (3.15). We assume that we have a block
encoded version of the linearized matrix $A$ from (2.60). This can be obtained
from oracle access to the matrices $F_{i}$ using Lemma 9. Our goal is to
create a block encoding of the solution $y(t)$, which is an $N\times M$
matrix. The algorithm consists of the following steps.

  1. 1.

Following [25], we first create a block encoded unitary that encodes the
linear system that produces the solution of the linear ODE from (2.60) and
invert the linear system to give a unitary $U_{1}$ (which is a
$(\alpha_{1},a_{1},\epsilon_{1})$ block encoding of $L^{-1}$). More precisely,
we create the block encoded matrix $L$ from [25], which is given by

$L=I-N\,,$ (5.57)

where

$\displaystyle N=M_{1}(I-M_{2})^{-1}$ (5.58) $\displaystyle
M_{1}=\sum_{j=0}^{k-1}\ket{j+1}\bra{j}\otimes\frac{Ah}{j+1}$ (5.59)
$\displaystyle M_{2}=\sum_{j=0}^{k}\ket{0}\bra{j}\otimes I\,.$ (5.60)

It was shown in [25] that $L^{-1}$ acts on the initial state $z^{\prime}_{in}$
and produces a state $z^{\prime}$ that encodes the solution of the linear ODE
i.e.,

$L^{-1}z^{\prime}_{in}=z^{\prime}\,,$ (5.61)

where

$z^{\prime}_{in}=\ket{0,0,x_{0}}+h\sum_{i=0}^{m-1}\ket{i,1,b^{\prime}}\,.$
(5.62)

In the equation above and in [25], $b^{\prime}$ is a vector which can be
prepared efficiently as a normalized quantum state.

  2. 2.

In our case, the initial state $x_{0}$ and the constant term in the linear ODE
$b$ are $(N+M)\times M$ matrices. Recall that $b$ is the matrix

$b=\begin{pmatrix}e\\\ f\end{pmatrix}=\begin{pmatrix}F_{1}w\\\
F_{2}w\end{pmatrix}\,,$ (5.63)

where $w$ can be chosen to be any time-independent matrix of the same
dimension as $x_{0}$.

To be able to use the algorithm from [25], we encode these as block unitary
matrices as follows i.e., we want to block-encode the matrix

$z_{in}=\ket{0,0}\bra{0,0}\otimes
x_{0}+h\sum_{i=0}^{m-1}\ket{i,1}\bra{0,0}\otimes b\,.$ (5.64)

Every column of $z_{in}$ when acted upon by $L^{-1}$ gives a column of the
solution matrix. To implement $z_{in}$ as a block unitary, we can use sparse
matrix encoding to first block-encode

$\ket{0,0}\bra{0,0}\otimes x_{0}+\sqrt{m}h\ket{0,1}\bra{0,0}\otimes b\,.$
(5.65)

Next we implement a unitary on the left taking

$\ket{0,0}\rightarrow\ket{0,0}\hskip 7.22743pt\text{and}\hskip
7.22743pt\ket{0,1}\rightarrow\frac{1}{\sqrt{m}}\sum_{i=0}^{m-1}\ket{i,1}\,.$
(5.66)

Denote this block encoded unitary $U_{2}$ with parameters
$(\alpha_{2},a_{2},\epsilon_{2})$.

  3. 3.

Next we use the block-encoding of $L^{-1}$ and a block encoding of $z_{in}$ to
get a block encoding of the solution i.e., $L^{-1}z_{in}$. Denote by $U_{3}$
the block encoding of $L^{-1}z_{in}$ obtained by multiplying the block
encodings of $L^{-1}$ and $z_{in}$. Let $U_{3}$ be a
$(\alpha_{3},a_{3},\epsilon_{3})$ block encoding of $L^{-1}z_{in}$.

  4. 4.

We now have a unitary that block encodes the solution $x_{T}$. This can be
written as

$x_{T}=\begin{pmatrix}u(T)\\\ v(T)\end{pmatrix}\,.$ (5.67)

Using matrix arithmetic, we can block encode $x_{T}+w$ as $U_{4}$ with
parameters $(\alpha_{4},a_{4},\epsilon_{4})$.

  5. 5.

Now conditioned on the lower block, we can invert $v$ and create a block
unitary $U_{5}$. Assume that it is a $(\alpha_{5},a_{5},\epsilon_{5})$ block
encoding of $v^{-1}$.

  6. 6.

Finally, we multiply $U_{4}$ and $U_{5}$ to create a block version
$U_{6}=U_{4}U_{5}$ of the solution $y_{T}=(u+w)v^{-1}$ with parameters
$(\alpha_{6},a_{6},\epsilon_{6})$.

We compute the block-encoding triples and analyze the algorithm below.

####  5.2.1 Analysis of the algorithm

In this subsection, we show that the algorithm above produces a solution for
the matrix Riccati equation. Specifically, we prove the following.

######  Theorem 11.

For the matrix Riccati equation defined in (3.15), let $\kappa_{V}$ be the
condition number of $V(t)$ coming from the Jacobi condition. Then, the above
algorithm produces a matrix encoded as a block unitary, which is $\epsilon$
close to the solution with query and gate complexity given by

$O(s\kappa_{V}\kappa_{L}\operatorname{polylog}(\frac{1}{\epsilon},s,\kappa_{V},\kappa_{L},M+N))\,.$
(5.68)

where

$\displaystyle C(A)\leq\begin{cases}C_{d}(1+C_{d}T\|F_{2}\|),&\text{if
$F_{0}=0$ and $F_{2}\neq 0$}\\\ C_{d}(1+C_{d}T\|F_{0}\|),&\text{if $F_{2}=0$
and $F_{0}\neq 0$}\\\ \exp(\mu+\|F_{0}\|+\|F_{2}\|),&\text{if both
nonzero}\end{cases}$ (5.69)
$\displaystyle\|A\|\leq\max(\|F_{0}\|+\|F_{1}\|,\|F_{2}\|+\|F_{3}\|)\,,$
(5.70)

where the log-norm of $F_{1}$ is denoted $\mu$ and

$C_{d}=\max\\{\sup_{t\in[0,T]}\|\exp(F_{1}t)\|,\sup_{t\in[0,T]}\|\exp(F_{3}t)\|\\}\,,$
(5.71)

and

$\kappa_{L}=O(C(A)T\|A\|)\,.$ (5.72)

###### Proof.

We now bound the query and gate complexities of the algorithm to produce a
state $\epsilon$ close to the solution. From [25], we have that the condition
number of $L$ is $\kappa_{L}=O(C(A)m)$, where $m=T\|A\|$.

  1. 1.

To implement $U_{1}$, which is a block encoding of $L^{-1}$, the query and
gate complexities are given by Lemma 4. If we block encode $L$ as a
$(s,\log(m(N+M)),\delta)$ block encoding, then $L^{-1}$ is a
$(\alpha_{1},a_{1},\epsilon_{1})$ block encoding, where

$\displaystyle\alpha_{1}=2\kappa_{L}$ (5.73) $\displaystyle
a_{1}=\log(m(N+M))+O(\log(\kappa_{L}^{2}\log 1/\epsilon))\,.$ (5.74)

We also need to pick $\delta$ according to Lemma 4 as

$\delta=o(\epsilon/\kappa_{L}^{2}\log^{3}(\kappa_{L}^{2}/\epsilon))\,.$ (5.75)

The gate complexity is

$T_{1}=O(s\kappa_{L}(\log(m(N+M))+T_{U})\log^{2}(\kappa_{L}^{2}/\epsilon))\,,$
(5.76)

where the cost of block encoding $L$ is

$T_{U}=O(\log(m(N+M))+\log^{2.5}(s^{2}/\delta))\,.$ (5.77)

  2. 2.

To implement $U_{2}$ as a $(\alpha_{2},a_{2},\epsilon_{2})$ block encoding of
$z_{in}$, we can use Lemma 3. The parameters are

$\displaystyle\alpha_{2}=s$ (5.78) $\displaystyle a_{2}=\log(m(N+M))\,.$
(5.79)

The query complexity is $O(1)$ and the gate complexity is
$O(\log(m(N+M))+\log^{2.5}(s^{2}/\epsilon_{2}))$.

  3. 3.

Multiplying two block encoded unitaries can be done using Lemma 5. The
complexity of implementing $U_{3}$ is $O(T_{1}+T_{2})$ and the parameters are
$(\alpha_{3},a_{3},\epsilon_{3})$, where

$\displaystyle\alpha_{3}=\alpha_{1}\alpha_{2}$ (5.80) $\displaystyle
a_{3}=a_{1}+a_{2}$ (5.81)
$\displaystyle\epsilon_{3}=\epsilon_{1}\alpha_{2}+\epsilon_{2}\alpha_{1}\,.$
(5.82)

  4. 4.

This step can be done using matrix arithmetic. A block encoding of $w$ with
parameters $(1,\log m(N+M),\epsilon)$ can be done in
$T_{w}=O(\log(m(N+M))+\log^{2.5}(1/\epsilon))$ gates. Therefore, the gate
complexity of block encoding $u+w$ is $T_{4}=O(T_{3}+T_{w})$.

  5. 5.

Inverting $v$ can be done using Lemma 4. We can get a
$(\alpha_{5},a_{5},\epsilon_{5})$ block encoding of $v^{-1}$ where

$\displaystyle\alpha_{5}=2\kappa_{V}$ (5.83) $\displaystyle a_{5}=O(\log
M+\log(\kappa_{V}^{2}\log(1/\epsilon_{5})))\,.$ (5.84)

The gate complexity is

$T_{5}=O(\alpha_{4}\kappa_{V}(\log
M+T_{4})\log^{2}(\kappa_{V}^{2}/\epsilon))\,.$ (5.85)

  6. 6.

Finally, we multiply the block encoded unitaries of $u+w$ and $v^{-1}$ to get
a $(\alpha_{6},a_{6},\epsilon_{6})$ block encoding of $y(T)$, where

$\displaystyle\alpha_{6}=\alpha_{5}\alpha_{4}$ (5.86) $\displaystyle
a_{6}=a_{5}+a_{4}$ (5.87)
$\displaystyle\epsilon_{6}=\alpha_{4}\epsilon_{5}+\alpha_{5}\epsilon_{4}\,.$
(5.88)

The gate complexity is $T_{6}=O(T_{5}+T_{4})$.

Putting it together, we get a

$(2\kappa_{V}(2\kappa_{L}s+1),\operatorname{polylog}(\kappa_{L},\kappa_{V},1/\epsilon)+\log(M+N),\epsilon)$
(5.89)

block-encoding of the solution with query complexity

$O(s\kappa_{V}\kappa_{L}\operatorname{polylog}(\frac{1}{\epsilon},s,\kappa_{V},\kappa_{L},M+N))\,.$
(5.90)

∎

###  5.3 Applications

####  5.3.1 Linear quadratic regulator

In this section, we describe how to implement the solution to the matrix
Riccati equation as a block encoded unitary. Recall from Proposition 1 that
the matrix Ricatti equation can be linearized in a larger dimension. More
specifically, we can define an $(N+M)\times(N+M)$ matrix

$A=\begin{pmatrix}F_{1}&F_{0}\\\ F_{2}&F_{3}\end{pmatrix}\,,$ (5.91)

such that if the solution of the linear ODE

$\dot{x}=Ax+b\,,$ (5.92)

is

$x(t)=\begin{pmatrix}u(t)\\\ v(t)\end{pmatrix}\,,$ (5.93)

then the solution of the matrix Riccati equation is given by

$y(t)=u(t)v^{-1}(t)\,.$ (5.94)

The initial condition for the linear ODE is

$x(0)=\begin{pmatrix}y_{0}\\\ I\end{pmatrix}\,,$ (5.95)

where $y_{0}$ is the initial condition of the Riccati equation. Note that $x$
is an $(N+M)\times M$ matrix and $y$ is an $N\times M$ matrix.

We now explain how to solve a boundary value problem (which is needed for the
LQR problem) using the above techniques. First, we need the following lemma.

######  Lemma 12.

For the linear ODE $\dot{x}=Ax+b$, if at time $t_{f}$, the solution is $x_{f}$
obtained from an initial condition $x_{0}$, then we can to $x_{0}$ from the
initial condition $x_{f}$ by time evolving the equation $\dot{x}=-Ax-b$.

###### Proof.

We have that

$x_{f}=\exp(At_{f})x_{0}+\int_{0}^{t_{f}}\exp(At)bdt\,.$ (5.96)

Using this it can be seen that

$x_{0}=\exp(-At_{f})x_{f}-\exp(-At_{f})\int_{0}^{t_{f}}\exp(At)bdt=\exp(-At_{f})x_{f}+\int_{0}^{t_{f}}\exp(-At)(-b)dt\,,$
(5.97)

where we have used the invariance of the integral under the change of variable
$t\rightarrow t_{f}-t$. ∎

Suppose we need to solve a boundary value problem given by the condition
$y_{t_{f}}=y_{f}$, then we can first create a matrix $x_{f}$ as follows.

$x_{f}=\begin{pmatrix}y_{f}\\\ I\end{pmatrix}\,.$ (5.98)

We can run the algorithm for a linear ODE by block encoding $-A$ and $-b$.
This will give us a block-encoding of $x_{0}$ by conditioning on the time
register. More specifically, we get the following block-encoded unitary

$L^{-1}z_{in}=\sum_{i=1}^{m}\ket{i,0}\bra{0,0}\otimes y_{i}\,,$ (5.99)

where $y_{i}$ is the block-encoding of the solution at different times. At
time $m=T\|A\|$, we have the block-encoding of the approximate solution
$y_{m}$. By applying a shift operator on the first register, we can bring it
into a standard block-encoding as in Definition 2. We can then use this block-
encoding of $y(0)$ to run the algorithm forwards i.e., now by block-encoding
$A$ and $b$ to the final time $T$ to get a block encoding of $y(T)$.

In the next subsection, we describe the algorithm to solve the Hamilton-
Jacobi-Bellman (HJB) equation using its connection to the matrix Riccati
equation, which produces a block encoding of a solution to the matrix Riccati
equation.

####  5.3.2 Quantum algorithm for the linear quadratic regulator problem

In this subsection, we will use the results on matrix Riccati equation from
above to give a quantum algorithm to solve the linear quadratic regulator
problem. For this problem, let us denote

$F_{0}=Q\,,F_{1}=F\,,F_{2}=GR^{-1}G^{\dagger}\,,F_{3}=F^{\dagger}\,,$ (5.100)

with their dimension being $N$ and sparsity $s$.

######  Theorem 12.

There exists a quantum algorithm that produces a quantum solution to the
linear quadratic regulator, which is a special case of the HJB equation
(2.90), with query and gate complexity

$O(s\kappa_{V}\kappa_{L}\operatorname{polylog}(\frac{1}{\epsilon},s,\kappa_{V},\kappa_{L},N))$
(5.101)

where (as before)

$\displaystyle C(A)\leq\begin{cases}C_{d}(1+C_{d}T\|F_{2}\|),&\text{if
$F_{0}=0$ and $F_{2}\neq 0$}\\\ C_{d}(1+C_{d}T\|F_{0}\|),&\text{if $F_{2}=0$
and $F_{0}\neq 0$}\\\ \exp(\mu+\|F_{0}\|+\|F_{2}\|),&\text{if both
nonzero}\end{cases}$ (5.102)
$\displaystyle\|A\|\leq\max(\|F_{0}\|+\|F_{1}\|,\|F_{2}\|+\|F_{3}\|)\,,$
(5.103)

where the log-norm of $F_{1}$ is denoted $\mu$ and

$C_{d}=\max\\{\sup_{t\in[0,T]}\|\exp(F_{1}t)\|,\sup_{t\in[0,T]}\|\exp(F_{3}t)\|\\}\,,$
(5.104)

and

$\kappa_{L}=O(C(A)T\|A\|)\,.$ (5.105)

###### Proof.

The quantum algorithm essentially constructs a block-encoding of the solution
$P$ of the matrix Riccati equation defined in (2.94) and then uses (2.95) to
construct the solution. We discuss these steps next.

First, using oracle access to the matrices $R,Q,F$ and $G$, we can construct
block-encoded unitaries corresponding to these matrices. As above, let us
denote

$F_{0}=Q\,,F_{1}=F\,,F_{2}=GR^{-1}G^{\dagger}\,,F_{3}=F^{\dagger}\,,$ (5.106)

and denote their block-encoded unitary versions $U_{0}$ through $U_{3}$. To
construct the matrix $A$ defined in (2.60), we conditionally apply each of the
unitaries $U_{i}$ i.e., first create

$V_{1}=\ket{0}\bra{0}\otimes U_{1}+\ket{1}\bra{1}\otimes U_{3}\,,$ (5.107)

then create the off-diagonals as follows.

$V_{2}=(X\otimes I)(\ket{0}\bra{0}\otimes U_{2}+\ket{1}\bra{1}\otimes
U_{0})=\ket{1}\bra{0}\otimes U_{2}+\ket{0}\bra{1}\otimes U_{0}\,.$ (5.108)

Adding these two block matrices i.e., $U=V_{1}+V_{2}$, we get a block encoded
version of the linear operator $A$. This takes $O(s)$ queries, where $s$ is
the maximum sparsity of the matrices $R,Q,F$ and $G$.

Next, we use block-encoded implementation of matrix Riccati algorithm from
Section 5.2 to create a block-encoded version of $P$ from (2.94).

Then, we create a block encoding of $R^{-1}$ and $G^{\dagger}$. Using these
and the oracle for $x(t)$, we can create the state

$\bar{u}(t)=-R^{-1}G^{\dagger}Px(t)\,.$ (5.109)

The query and gate complexities follow from the complexity of the matrix
Riccati algorithm and the complexity of matrix arithmetic. ∎

##  6 Conclusions and open questions

We have presented quantum algorithms to simulate realistic mechanical systems
such as ones that satisfy the following differential equation.

$M\ddot{q}+R\dot{q}+Vq=s\,,$ (6.1)

where $q(t)$ is the position vector of the mechanical system. This includes
damped coupled oscillators extending the work of [33]. We have shown that when
the damping strength is inverse polynomial in the number of qubits, then the
problem of estimating the kinetic energy of the system is ${\mathsf{BQP}}$
hard. This means that unless ${\mathsf{BPP}}={\mathsf{BQP}}$, no classical
algorithm can solve this problem efficiently.

We then give quantum algorithms to simulate the vector and matrix Riccati
equation. Using the algorithm for a matrix Riccati equation, one can solve the
linear quadratic regulator problem. The LQR problem is a fundamental problem
in control theory and involves a quadratic nonlinearity. Our quantum algorithm
can efficiently solve a quadratic nonlinear differential equation for a large
nonlinearity ($R\gg 1$ in the language of [34]). To our knowledge, this is the
first work that gives an efficient algorithm to solve a differential equation
with $R>1$ with rigorous performance bounds. Aside from the LQR problem, the
Riccati equation is used in predicting the stability of magnetohydrodynamics
(MHD) [41].

We list a few open problems that arise from this work. First, applying the
algorithm to estimate kinetic energy of damped coupled oscillators to estimate
other quantities of interest would be useful. Second, the quantum algorithm
for damped coupled oscillators is efficient for damping strengths that are
$\operatorname{polylog}(d)$, where $d$ is the dimension. However, the hardness
result requires that the strength be less than some inverse polynomial. Can it
be improved to a constant damping strength? A related question is if there are
classical algorithms that can estimate the kinetic energy in the presence of
large enough damping strength. Third, are there other end-to-end applications
of the quantum algorithm for the matrix Riccati equation? Another important
direction would be to generate accurate resource for end-to-end applications
like the stability of magnetohydrodynamics for instances at utility scale.

##  7 Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Fusion Energy Sciences, under Award Number DE-
SC0020264 and work supported by the Defense Advanced Research Projects Agency
(DARPA) under Contract No. HR001122C0063.

## Appendices

##  Appendix A Auxiliary lemmas

######  Lemma 13.

When $b=0$ in (3.14), the quantum algorithm from Section 4 outputs
$\hat{\dot{q}}$ such that we have

$|\dot{q}_{i}-\hat{\dot{q}}_{i}|\leq\epsilon\dot{q}_{i}\,,$ (A.1)

for some constant $\epsilon$.

###### Proof.

Let $x=(q,\dot{q})^{T}$ and let $\Pi_{i}$ be the projector onto the $i^{th}$
component of $x$. From the proof of Theorem 3 of [25], we have

$\|\Pi_{i}(x(T)-\hat{x})\|\leq\frac{me^{3}}{(k+1)!}\|\Pi_{i}x(T)\|\leq\epsilon
x_{i}\,,$ (A.2)

where we assume that we have chosen $k$ as in [25] such that

$\frac{me^{3}}{(k+1)!}\leq\epsilon\,.$ (A.3)

By choosing $i$ to correspond to $\dot{q}_{i}$, we get the statement in the
lemma. ∎

## References

  * [1] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, “Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics,” in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019, (New York, NY, USA), p. 193–204, Association for Computing Machinery, 2019.  https://doi.org/10.1145/3313276.3316366.
  * [2] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, “Grand unification of quantum algorithms,” PRX Quantum, vol. 2, p. 040203, Dec 2021.
  * [3] G. H. Low and Y. Su, “Quantum eigenvalue processing,” 2024.
  * [4] S. Lloyd, “Universal quantum simulators,” Science, vol. 273, no. 5278, pp. 1073–1078, 1996.  https://doi.org/10.1126/science.273.5278.1073.
  * [5] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, “Efficient quantum algorithms for simulating sparse Hamiltonians,” Communications in Mathematical Physics, vol. 270, p. 359–371, 2007.  https://doi.org/10.1007/s00220-006-0150-x.
  * [6] G. H. Low and I. L. Chuang, “Optimal hamiltonian simulation by quantum signal processing,” Phys. Rev. Lett., vol. 118, p. 010501, Jan 2017.  https://doi.org/10.1103/PhysRevLett.118.010501.
  * [7] G. H. Low and I. L. Chuang, “Hamiltonian Simulation by Qubitization,” Quantum, vol. 3, p. 163, July 2019.  https://doi.org/10.22331/q-2019-07-12-163.
  * [8] S. Chakraborty, A. Gilyén, and S. Jeffery, “The Power of Block-Encoded Matrix Powers: Improved Regression Techniques via Faster Hamiltonian Simulation,” in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019) (C. Baier, I. Chatzigiannakis, P. Flocchini, and S. Leonardi, eds.), vol. 132 of Leibniz International Proceedings in Informatics (LIPIcs), (Dagstuhl, Germany), pp. 33:1–33:14, Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, 2019.  https://doi.org/10.4230/LIPIcs.ICALP.2019.33.
  * [9] M. Kieferová, A. Scherer, and D. W. Berry, “Simulating the dynamics of time-dependent hamiltonians with a truncated dyson series,” Phys. Rev. A, vol. 99, p. 042314, Apr 2019.
  * [10] D. W. Berry, A. M. Childs, Y. Su, X. Wang, and N. Wiebe, “Time-dependent Hamiltonian simulation with $L^{1}$-norm scaling,” Quantum, vol. 4, p. 254, Apr. 2020.  https://doi.org/10.22331/q-2020-04-20-254.
  * [11] A. Kalev and I. Hen, “Quantum Algorithm for Simulating Hamiltonian Dynamics with an Off-diagonal Series Expansion,” Quantum, vol. 5, p. 426, Apr. 2021.
  * [12] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” Physical Review Letters, vol. 103, no. 15, p. 150502, 2009\.  https://doi.org/10.1103/PhysRevLett.103.150502.
  * [13] D. Dervovic, M. Herbster, P. Mountney, S. Severini, N. Usher, and L. Wossnig, “Quantum linear systems algorithms: a primer,” 2018.
  * [14] A. Ambainis, “Variable time amplitude amplification and quantum algorithms for linear algebra problems,” in 29th International Symposium on Theoretical Aspects of Computer Science (STACS 2012) (C. Dürr and T. Wilke, eds.), vol. 14 of Leibniz International Proceedings in Informatics (LIPIcs), (Dagstuhl, Germany), pp. 636–647, Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, 2012.  https://doi.org/10.4230/LIPIcs.STACS.2012.636.
  * [15] A. M. Childs, R. Kothari, and R. D. Somma, “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision,” SIAM Journal on Computing, vol. 46, no. 6, pp. 1920–1950, 2017.  https://doi.org/10.1137/16M1087072.
  * [16] P. C. Costa, D. An, Y. R. Sanders, Y. Su, R. Babbush, and D. W. Berry, “Optimal scaling quantum linear-systems solver via discrete adiabatic theorem,” PRX Quantum, vol. 3, p. 040303, Oct 2022.  https://doi.org/10.1103/PRXQuantum.3.040303.
  * [17] S. Lloyd, S. Garnerone, and P. Zanardi, “Quantum algorithms for topological and geometric analysis of data,” Nature Communications, vol. 7, no. 1, p. 10138, 2016.
  * [18] D. W. Berry, Y. Su, C. Gyurik, R. King, J. Basso, A. D. T. Barba, A. Rajput, N. Wiebe, V. Dunjko, and R. Babbush, “Analyzing prospects for quantum advantage in topological data analysis,” PRX Quantum, vol. 5, p. 010319, Feb 2024.
  * [19] S. Lloyd, M. Mohseni, and P. Rebentrost, “Quantum principal component analysis,” Nature Physics, vol. 10, no. 9, pp. 631–633, 2014.
  * [20] J. van Apeldoorn, A. Gilyén, S. Gribling, and R. de Wolf, “Quantum SDP-Solvers: Better upper and lower bounds,” Quantum, vol. 4, p. 230, Feb. 2020.  https://doi.org/10.22331/q-2020-02-14-230.
  * [21] C. Sünderhauf, E. Campbell, and J. Camps, “Block-encoding structured matrices for data input in quantum computing,” Quantum, vol. 8, p. 1226, Jan. 2024.
  * [22] B. D. Clader, A. M. Dalzell, N. Stamatopoulos, G. Salton, M. Berta, and W. J. Zeng, “Quantum resources required to block-encode a matrix of classical data,” IEEE Transactions on Quantum Engineering, vol. 3, pp. 1–23, 2022\.
  * [23] D. W. Berry, “High-order quantum algorithm for solving linear differential equations,” Journal of Physics A: Mathematical and Theoretical, vol. 47, no. 10, p. 105301, 2014.  https://doi.org/10.1088/1751-8113/47/10/105301.
  * [24] 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.  https://doi.org/10.1007/s00220-017-3002-y.
  * [25] H. Krovi, “Improved quantum algorithms for linear and nonlinear differential equations,” 2022.  arXiv:2202.01054.
  * [26] D. Fang, L. Lin, and Y. Tong, “Time-marching based quantum solvers for time-dependent linear differential equations,” 2022.  https://doi.org/10.48550/arXiv.2208.06941.
  * [27] D. W. Berry and P. C. S. Costa, “Quantum algorithm for time-dependent differential equations using dyson series,” 2022.  https://doi.org/10.48550/arXiv.2212.03544.
  * [28] D. An, A. M. Childs, and L. Lin, “Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters,” 2023.  https://doi.org/10.48550/arXiv.2312.03916.
  * [29] D. An, J.-P. Liu, and L. Lin, “Linear combination of hamiltonian simulation for nonunitary dynamics with optimal state preparation cost,” Phys. Rev. Lett., vol. 131, p. 150603, Oct 2023.  https://link.aps.org/doi/10.1103/PhysRevLett.131.150603.
  * [30] D. An, J.-P. Liu, D. Wang, and Q. Zhao, “A theory of quantum differential equation solvers: limitations and fast-forwarding,” 2022.  https://doi.org/10.48550/arXiv.2211.05246.
  * [31] S. Jin, N. Liu, and Y. Yu, “Quantum simulation of partial differential equations via schrodingerisation,” 2022.  https://doi.org/10.48550/arXiv.2212.13969.
  * [32] A. M. Dalzell, S. McArdle, M. Berta, P. Bienias, C.-F. Chen, A. Gilyén, C. T. Hann, M. J. Kastoryano, E. T. Khabiboulline, A. Kubica, G. Salton, S. Wang, and F. G. S. L. Brandão, “Quantum algorithms: A survey of applications and end-to-end complexities,” 2023.
  * [33] R. Babbush, D. W. Berry, R. Kothari, R. D. Somma, and N. Wiebe, “Exponential quantum speedup in simulating coupled classical oscillators,” 2023.
  * [34] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, “Efficient quantum algorithm for dissipative nonlinear differential equations,” Proceedings of the National Academy of Sciences, vol. 118, no. 35, 2021.  https://doi.org/10.1073/pnas.2026805118.
  * [35] C. Xue, Y.-C. Wu, and G.-P. Guo, “Quantum homotopy perturbation method for nonlinear dissipative ordinary differential equations,” New Journal of Physics, vol. 23, p. 123035, dec 2021.  https://doi.org/10.1088/1367-2630/ac3eff.
  * [36] S. Lloyd, G. De Palma, C. Gokler, B. Kiani, Z.-W. Liu, M. Marvian, F. Tennie, and T. Palmer, “Quantum algorithm for nonlinear differential equations,” 2020\.  https://doi.org/10.48550/arXiv.2011.06571.
  * [37] A. M. Childs, J.-P. Liu, and A. Ostrander, “High-precision quantum algorithms for partial differential equations,” Quantum, vol. 5, p. 574, Nov. 2021\.  https://doi.org/10.22331/q-2021-11-10-574.
  * [38] S. Jin and N. Liu, “Quantum algorithms for computing observables of nonlinear partial differential equations,” 2022.
  * [39] N. Linden, A. Montanaro, and C. Shao, “Quantum vs. classical algorithms for solving the heat equation,” Communications in Mathematical Physics, vol. 395, no. 2, pp. 601–641, 2022.
  * [40] D. Lewis, S. Eidenbenz, B. Nadiga, and Y. Subaşı, “Limitations for quantum algorithms to solve turbulent and chaotic systems,” 2023.
  * [41] A. S. Glasser, A. H. Glasser, R. Conlin, and E. Kolemen, “An ideal MHD W stability analysis that bypasses the Newcomb equation,” Physics of Plasmas, vol. 27, p. 022114, 02 2020.
  * [42] I. Gelfand, S. Fomin, and R. Silverman, Calculus of Variations.  Dover Books on Mathematics, Dover Publications, 2000.
  * [43] D. Liberzon, Calculus of Variations and Optimal Control Theory: A Concise Introduction.  Princeton University Press, 2012.
  * [44] M. A. De Gosson, Symplectic geometry and quantum mechanics, vol. 166.  Springer Science & Business Media, 2006.
  * [45] H. Goldstein, C. Poole, and J. Safko, Classical Mechanics.  Addison Wesley, 2002.
  * [46] M. Friswell, Dynamics of Rotating Machines.  Cambridge Aerospace Series, Cambridge University Press, 2010.
  * [47] P. Lancaster, “Stability of linear gyroscopic systems: A review,” Linear Algebra and its Applications, vol. 439, no. 3, pp. 686–706, 2013.  Special Issue in Honor of Harm Bart.
  * [48] F. H. Clarke and V. Zeidan, “Sufficiency and the jacobi condition in the calculus of variations,” Canadian Journal of Mathematics, vol. 38, no. 5, p. 1199–1209, 1986.
  * [49] J. Schiff and S. Shnider, “A natural approach to the numerical integration of riccati differential equations,” SIAM Journal on Numerical Analysis, vol. 36, no. 5, pp. 1392–1413, 1999.
  * [50] L. Hlavaty, S. Steinberg, and K. B. Wolf, “Riccati equations and lie series,” Journal of Mathematical Analysis and Applications, vol. 104, no. 1, pp. 246–263, 1984.
  * [51] W. Reid, Riccati Differential Equations.  ISSN, Elsevier Science, 1972.
  * [52] B. Anderson and J. Moore, Optimal Control: Linear Quadratic Methods.  Dover Books on Engineering, Dover Publications, 2007.
  * [53] R. Bhatia, Matrix Analysis.  Graduate Texts in Mathematics, Springer New York, 1996.  https://doi.org/10.1007/978-1-4612-0653-8.
  * [54] T. Kato, Perturbation theory for linear operators.  Grundlehren der mathematischen Wissenschaften, Springer Berlin Heidelberg, 2013.
  * [55] Y. R. Sanders, G. H. Low, A. Scherer, and D. W. Berry, “Black-box quantum state preparation without arithmetic,” Phys. Rev. Lett., vol. 122, p. 020502, Jan 2019.
  * [56] D. W. Berry, M. Kieferová, A. Scherer, Y. R. Sanders, G. H. Low, N. Wiebe, C. Gidney, and R. Babbush, “Improved techniques for preparing eigenstates of fermionic hamiltonians,” npj Quantum Information, vol. 4, no. 1, p. 22, 2018.
  * [57] Y. R. Sanders, D. W. Berry, P. C. Costa, L. W. Tessler, N. Wiebe, C. Gidney, H. Neven, and R. Babbush, “Compilation of fault-tolerant quantum heuristics for combinatorial optimization,” PRX Quantum, vol. 1, p. 020312, Nov 2020\.