text:
\section{Introduction}%
\label{sec:introduction}

In equation discovery we try to find a mathematical expression that describes
either a dataset or a model in order to get fundamental theoretical insight into
a problem that goes beyond being able to predict it.
Here we focus on finding differential equations such as
\begin{equation}
  \label{eq:pde}
  \bm u_t = F(\bm u, \bm u_x, \bm u_{xx}, \dots, \bm x, t).
\end{equation}
from data.
One way of approaching this problem is from the angle of the recently popularized
\emph{Neural Differential Equations} which involves fitting a model to represent
the RHS of Eq.~\ref{eq:pde}:
\begin{equation}
  \bm x_t = \text{NeuralDE}(\bm x, \bm\theta).
\end{equation}
\todo{explain further how to do eq discovery with NeuralDE or add citation.}
The parameters of the NeuralDE can be found by solving it with a classic DE
solver and some trial parameters $\bm\theta$ and computing the MSE between the
solution and data without needing labeled $(\bm x,\bm x_t)$ pairs.  However,
this approach comes with all the problems that typically arise during RNN
training such as bifurcations in the parameters $\bm\theta$ which complicates
convergence.

In this paper we apply a different approach which builds on
a combination of \emph{Physics Informed Neural Networks} (PINN) and
sparse identification of dynamical systems (SINDy) to discover governing
equations. This combination has been suggested by XXX et al (DeepMoD).\todo{citation}
DeepMoD comes with a couple of advantages compared to plain SINDy or NeuralDEs:
\begin{enumerate}
  \item SINDy needs access to derivatives of measurements, which is avoided by DeepMoD.
  \item DeepMoD can work with much sparser data than SINDy (again because no
    derivates have to be computed for the measurements)
  \item DeepMoD does not need to solve the DEs which avoids all the RNN training problems of NeuralODEs.
\end{enumerate}

\paragraph{Physics Informed Neural Networks}%
\label{par:physics_informed_neural_networks}
(PINNs) can solve a given partial differential equation (PDE) optionally taking into
account some measurement data. Assume we have a PDE of the form of Eq.~\ref{eq:pde}.
PINNs then fit a neural network $\hat{\bm u}$ that approximates the solution
to Eq.~\ref{eq:pde}:
\begin{equation}
  \hat{\bm u}(\bm x, t, \bm \theta) \approx \bm u(\bm x, t),
\end{equation}
where $\bm\theta$ are the parameters of the network.  To achieve this, PINNs
are based on a loss function that is composed of three parts.
\begin{equation}
  \mathcal L = \mathcal L_I + \mathcal L_B + \mathcal L_D
\end{equation}
The first part forces the solution $\hat{\bm u}$ to obey Eq.~\ref{eq:pde} on a pre-specified
spatio-temporal grid
\begin{equation}
  \label{eq:internal_loss}
  \mathcal L_I = \frac{1}{M}\sum_{i=1}^M || \hat{\bm u}_t - F(\hat{\bm u},\hat{\bm u}_x,\hat{\bm u}_{xx},\dots,\bm x, t) ||_2
\end{equation}
The second term $\mathcal L_B$ represents the boundary conditions of the PDE
and the third data term $\mathcal L_D$ is a classic MSE between of
the solution $\hat{\bm u}$ with inputs $\{\bm x_i,t_i\}_{i=1}^N$ and labels $\{\bm y_i\}_{i=1}^N$
\begin{equation}
  \mathcal L_D = \frac{1}{N}\sum_{i=1}^N || \hat{\bm u}(\bm x_i, t_i,\bm\theta) - \bm y_i ||_2
\end{equation}
Boundary and data loss can both be used but only one or the other is necessary
for the PINN to find a valid solution.

\paragraph{Sparse Identification of Non-linear Dynamics}%
\label{par:sparse_identification_of_non_linear_dynamics}
(SINDy\todo{cite}) enables PINNs to not just \emph{solve} a PDE, but to \emph{discover}
it just from measurements. By replacing the RHS of the PDE with an instantiation of
SINDy we can learn the
RHS from data alone.
The SINDy approach assumes that we are given state measurements $\bm u$ as well
as their derivatives $\bm u_t$ and performs sparse regression on a library
of basis functions. Given a (problem specific) library of basis functions.
The library has to be chosen manually, but can be very large. It typically
contains polynomials up to a degree $N$, time and/or spatial derivatives, and
other functions that might seem beneficial for the problem at hand.  By
collecting all basis functions in a function $\bm\varphi(\bm u) : N\rightarrow M$
\begin{equation}
  \bm\varphi(\bm u) = [u^1,\, u^2,\, u^1_x,\, u^1 u^2_x \dots]^T,
\end{equation}
we can represent the RHS of Eq.~\ref{eq:pde} by a linear combination of $\bm\varphi$
if the basis functions contain all necessary terms
\begin{equation}
  \bm u_t = \bm W \bm\varphi(\bm u),
\end{equation}
where we omit the implicit dependency of $\bm\varphi$ on $\bm u$ and $t$ which
comes from $\bm u(\bm x, t)$.
The result is obtained by performing sparse regression on the coefficient matrix
$\bm W$. This approach alone requires access to the derivatives of the measurements
which are often hard to obtain. There are methods to compute them in a pre-processing
step, but methods (such as the Savitzky-Golay Filter) require a high sample density.
In the next paragraph we describe how a combination of PINNs and SINDy (called
DeepMoD) can be used to circumvent this issue.

\section{Method}%
\label{sec:method}

\paragraph{DeepMoD}%
\label{par:deepmod}
ties together SINDy and PINNs by replacing the known RHS $F$ of the PDE
by a linear combination of a library of basis functions.
Essentially, it replaces the measurements $\bm u$ by the output of a universal
function approximator $\hat{\bm u}$ such that the library becomes
\begin{equation}
  \bm\varphi(\hat{\bm u}) = [\hat u^1,\, \hat u^2,\, \hat u^1_x,\, \hat u^1\hat u^2_x \dots]^T.
\end{equation}
Then we can rewrite Eq.~\ref{eq:internal_loss}:
\begin{equation}
  \label{eq:deepmod_loss}
  \mathcal L_I = \frac{1}{N}\sum_{i=1}^N || \hat{\bm u}_t(\bm x_i, t_i, \bm\theta) - \bm W \bm\varphi(\hat{\bm u}) ||_2,
\end{equation}
where $\bm W$ is a matrix of coefficients choosing from the vector of basis
functions.  Applying $\ell_1$-regularization to the basis coefficients results
in a sparse matrix which effectively discovers a minimal set of functions that
are needed to fit the given dataset.
Together with the standard PINN data loss the DeepMoD loss is
\begin{equation}
  \mathcal L = \mathcal L_I + \mathcal L_D + \lambda || \bm W ||_1,
\end{equation}
where $\lambda$ is the $\ell_1$-regularization hyper-parameter.  Above, the
boundary condition is removed because boundary conditions are generally not
known if a system is just measured and the underlying PDE is unknown.

To further generalize the approach of DeepMoD we rewrite Eq.~\ref{eq:deepmod_loss}
in terms of two functions $\text{LHS}$ and $\text{RHS}$
\begin{equation}
  \mathcal L_I = \frac{1}{N}\sum_{i=1}^N || \text{LHS}(\hat{\bm u}) - \text{RHS}(\hat{\bm u}, \bm W) ||_2,
\end{equation}
where $\text{LHS}$ represents our knowledge about the PDE and can contain any number of
terms for which we already know the coefficients and the function $\text{RHS}$ represents
the terms we want to discover. Both $\text{LHS}$ and $\text{RHS}$ of course implicitly depend
on $\bm x$, $t$, and $\bm\theta$.

Both $\text{RHS}$ and $\text{LHS}$ can contain derivatives with respect to the input $\bm x$
and/or $t$ which means that in the final optimization w.r.t. $\bm\theta$ and
$\bm W$ needs to handle second order gradient computations.
We employ forward differentiation and the package \emph{GalacticOptim.jl\todo{cite}} to
achieve this. All gradients are computed exactly via automatic differentiation.

The approach is flexible and incorporate prior knowledge about the system.
For example, given data of a system that follows the generalized Burger's equation
we can identify the form of the non-linear term only.

\newpage
\section{Experiments}%
\label{sec:experiments}

\subsection{Damped Harmonic Oscillator}%
\label{sub:damped_harmonic_oscillator}

The damped, harmonic oscillator is a linear, second order ODE
\begin{equation}
  x_{tt} = -\omega^2 x - \delta x_t
\end{equation}
with position $x$, velocity $x_t$, acceleration $x_{tt}$, frequency $\omega$
and damping coefficient $\delta$.
We obtain training data by numerically solving the true ODE and train on a grid
of 25 points ($0:0.393:3\pi$). We use only position measurements (no velocities) to learn
the generating ODE of the damped oscillator.

Typically, this second order system would now be converted to a system of first
order equations
\begin{align}
  \frac{dx_1}{dt} &= x_2 \\
  \frac{dx_2}{dt} &= -\omega^2 x_1 - \delta x_2
\end{align}
In this case, the model would have to learn by itself that the derivative of
first variable $x_1$ (i.e. the position) is simply the velocity. In order to
reliably converge to this solution the model would need velocity data.
We can avoid this additional data collection by directly learning the second
order ODE. To achieve this we define the $\text{LHS}$ as
\begin{equation}
  \text{LHS}(\hat x) = \hat x_{tt}
\end{equation}
and provide the following basis functions for the $\text{RHS}$
\begin{equation}
  \bm\varphi(\hat x) = [\hat x,\,\hat x^2,\,\hat x^3,\,\hat x_t]^T.
\end{equation}
We therefore want to solve the problem of finding a network $\hat x : \mathbb
R\rightarrow\mathbb R$ and sparse coefficients $\bm w$
\begin{equation}
  \hat x_{tt} \approx \bm w \cdot \bm\varphi(\hat x),
\end{equation}
where $\bm w$ is just a vector in this case because we have only one state variable $x$.

\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{\input{damped-harmonic.tikz}}
  \caption{Damped harmonic oscillator discovery. $\bm w = [-2.20,\, 0.08,\, 10^{-5},\,
  -0.51]$. Target parameters $\omega^2=2.25$ and $\delta=0.5$.}%
  \label{fig:damped-harmonic}
\end{figure}

\subsection{Lotka-Volterra}%
\label{sub:lotka_volterra}

The Lotka-Volterra (LV) equations are a pair of non-linear, ordinary differential equations
\begin{align}
  \label{eq:lotka-volterra}
  x_t &= \alpha x - \beta xy\\
  y_t &= \delta xy - \gamma y
\end{align}
that describe a predator-prey system where $x$ is the number of prey and $y$
the number of predators. The parameters $\alpha$, $\beta$, $\gamma$, and $\delta$
define the interaction of the two species.

We obtain training data by solving the true LV equations with
\emph{DifferentialEquations.jl} and define a grid of 100 points (-8:0.16:8)
from which we sample batches of 25 points in each training step.

The $\text{LHS}$ is defined as $\hat{\bm u}_t = [x_t, y_t]^T$ and the $\text{RHS}$ contains
the basis functions
\begin{equation}
  \bm\varphi(\hat{\bm u}) = [x,y,x^2,y^2,xy,x^3,y^3]^T.
\end{equation}
We use a separate network for each one of the two state variables $x$ and $y$
to enforce the two variables to be independent of each other. This way the only
way of the two variables interacting is via the terms in the $\text{RHS}$ and not
within the network $\hat{\bm u}$ itself. This leads to more reliable convergence to
the basis function coefficients we are seeking.

\begin{figure}
  \centering
  \input{lotka-volterra.tikz}
  \caption{Data and solution of the LV equations as obtained by the PINN that is
  part of DeepMoD in the left column, coefficients of the used basis functions
  during training in the right column.}%
  \label{fig:lotka-volterra}
\end{figure}\todo{x-axis does not show training steps but time...}

\begin{table}
  \centering
  \caption{(Hyper-)parameters to obtain the results shown in Fig.~\ref{fig:lotka-volterra}.}
  \label{tab:lotka-volterra}
  \begin{tabular}{c}
    please fill out...
  \end{tabular}
\end{table}

\subsection{Duffing Oscialltor}%
\label{sub:duffing_oscialltor}

The Duffing Oscillator is a non-linear, second order, differential equation which
includes both damping and a forcing term
\begin{equation}
  x_{tt} + \delta x_t + \alpha x + \beta x^3 = \gamma \cos(\omega t),
\end{equation}
where $\delta$ controls the damping, $\alpha$ controls the linear
stiffness, $\beta$ controls the non-linearity in the restoring force, $\gamma$
is the amplitude of the forcing and $\omega$ is its frequency.
The Duffing equation exhibits multiple period cycles and even chaotic behaviour
for different coefficients. We fix $\alpha=1$, $\beta=-1$, $\delta=0.3$, and $\omega=1.2$.
The forcing coefficient $\gamma$ is varied to obtain different period cycles and
chaotic solution as shown in Fig.~\ref{fig:duffing-period-3} and Fig.~\ref{fig:duffing-chaos}.

For the discovery of this equation we assume that the frequency $\omega$ of the
forcing is known so we aim to find the coeffcients $\alpha$, $\beta$, $\gamma$,
and $\delta$ from the basis
\begin{equation}
  \bm\varphi(\hat x) = [\hat x,\, \hat x^2,\, \hat x^3,\, \hat x^4,\, \hat x_t,\, \cos(\omega t)]^T.
\end{equation}

\todo{}

\begin{figure}
  \centering
  \input{duffing-period-3.tikz}
  \caption{Duffing oscillator with forcing coefficient $\gamma=0.29$ resulting
  in a period-3 cycle. Discovered $\bm w = [1.001,\, 10^{-5},\, -1.010,\, 0.008,\,  -0.299,\,  0.290]^T$.
  Target parameters $[1,\, 0,\, -1,\, 0,\, -0.3,\, 0.29]^T$.}%
  \label{fig:duffing-period-3}
\end{figure}

\begin{figure}
  \centering
  \input{duffing-chaos.tikz}
  \caption{Duffing oscillator with forcing coefficient $\gamma=0.5$ resulting
  in chaotic behaviour. Discovered $\bm w = [0.964,\, -0.019,\, -0.974,\,  0.007,\,  -0.298,\,  0.484,\,]^T$.
  Target parameters $[1,\, 0,\, -1,\, 0,\, -0.3,\, 0.5]^T$.}%
  \label{fig:duffing-chaos}
\end{figure}

\newpage
\subsection{Burgers' Equation}%
\label{sub:burgers_equation}

Burgers' Equation is a simplification of the Navier-Stokes equation. It can be
investigated in one spatial and one temporal dimension and hence makes for a good
example of a low dimensional PDE
\begin{equation}
  \label{eq:burgers}
  u_t = -uu_x + \nu u_{xx}
\end{equation}
with the speed of the fluid $u$, spatial coordinate $x$, temporal coordinate $t$,
and viscosity $\nu$.
The solution to Eq.~\ref{eq:burgers} with sinusoidal initial conditions is shown
in Fig.~\ref{fig:burgers-sin}.
\begin{figure}
  \centering
  \includegraphics[width=0.4\linewidth]{burgers-sin.pdf}
  \caption{Solution to Burgers' equation with initial condition $u(x,0)=\sin(2\pi x)$ and $\nu=0.005$.
    Note the shock-wave that is forming around $x=0.5$.}%
  \label{fig:burgers-sin}
\end{figure}

\begin{figure}
  \centering
  \input{burgers.tikz}
  \caption{Discovery of the two relevant terms in Burgers' equation.
  Discovered coefficents $\bm w = [-0.990,\,  0.005,\,  <10^{-5},\,  <10^{-5},\, <10^{-5},\, <10^{-5},\, -0.003,\, <10^{-5}]^T$. Target coefficients $[1, 0.005, 0, 0,...]$}%
  \label{fig:burgers-discovery}
\end{figure}

The generalized Burgers equation extends the linear term in front of $u_x$ to
an arbitrary function $c(u)$:
\begin{equation}
  \label{eq:general-burgers}
  u_t - \nu u_{xx} = -c(u)u_x,
\end{equation}
which we will use to demonstrate how to encode addition prior knowledge into
DeepMoD that goes beyond a single-term LHS.
Eq.~\ref{eq:general-burgers} is already rearranged such that we only have to
define the LHS of DeepMoD as
\begin{equation}
  \text{LHS}(\hat u) = \hat u_t - \nu\hat u_{xx}.
\end{equation}
In this case we choose the basis functions as polynomials up to degree three
elementwise multiplied by the spatial derivative
\begin{equation}
  \bm\varphi(\hat u) = [1,\, \hat u,\, \hat u^2,\, \hat u^3]^T \odot \hat u_x
\end{equation}
to discover the function $c(u)$ from an expression $c(u) = (w_1 + w_2u + w_3u^2 + w_4u^3)u_x$.

\begin{figure}
  \centering
  \input{burgers-general.tikz}
  \caption{Discovery $c(u)$ in Burgers' equation.}%
  \label{fig:burgers-general-discovery}
\end{figure}

\end{document}

\section{Equation discovery with Physics Informed Neural Networks}%

Assuming a PDE of the form
\begin{equation}
  \frac{d\bm u}{dt} = F(\bm u, \partial_x u, \partial^2_x u, \,\dots\,, \bm x, t),
\end{equation}
with a vector of dependent variables $\bm u = \{u_m(x_n,t)\}_{mn}$, partial
derivatives $\partial_x u=\frac{d\bm u}{dx}$, we want to discover the form $F$ of the
RHS.
\emph{Physics Informed Neural Networks} (PINNs) can be used to \emph{solve}
(partial) differential equations \emph{given} the form $F$ of the PDE. They are
trying to find a neural network $\hat u$ that approximates the solution to the PDE
\begin{equation}
  \hat u(\bm x, t, \bm\theta) \approx \bm u(\bm x, t)
\end{equation}
by minimizing a loss function that is made of three different components: An
\emph{internal} loss, a \emph{boundary} loss, and an optional \emph{data} loss.
The internal loss $\mathcal L_F$ enforces solutions that follow the form of the PDE $F$
\begin{equation}
  \label{eq:pinn_f}
  \mathcal L_F(\bm\theta) = \sum_{ij} \norm{\frac{d\hat u(\bm x_j,t_i,\bm\theta)}{dt}
                           - F(\hat u,\partial_x\uh, \,\dots\,,\bm x_j, t_i)}_2
\end{equation}
which is evaluated on a spatio-temporal grid (indicated by the indices $j$ and $i$).
The boundary loss $\mathcal L_B$ is a sum of a collection of boundary conditions
(e.g. $\hat u(x,0,\bm\theta)-\sin(x)\approx 0$).
The data loss $\mathcal L_D$ can be optionally added to the PINN loss:
\begin{equation}
  \mathcal L_D(\bm\theta) = \sum_{ij} \norm{\hat u(\bm x_j,t_i,\bm\theta) - \bm u_{ij}}_2,
\end{equation}
which results in the final loss $\mathcal L$
\begin{equation}
  \mathcal L(\bm\theta) =  \lambda_1\mathcal L_F(\bm\theta)
                         + \lambda_2\mathcal L_B(\bm\theta)
                         + \lambda_3\mathcal L_D(\bm\theta)
\end{equation}

In order to \emph{discover} an unknown $F$ from a set of measurements we can
construct an $\hat F$ from a set of library functions $\varphi_i$
\begin{equation}
  \hat F(\hat u,\partial_x\uh, \,\dots\,, \bm x, t) = \bm W \bm\varphi(\hat u,\bm x, t)
\end{equation}
where $\bm W$ is a matrix of coefficients for the vector library of library
functions $\bm\varphi$.
Eq.~\ref{eq:pinn_f} therefore becomes
\begin{equation}
  \mathcal L_F(\bm\theta) = \sum_{ij} \norm{\frac{d\hat u(\bm x_j,t_i,\bm\theta)}{dt}
                                          - \bm W\bm\varphi(\hat u,\bm x_j, t_i)}_2.
\end{equation}
In order to obtain only relevant terms from the library $\bm\varphi$ we have to
introduce an additional regularization term to the loss
\begin{equation}
  \mathcal L(\bm\theta) =  \lambda_1\mathcal L_F(\bm\theta)
                         + \lambda_2\mathcal L_B(\bm\theta)
                         + \lambda_3\mathcal L_D(\bm\theta)
                         + \lambda_4\norm{\bm W}_1
\end{equation}

\subsection{Discovering the Lotka-Volterra System}%
\label{sub:discovering_the_lotka_volterra_system}

\begin{figure}
  \centering
  \includegraphics[width=0.6\linewidth]{figures/lotka-volterra.pdf}
  \caption{Lotka-Volterra discovery. Basis functions: $x, y, x^2, y^2, xy$}%
  \label{fig:lotka-volterra}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.6\linewidth]{figures/lotka-volterra-badrun.pdf}
  \caption{Bad lotka volterra run with additional basis function:  Basis
  functions: $1, x, y, x^2, y^2, xy$}%
  \label{fig:figures/lotka-volterra-badrun}
\end{figure}

\subsection{Discovering the a Damped Harmonic Oscillator}%
\label{sub:discovering_the_lotka_volterra_system}

\begin{figure}
  \centering
  \includegraphics[width=0.6\linewidth]{figures/damped-harmonic.pdf}
  \caption{Damped harmonic oscillator discovery}%
  \label{fig:damped-harmonic}
\end{figure}

\section{Solving Navier Stokes with PINNs}%
\label{sec:navier_stokes}

\begin{itemize}
  \item NSFnets~\cite{Jin2020} solve the Navier-Stokes equations either in the
    velocity-pressure formulation or the vorticity-velocity formulation with
    networks that learn the mapping $(t,x,y,z)\rightarrow(u,v,w,p)$
    (or $(t,x,y,z)\rightarrow(u,v,w,\xi)$). They employ \emph{dynamic weights}
    to determine good coefficients for the different loss contributions.

  \item Conservative PINNs (cPINN) comes with additional \emph{interface
    conditions} and solve a viscous Burgers equation (a form of Navier-Stokes)
    \cite{Jagtap2020}.  Interface conditions are additional loss terms that
    enforce e.g. flux continuity conditions. Additionally they divide the domain
    in subdomains which each get a separate network (e.g. more complex networks
    for harder subdomains). Networks are joined by interface conditions

  \item SINDy: Discovered \emph{mean-field model} of vortex shedding. This was
    done by reducing the dimension of flow simulation data to three coordinates
    which were then discovered by SINDy \cite{Rudy2017}. As far as I understand
    they were not learning from the full simulation data but from the reduced
    system.
\end{itemize}

\section{Literature Review}

Our problem can be split into two classes of problems. The first one,
which is often called \eemph{system identification} is concerned with estimating
parameters of a model for which we know the functional form $F$. Consider, for
example, trying to identify the frequency of a harmonic oscillator
\begin{equation}
  \label{eq:osc}
  \ddot{x} = -\omega^2 x = F(x,\dot x,\omega,t)
\end{equation}
from measurements.
The second class of problems which we will refer to as \eemph{equation
discovery} assumes only partial/no knowledge of the functional form of the
problem (making it much harder).
When it comes to doing sparse regression the borders between the two approaches
become blurry, because the choice of the model (or the basis functions) already
restricts the space of $F$. As an example, we can prescribe $F$ to be a linear
combination. Pruning as many parameters as possible could be viewed as discovering
a linear equation but ultimately it then is just a fit of the restricted $F$ to
the date (i.e. system identification).\\

In our case we also have to distinguish between finding an expression for $F$
(\eemph{discovering the ODE} Sec.~\ref{sub:neural_odes}) or finding the
solution $x$ (\eemph{discovering the solution}
Sec.~\ref{sub:physics_informed_neural_networks}).\\

Equation discovery can be split into two different approaches. The first
one being \eemph{genetic programming}, which tries to fit symbolic expressions to given
datasets as described e.g. by \cite{Schmidt2009} and implemented
in the Eureqa package \citep{Praksova2011}.
Sections~\ref{sec:discovering_symbolic_models} and \ref{sec:ai_feynman} are
extensions to this approach.

The second class of approaches focusses on \eemph{sparse regression}. For the purpose of
equation discovery they all essentially want to solve $\ell_0$-regularization
which is a combinatorial, not differentiable problem. Hence most approaches
relax the problem and go with $\ell_1$-regularization instead.  PINNs
(Sec.~\ref{sub:physics_informed_neural_networks}) and Neural ODEs
(Sec.~\ref{sub:neural_odes}) can be used here and SINDy (Sec.~\ref{sec:sindy})
and EQL (Sec.~\ref{sec:equation_learner}) are probably the most notable
instantiations of this approach.

\subsection{Differential Equations}%
\label{sec:differential_equations}

\emph{Not directly related to equation discovery}.
Quick recap on two fundamental scientific machine learning methods.  A good
place to start reading about them is
\href{https://mitmath.github.io/18337/lecture3/sciml.html}{Rackauckas' lecture}
on scientific
computing.

\subsubsection{Physics Informed Neural Networks (PINNs)}%
\label{sub:physics_informed_neural_networks}

Given a system of ordinary differential equations
\begin{equation}
  \dot{\bm x} = F(\bm x, t)
\end{equation}
with $t\in[0,1]$ and known initial condition $\bm x(0)=\bm x_0$ we can solve
the ODE by fitting a neural network (NN):
\begin{equation}
  \text{NN}(\bm\theta,t) \approx \bm x(t)
\end{equation}
by training with a loss that just enforces the ODE and its boundary conditions
\begin{equation}
  \text{MSE}_F = \sum_i \left(\frac{\text{NN}(\bm\theta,t_i)}{dt}
                        - F(\text{NN}(\bm\theta,t_i),t_i) \right)^2
                        + (NN(0) - \bm x_0)^2.
\end{equation}
PINNs~\citep{Raissi2019} can incorporate measurements
with an additional term in the loss function
\begin{align}
  \text{MSE}_x &= \sum_i||\text{NN}(\bm\theta, t_i) - \bm x_i||_2\\
  \mathcal L(\bm\theta) &= \text{MSE}_x + \text{MSE}_h
\end{align}
The information encoded in the loss also does not have to be the full ODE but can
be just an additional constraint that the solution should take into account.
In the simplest case: $h(t)=\dot{\bm x}-F(\bm x,t)$.
Hence, a PINN no longer finds just \emph{any} solution that fits the data,
but one that agrees with data \emph{and} the given physical constraint.
In order to use PINNs for system identification the network can be limited in
its size and/or regularizations can be added to $\mathcal L$.
By replacing the NN in the above equation with a Chebyshev basis one essentially
ends up at the approach of solving DEs on DQCs by~\cite{Kyriienko2020}.

\subsubsection{Neural ODEs}%
\label{sub:neural_odes}

Neural ODEs (as opposed to PINNs) do not fit a solution of the ODE but the RHS
of the ODE itself:
\begin{equation}
  \dot{\bm x} = f(\bm x, \bm\theta, t)
\end{equation}
This is achieved by solving the ODE with the current parameters $\bm\theta$
and computing the distance to the data
\begin{equation}
  \mathcal L(\bm\theta) = ||\text{solve}(\text{NN}, \bm\theta, \bm t) - \bm X||_2,
\end{equation}
\todo{fix notation neural ode loss}
where $\bm t$ are the times at which $\bm X$ is sampled and $\text{solve}$
refers to some method that can solve a given differential equation (this is
also called a shooting method). This training setup forces the NN to have very
little drift in early time-steps as opposed to training simply with examples
of $\bm X$ and $\dot{\bm X}$ but is much trickier to train due to the classic
RNN training problems (i.e. vanishing/exploding gradients and bifurcations).

\eemph{Universal differentitial equations}
(UDE)~\citep{Rackauckas2020} extend Neural ODEs by hard-coding
certain parts of the ODE model which are known.
\begin{equation}
  \dot{\bm x} = f(\bm x, t, \text{NN}(\bm\theta,t))
\end{equation}
The known parts of $f$ do not have to be learned and the model can focus on the
parts that are missing.

\subsection{Sparse Identification of Nonlinear Dynamics (SINDy)}%
\label{sec:sindy}

SINDy \citep{Brunton2016, Rudy2017, Champion2019} tries to find the RHS of
\begin{equation}
  \dot{\bm x} = F(\bm x, t)
\end{equation}
that fits snapshots of data $\bm X = [\bm x_1 \, \bm x_2 \, \dots \, \bm
x_m]^T$ and $\dot{\bm X} = [\dot{\bm x}_1 \, \dot{\bm x}_2 \, \dots \,
\dots{\bm x}_m]^T$ by providing a library of basis functions $\bm\Theta(\bm X)
= [\bm\theta_1(\bm X)\, \dots \, \bm\Theta(\bm X)]^T$ on which sparse
regression is applied:
\begin{equation}
  \mathcal L(\bm\theta) = ||\dot{\bm X} - \bm\Theta(\bm X)\bm\Xi||_2 + \lambda||\bm\Xi||_1,
\end{equation}
where $\bm\Xi$ can pick out any combination of inputs $\bm x_i$ and basis
functions $\bm\theta_i$. By using an auto-encoder to first compress the given
dataset SINDy can by used for high-dimensional datasets as well.
SINDy needs access to both measurements $\bm X$ as well as their derivatives
$\dot{\bm X}$. The derivatives can be computed from $\bm X$ e.g. via the
Savitzky-Golay filter (which needs quite dense measurements, \cite{Schafer}).
Instead of using a filter to estimate the derivative, SINDy can be used in
conjunction with UDEs to estimate $\dot{\bm X}$ and discover unknown parts of
the given dataset (briefly described in~\cite{Rackauckas2020}).

\subsection{Koopman Operator Theory}%
\label{sub:koopman_operator_theory}

\emph{Not directly related to equation discovery}. Koopman Operator Theory
(extensively described in the notes of~\cite{Brunton2019Notes}) provides a way
of linearising highly non-linear dynamical systems of the form
\begin{equation}
  \bm x_{t+1} = \bm F(\bm x_t)
\end{equation}
by finding an \eemph{infinite-dimensional} but \eemph{linear} operator $\mathcal K$.
Consider a measurement function $g : \bm M \rightarrow \mathbb R$ which represents
an observation of the state $\bm x_t$. The Koopman operator $\mathcal K$ defines
a linear, infinite-dimensional dynamical system which advances an observation
to the next time step:
\begin{equation}
  g(\bm x_{t+1}) = \mathcal K g(\bm x_t)
\end{equation}
Being infinite-dimensional, the Koopman operator is somewhat problematic to use
in practice, which is why it is typically approximated by a finite set of its
eigenfunctions $\varphi$ (with corresponding eigenvalue $\lambda$):
\begin{equation}
  \varphi(\bm x_{t+1}) = \mathcal K \varphi(\bm x_t) = \lambda \varphi(\bm x_t)
\end{equation}
With these eigenfunctions we can approximate any measurement function
\begin{equation}
  g(\bm x_t) = \sum_j^\infty v_j\varphi_j(\bm x_t)
\end{equation}
\todo{$v_j$ are Koopman modes. what does that mean?}

To find the Koopman eigenfunctions we can take the time derivative of the continuous
time formulation of the Koopman operator:
\begin{equation}
  \frac{d}{dt}\varphi(\bm x) = \mathcal K\varphi(\bm x) = \lambda \varphi(\bm x)
\end{equation}
which results in a partial differential equation for the eigenfunction
\begin{equation}
  \nabla \varphi(\bm x) \cdot \bm f(\bm x) = \lambda \varphi(\bm x).
\end{equation}
Solving this differential equation for a number of important eigenfunctions
amounts to finding an approximation of the Koopman operator.

Another way of finding eigenfunctions $\bm\varphi(\bm x)$ is by leveraging
auto-encoders. The each variable in latent space of the auto-encoder can be
trained to represent an eigenfunction by additional terms in the loss function
that enforce the linearity and the prediction constraint:
\begin{equation}
  \mathcal L = ||\bm x - \bm\varphi^{-1}(\bm\varphi(\bm x))||
             + \alpha||\bm\varphi(\bm x_{t+1}) - \bm K\bm\varphi(\bm x_t)||
             + \beta||\bm x_{t+1} - \bm\varphi^{-1}(\bm K\bm\varphi(\bm x_t))||.
\end{equation}
The correct loss function used in~\cite{Lusch2017} contains some additional
terms like a $\ell_2$-regularization.

\subsection{Equation Learner}%
\label{sec:equation_learner}

The Equation Learner (EQL, \cite{Sahoo2018}) could be regarded as a
layered SINDy approach. The model consists of linear transformations $\bm y=\bm
W\bm x+\bm b$ followed by a layer of \emph{inhomogeneous} functions.  These
functions include unary functions ($\sin$, $\cos$, etc.) and binary functions
like multiplication and division.  The divergence problems with division are
approached by only including division in the last layer and by using curriculum
learning to slowly progress from a strongly regularized division function
\begin{equation}
  h^\theta(a,b) = \begin{cases}
    \frac{a}{b} \qquad \text{if}\quad b>\theta \\
    0 \qquad \text{otherwise}
  \end{cases}
\end{equation}
to the unregularized one. Further, the assumption is that nature does not
diverge, so one branch of division ($b>0$) should suffice as a basis function
which is enforced by a penalty term $P^\theta$. The network is trained in three
phases: 1) train without regularization 2) train with $\ell_1$-regularization
3) prune weights below a threshold $w<0.001$ and train without regularization
(which they call fixed $\ell_0$).

\subsection{Sparse Symplectically Integrated Neural Networks (SSINNs)}%
\label{sec:sparse_symplectically_integrated_neural_networks_ssinns}

Inspired by \emph{Hamiltonian Neural Networks} (HNNs; \cite{Greydanus2019})
SSINNs~\citep{Dipietro2020} use two networks $V_{\theta_1}$ (potential energy)
and $T_{\theta_2}$ (kinetic energy) and symplectically integrate them to obtain
surrogates for separable Hamiltonians. They use $\ell_1$-regularization to get
sparse models. The two networks are restricted to a single linear layer in the
paper, but they show how the approach could be extended to more complicated
models.

\subsection{Discovering Symbolic Models}%
\label{sec:discovering_symbolic_models}

Genetic programming has been used to find symbolic expressions for
physical laws from experimental data~\citep{Schmidt2009, Zaremba2014}.
The approach unfortunately does not scale well,
which is why \citep{Cranmer2020} introduced a hybrid approach of
genetic programming and deep learning. They fit a \eemph{graph neural network}
GNN to a dataset (e.g. a system of particles) and then try to find symbolic
expressions for smaller subparts of the network.

\subsection{AI Feynman}%
\label{sec:ai_feynman}

This approach \citep{Udrescu2020} is a combination of genetic programming
and sparse regression (which the authors call symbolic regression) which also
build on recursively breaking down the problem into simpler and simpler subtasks
before fitting symbolic expressions to them.

\subsection{Disentangling Autoencoders}%
\label{sec:disentangling_autoencoders}

This is not directly equation discovery related, but \cite{Iten2020} use
disentangling autoencoders to learn compressed representations of time-series.
They prescribe the size of the latent vector that they are using to encode
their data and then analyze the latent space. For their examples they show that
the latent space represents i.e. conserved quantities (momentum in particle
collisions), or other relevant quantities like frequency and damping factor of
a pendulum. It seems like they are not doing much more than prescribing the
latent size and enforcing the latent variables to independent though.

\subsection{$\ell_0$-regularization}%
\label{sec:l0_regularization}

Most (if not all) the sparse regression problems only use $\ell_1$-regularization
because it can be used as a differentiable surrogate for $\ell_0$.
Hence, solving $\ell_0$-regularization on a QC would be of interest to many fields
far beyond equation discovery.

Given a linear regression problem with parameters $\bm\theta$ and assuming we are
interested in pushing as many parameters in that model to zero (i.e. we want
sparsity). The state of the art for doing this is $\ell_1$-regularization:
\begin{equation}
  \mathcal L(\bm\theta) = ||\bm X\bm\theta - \bm y||_2 + \lambda||\bm\theta||_1
\end{equation}
This is computationally efficient and works well in many cases, but $\ell_1$
tends to \emph{over-shrink} the parameters (i.e. bias them too much towards zero).
For equation discovery (and many other problems) we really want to avoid this
shrinkage. The $\ell_0$-regularization is the much more natural objective, because
it encourages parameters to be either exactly zero or some other value:
\begin{equation}
  \mathcal L(\bm\theta) = ||\bm X\bm\theta - \bm y||_2 + \lambda||\bm\theta||_0
    = ||\bm X(\bm\theta\odot\bm z) - \bm y||_2 + \lambda\sum_i z_i,
\end{equation}
where $\bm z \in \{0,1\}^n$.
Unfortunately this objective is not differentiable. It essentially poses a
combinatorial problem combined with linear regression \citep{Chen2020} for
which some relaxations exist like \citep{Louizos2017}. A comparison of
different approaches to $\ell_0$-regularization
including their computational aspects is given in \cite{Hazimeh2020} (easier to
read \href{https://www.mit.edu/~rahulmaz/L0Reg.pdf}{version}).

Formally the objective is
\begin{equation}
  \arg\max_{\bm\theta\in\mathbb R^n} ||\bm X\bm\theta - \bm y||_2
  \quad \text{s.t.} \quad
  ||\bm\theta||_0 \leq k
\end{equation}
where $k$ controls the model size.
The most established way of doing this seems to be via the \texttt{leaps}
function of \texttt{R} (feasible for less than 30 parameters in $\bm\theta$).
Currently \eemph{mixed-integer optimization} (MIO) by \citep{Bertsimas2016}
seems to be fastest at solving this problem (for up to 1000 parameters).
Intro to $\ell_0$-regularization with
\href{https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/linear_regression/l0_regression_gcl.ipynb}{Gurobi}
which links to \href{https://arxiv.org/abs/1507.03133}{this paper by Bertsimas}.

The \eemph{Quantum Approximate Optimization Algorithm}
(\href{https://qiskit.org/textbook/ch-applications/qaoa.html}{QAOA}) can solve
the problem
\begin{equation}
  \arg\max_{x\in\{0,1\}^n} \mathcal L(x)
\end{equation}

Can we split our problem into a binary part (solved on DQC) and the regression
part (solved on CC) as suggested here:
\href{https://arxiv.org/pdf/2010.07852.pdf}{MIO on QC}?
\href{https://arxiv.org/pdf/1508.04212.pdf}{Constrained Quantum Annealing} might
be relevant as well.

It could be possible to get at this problem with a more \emph{quantum} approach
than what is currently done on classical computers - with ancilliary bits that
exploit superpositions of $\ell_0$-regularized terms\todo{Niklas: sorry for the
fuzzy formulation. I still have to understand how this would be done exactly}.

\subsection{Maybe}%
\label{sec:maybe}

\begin{itemize}
  \item \href{https://github.com/AllanSCosta/occam-net}{OccamNet}
  \item Dynamic Mode Decomposition \& Koopman Operator
  \item \href{https://arxiv.org/pdf/1812.04426.pdf}{PDE-Net} /
    \href{https://arxiv.org/pdf/2006.02377.pdf}{RODE-Net}
  \item \href{https://export.arxiv.org/pdf/1907.00821}{Process Based Modelling}
\end{itemize}

\printbibliography

\end{document}

\section{Introduction}%
\label{sec:introduction}

Recent advances in automatic differentiation in classical computers have led to novel machine learning based approaches to scientific computing, often referred to as scientific machine learning.
For example, incorporating neural networks into classical ordinary differential equation solvers \cite{Chen2018a}, or manipulating neural network derivatives towards prescribed differential equations \cite{Raissi2019}.
The applications of such techniques have included inference of unobserved states \cite{Mao2020}, inference of differential equation parameters \cite{Tartakovsky2018} and approaches to solve differential equations \cite{Raissi2019}, or discover differential equations from observated data \cite{Brunton2016a}.
Such methods depend on machine learning techniques and are often end-to-end differentiable.

Recent progress in automatic differentiation on quantum computers \cite{Mitarai2018, vidal2018calculus, Crooks2019a, killoran2019, Mari2021, Banchi2021measuringanalytic, Cerezo2021, hubregtsen2021singlecomponent, hubregtsen2021singlecomponent, kyriienko2021} introduce the prospect of extending recent classical results in scientific machine learning to a quantum setting.
These techniques have already led to quantum methods of solving differential equations \cite{Kyriienko2020} which are analogous to methods involving classical neural networks \cite{Raissi2019}.
However, the existence of quantum automatic differentiation opens the door to further cross pollination between classical scientific machine learning and quantum computing.

In this work, we target an important task in scientific machine learning and modelling in general; while many systems can be modelled using differential equations, and their solving constitutes an important task in phenomenological understanding, in many practical settings these equations are not (fully) known. In such cases, one may have access to observations on some target system, and some initial idea of the dynamics.

Parameter inference in differential equations is an essential part of combining theoretical models with empirical observations.
Oftentimes a given mechanism might be speculated for an observed process in the form of a differential equation, but whose coefficients must be inferred from empirical observations.

Model discovery (often referred to as equation discovery) tries to find a mathematical expression
that describes a given dataset. We specifically tackle the problem of
discovering a differential equation only from observations of a given system.
By \emph{discovery} we mean an interpretable (i.e. human readable) model in
symbolic form, which can be converted directly into a differential equation.

This task is generally not suited for treatment by finite-differencing techniques, since such applications require many repeated simulations. Instead, to treat such a case of mixture of equation solving and model discovery efficiently, we generalize two techniques from the realm of classical scientific machine learning to a quantum setting: differential equation parameter inference and differential equation model discovery.
We also highlight conceptual similarities between scientific machine learning and quantum computing that might point to further research avenues.

\begin{figure}[ht]
  \centering
  \resizebox{\textwidth}{!}{\input{qmd.tex}}
  \caption{Schematic of the quantum model discovery approach (QMoD\xspace). The solution to the differential equation is represented by a (quantum)
    universal function approximator (UFA) with variational angles $\bm\theta$.
    Both the forward pass and the gradients w.r.t. parameters $\bm\theta$ \emph{as well as}
    gradients w.r.t. inputs (in this example $x$ and $t$) are computed on the DQC.
    QMoD\xspace learns a sparse $\bm w$ which chooses only a few functions from the library of basis functions $\bm\varphi$.
    Gradients w.r.t. the basis function coefficients $\bm w$ are computed classically.
    The figure depicts how QMoD\xspace learns the solution $\hat u$ and a sparse $\bm w$ that only
    contains the two basis functions that are necessary to represent the input data from Burgers' equation.
    }%
  \label{fig:qmd_schematic}
\end{figure}

We first recap how quantum computers can be used to represent solutions to differential equations via variational procedures.
Then we introduce the broad field of scientific machine learning and propose how quantum computers can be used for these nascent developments.

\paragraph{Universal Function Approximation on Quantum Devices.}%
\label{par:universal_function_approximation_on_quantum_devices}

For the past several decades, the number of transistors on a similarly-sized integrated chip doubles roughly every two years of so; this is known as Moore's law \cite{Moore1965} and is an observation rather than a guarantee. In the past 10 years or so, a shift has started where more computational power can only be reached by parallelization rather than minimization or increasing clock speeds \cite{theis2017}; this is due to the fundamental limitations of electronics. While alternative classical hardware proposals are interesting, such as photonic \cite{krishnamoorthy2009} or biological \cite{Bonnet2013} computers, these will quickly hit a wall of unsustainable exponential speed increase too. \\

A completely different compute paradigm, quantum computing, has existed for over 4 decades but has only recently emerged as a truly promising next stage of computational hardware developments. Large corporate players including IBM \cite{Kandala2017b} and Google \cite{Arute2019a}, as well as a range of startups including IonQ, Pasqal, Xanadu and Rigetti, are investing large amounts of funding for such hardware, and with that the development of quantum computational algorithms has gotten a resurgence from the initial fundamental algorithms like Shor's \cite{Shor1995} or Grover's \cite{Grover1996} algorithms. Some of the promising areas where universal fault-tolerant quantum computing hardware could show a possible quantum advantage include quantum chemistry simulations \cite{Cao2018e, elfving2020quantum} and quantum-enhanced machine learning \cite{Biamonte}.\\

Recent algorithm developments have also focused on compatibility with current noisy intermediate scale quantum (NISQ) devices \cite{Preskill2018c}. Such algorithms are typically variational in nature, and are being applied to fields such as quantum chemistry \cite{Kandala2017b, Elfving2020b} and heuristic optimization \cite{Farhi2014, Borle2020}. \\

Furthermore, such variational quantum algorithms also offer a promising platform for the field of Machine Learning, which largely depends on training models to fit datasets and making predictions. Quantum Machine Learning (QML) algorithms have been developed for solving classification \cite{Havlicek2019SupervisedSpaces} and regression \cite{Mitarai2018} tasks; the techniques include quantum versions of kernel methods \cite{Schuld2018a, SchuldSQMLareKernel}, support vector machines \cite{Rebentrost2013, Liu2020f}, generative adversarial networks \cite{DallaireDemers2018, Lloyd2018}, Born machines \cite{coyle2020born}, Boltzmann machines \cite{Amin2016}, quantum neural networks \cite{Mitarai2018, Benedetti2019, Abbas2021} and more \cite{Biamonte}.\\

One regression task, fitting a continuous function $f(x)$ to a dataset of scalar-valued points $(x_i, y_i)$ in k-dimensional space, was analyzed for example in \cite{Mitarai2018}. A quantum neural network, consisting of a quantum feature map circuit, a variational quantum circuit and a Hamiltonian operator was used to represent a parameterizable function $f(x)$ that could be evaluated after training using gradient descent on analytically computed gradients of the loss function with respect to the variational circuit parameters.\\

In Ref. \cite{Kyriienko2020}, Differentiable Quantum Circuits (DQCs) paradigm was introduced. A DQC
implements a quantum neural network that can be analytically differentiated with respect to variational parameters \emph{and} with respect to the network inputs.
With automatic differentiation w.r.t. inputs it became possible to solve differential equations in a variational approach. In the current contribution, we will briefly summarize how this is done and we draw more parallels to classical scientific machine learning.\\

\paragraph{Physics-Informed Machine Learning}%
\label{par:physics_informed_machine_learning}

Physics informed neural networks (PINNs) comprise a a recent technique in scientific machine learning of particular relevance to this current work, since they bear strong conceptual similarities to DQCs.
We see both DQCs and PINNs as universal function approximators (UFAs) which can be differentiated with respect to their inputs. Furthermore, both can be optimized so that such derivatives can be manipulated.
Consequently, we expect that future applications of PINNs may well translate naturally to a DQC setting. Indeed, part of this current work is to demonstrate how such cross-fertilization between the fields might take place.

Introduced in \cite{Lagaris1997}, the essence of PINNs is to use (higher-order) derivatives of neural network outputs with respect to neural network inputs as terms with which to train neural networks.
Consequently, neural networks can be trained to satisfy a given differential equation.
Recently published work \cite{Raissi2019} has cause a surge in popularity of PINNs.
For example, applications of PINNs have been found in diverse fields such as aerodynamic flow modelling \cite{Mao2020} and cardiac activation mapping \cite{Costabal2020}.
Apart from solving forward problems, they have also found use in solving inverse problems \cite{Mao2020} and optimal control \cite{Wang2021}.
See \cite{KarniadakisNature2021} for a more complete recent review of PINNs and their applications.

Due to the conceptual similarity betweens PINNs and DQCs, we anticipate ample scope for cross-pollination of ideas and techniques between these lines of research.
Of particular relevance to this current work is the application of PINNs to equation discovery \cite{Both2021}, which we cover in more detail in subsequent sections.

\paragraph{Classical Equation Discovery}%
\label{par:classical_equation_discovery}
A popular technique in equation discovery is to overparameterize a candidate differential equation solution in terms of a library of potential terms, and then fit to data while encouraging sparsity in differential equation coefficients \cite{Brunton2016a} through an appropriate regularization term.
The resulting non-zero coefficients provide an interpretable differential equation representation of a process which fits the data.
While popularized in the context of ordinary differential equations (ODEs) \cite{Brunton2016a}, such techniques have been extended to the context of partial differential equations (PDEs) \cite{Rudy2017}.

Such sparsity-driven equation discovery schemes have recently been introduced to PINNs \cite{Both2021}.
We are not aware of any approaches to equation discovery which levereage quantum computation.
Given the conceptual similarity between PINNs and DQCs, it is natural to extend such equation discovery techniques to a quantum setting, which is the purpose of this work.

Before we introduce our approach to \emph{quantum model discovery} (QMoD\xspace \cite{patent}) in Sec.~\ref{sub:discovering_diffeq} we briefly introduce the DQC strategy that uses \textit{quantum} neural networks as universal function approximators in Sec.~\ref{sub:quantum}.

\paragraph{Our Contributions.}%
\label{par:our_contributions}

Our main contribution is to extend DQCs to wider applications in scientific machine learning, namely:
\begin{itemize}
    \item parameter inference in differential equations
    \item discovery of unknown differential equations from empirical measurements
\end{itemize}

We demonstrate our approaches on two exemplary ordinary differential equations (ODE) and on one partial differential equation (PDE).

\section{Methods}%
\label{sec:method}
Differential equations lie at the heart of many engineering disciplines.
However, using traditional, mesh-based methods such as finite-differencing or
finite-elements do not lend themselves well to several important applications, such as optimal control or solving inverse problems, since such applications require many repeated simulations.
Therefore, it is desirable to find cost-effective alternatives such as
data-driven differentiable surrogate models which can be optimized via gradient-based means.
On classical hardware one way of obtaining differentiable surrogates that respect physical prior
information is via \emph{physics-informed neural networks} (PINNs).
They can be used to solve PDEs given an equation, boundary conditions, and (optionally) some measurements.
Our work builds on top of the quantum PINN analogue \citep{Kyriienko2020} which
we briefly describe in Sec.~\ref{sub:solving_diffeq}.

We then demonstrate how further scientific machine learning techniques can be extended to a quantum setting.
In Sec.~\ref{sub:inferring_diffeq}, we introduce quantum methods for parameter inference in differential equations. Sec.~\ref{sub:discovering_diffeq} describes how to \emph{discover} equations from measurements on quantum computers. By discovering we mean finding a unknown differential equation only
from a given dataset and a predetermined (large) library of possible terms.

\subsection{Differentiable Quantum Universal Function Approximators}\label{sub:quantum}
In DQC, quantum neural networks act as differentiable universal function approximators. These are
quantum circuits consisting of three components: a quantum feature map, a variational ansatz, and a Hamiltonian operator for readout/encoding of the function.
These elements act on an array of qubits and alter the system wavefunction in a specific way.

The feature map encodes the input variable $x$ using a predefined non-linearity $\varphi$
to the amplitudes of a quantum state $\hat{\mathcal U}_\varphi(x)\ket{0}^n$, where $\ket{0^n}$ is some initially
prepared state. In the simplest case this can be implemented with a single layer
of rotations
\begin{equation}
  \hat{\mathcal U}_\varphi(x) = \bigotimes_{j=1}^n \hat R_{\alpha,j}(\varphi_j(x)).
\end{equation}
Furthermore, as an example, choosing the non-linearity as
\begin{equation}
  \varphi_j(x) = j \cdot \text{arccos}(x)
\end{equation}
results in a Chebyshev basis in the amplitudes of $\hat{\mathcal U}_\varphi(x)\ket{0^n}$. The state now represents a basis
of $n$ Chebyshev polynomials which can be read out by measuring for example in the Z basis on each qubit.

With a subsequent variational circuit (with parameters $\bm\theta$)
the Chebyshev basis functions are combined; in the limit of large and controlled entanglement, such a setup allows access to up to $\mathcal{O}(2^n)$
Chebyshev polynomials, due to the chaining and nesting properties of products of these polynomials.
Similar effects occur for Fourier-type feature maps and other universal basis function sets.
A popular choice for $\hat{\mathcal U}_\theta$ is the so-called hardware efficient ansatz (HEA) with parameterized rotations
and CNOT gates between neighboring qubits on the device.

The output $\hat u(x,\bm\theta)$ of the circuit is then computed with a Hamiltonian operator of choice, $\hat{\mathcal C}$,
such that the final universal approximator for a scalar input variable $x$ and variational
parameters $\bm\theta$ then reads
\begin{equation}
  \label{eq:ufa}
  \hat u(x,\bm\theta) = \bra{0^n}\hat{\mathcal U}_\varphi^\dagger(x)\hat{\mathcal U}_\theta^\dagger\hat{\mathcal C}\hat{\mathcal U}_\theta\hat{\mathcal U}_\varphi(x)\ket{0^n}
\end{equation}
Under some conditions, equation~\ref{eq:ufa} is fully analytically differentiable via the parameter
shift rule~\citep{Mitarai2018} which requires two evaluations of the same circuit
with shifted angles for each circuit-level parameter. As a particular example, where the feature map
consists of single-qubit rotations, we can write the derivative w.r.t. $x$ as follows
\newcommand{\hat{\mathcal U}_{\varphi,j,j'}^\pm(x)}{\hat{\mathcal U}_{\varphi,j,j'}^\pm(x)}
\newcommand{\hat{\mathcal U}_{\varphi,j,j'}^{\pm,\dagger}(x)}{\hat{\mathcal U}_{\varphi,j,j'}^{\pm,\dagger}(x)}
\begin{align}
  \label{eq:paramshift}
  \frac{d}{dx}\hat u(x,\bm\theta) &= \frac{1}{2}\frac{d\varphi}{dx} (\braket\hat{\mathcal C}^+ - \braket\hat{\mathcal C}^-),\\
  \braket\hat{\mathcal C}^\pm &= \sum_{j'} \bra{0^n} \hat{\mathcal U}_{\varphi,j,j'}^{\pm,\dagger}(x)  \hat{\mathcal U}_\theta^\dagger  \hat{\mathcal C}  \hat{\mathcal U}_\theta  \hat{\mathcal U}_{\varphi,j,j'}^\pm(x) \ket{0^n},\\
  \hat{\mathcal U}_{\varphi,j,j'}^\pm(x) &=\bigotimes_{j=1}^n  R_{y,j}(\varphi(x)\pm\frac{\pi}{2}\delta_{j,j'}).
\end{align}

Please note that we so-far considered a single-dimensional independent variable $x$. However, the DQC technique can straightforwardly be extended to multiple variables $x, y, t$ etc by using multiple feature maps, one for each variable, and differentiating only those feature maps relevant to the variable differentiation of interest. An example application where this is shown in Ref. \cite{paine2021}.

Furthermore, the example shows the regular parameter shift rule, which only applies to involutory
quantum generators in the feature map. Generalized parameter shift rules exist for arbitrary generators, and hence for arbitrary feature map circuits \cite{kyriienko2021}. In some cases, more intricate feature maps can be highly beneficial for increased expressivity of the quantum universal approximator.

\subsection{Solving Differential Equations}%
\label{sub:solving_diffeq}
We begin by recapping how to use automatically differentiable universal function approximators (either quantum or classical) to solve partial differential equations.
This section unifies existing treatment of using neural networks to solve differential equations \cite{Raissi2019} and analogous quantum methods \cite{Kyriienko2020}.
We write a partial differential equation in the form
\begin{equation}
  \label{eq:de}
  F(u, \partial_t u, \partial_x u, \partial^2_x u,\,\dots\,, x, t) = 0,
\end{equation}
where $x$ and $t$ are dependent variables and $u(x,t)$ is the solution we seek.
For simplicity we will further refer to the LHS of equation~\ref{eq:de} as
$F(u,x,t)$ omitting terms like $\partial_x u$ or $\partial_t u$ which in principle are operators acting on $u$.
Let's imagine we had access to a trial function $u(x,t)$ and we would like to ask how well
it solves the differential equation system. To find the quality of an approximate solution
$\hat u(x,t,\bm\theta)$ parameterized by $\bm\theta$ we can formulate a loss function
\begin{equation}
\label{eq:diffeqloss}
  \mathcal L(\bm\theta) = \mathcal L_F + \mathcal L_\Omega,
\end{equation}
where $\mathcal L_F$ is derived from the differential equation \ref{eq:de}
\begin{equation}
  \mathcal L_F = \frac{1}{M} \sum_{i=1}^M
                 \norm{F(\hat u,x_i,t_i)}_2.
\end{equation}
and is thus minimized $F(u,x,t) = 0$, i.e. when $\hat u(x,t,\bm\theta)$ satisfies the partial differential equation dynamics.
The second loss term $\mathcal L_\Omega$ includes all boundary conditions which
have to be met in order to fully specify the DE.\\

In order to effectively train such a model we ideally would like a universal approximator $\hat u$
and access to its higher order derivatives both with respect to inputs $x$ and $t$,
as well as parameters $\bm\theta$. Recent advances in automatic differentiation~\citep{Gune2018}
have facilitated this using neural networks on classical hardware.
On quantum hardware we leverage \emph{differentiable quantum circuits} (DQC)~\citep{Kyriienko2020} which implement a universal approximator based on a quantum featuremap $\hat{\mathcal U}_\varphi(x)$ and a variational circuit $\hat{\mathcal U}_\theta$.

\subsection{Inferring Differential Equation Parameters}
\label{sub:inferring_diffeq}
As a first demonstration of how scientific machine learning can translate to novel quantum algorithms, we consider parameter inference in differential equations.
In many situations, an overall structure of a differential equation is known, but specific coefficients must be inferred from data.
For example, one might parameterize a given differential equation of some scalar variable $u$ with some parameters $\bm w$ to be identified:
\begin{equation}
  \label{eq:inferdewithweights}
  F_{\bm w}(u,x,t) = 0,
\end{equation}
where $F=\bm{w}_{1} \partial^2_x u(x,t) + \bm{w}_2 x$ is one example, and the coefficients $\bm w$ must be inferred from some given data set. We denote a given set of $N$ data points with $\{(x_i, t_i, y_i): i=1,2\dots N \}$, where $x_i$ and $t_i$ denote independent variables and $y_i$ is a scalar dependent variable.   We assume that $F$ is constructed so that there is no $\bm w$ such that $F=0$ for every choice of $u$ to avoid trivial solutions.

For example, to infer a one-dimensional diffusion coefficient one might define $F$ as follows:

\begin{equation}
  F_{\bm w} = \bm w \partial^2_x u - \partial_t u,
\end{equation}

where we have been mindful not to introduce a trainable coefficient to $\partial_t u$, which would allow \eqref{eq:inferdewithweights} to be satisfied trivially with coefficients as zero.

To infer parameters in such settings, we introduce a loss function

\begin{equation}
  \label{eq:inferdeloss}
  \mathcal L(\bm\theta, \bm w) = \mathcal L_F + \mathcal L_d.
\end{equation}

Note that unlike \eqref{eq:diffeqloss}, we do not include a term analagous to $\mathcal L_\Omega$ to represent boundary conditions.
Although such terms could be introduced if boundary conditions are known, in principle the dynamics of the interior of the domain are inferrable from the data alone, and consequently boundary effects are not necessarily important to include.

In this loss function,

\begin{equation}
  \mathcal L_F = \frac{1}{M} \sum_{i=1}^M
                 \norm{F_{\bm w}(\hat u,x_i,t_i)}_2
\end{equation}

ensures that the relevant differential equation dynamics are captured and

\begin{equation}
  \mathcal L_d = \frac{1}{N} \sum_{i=1}^N
                 \norm{\hat u(x_i,t_i,\bm\theta) - y_i}_2
\end{equation}

provides incentive for the empirical dataset to be fit.

Note that we treat the loss function as a function of both the variational parameters, $\bm\theta$, as well as the inferred parameters, $\bm w$.
In our experiments, we optimize both these parameters on an equal footing, i.e. they are trained with the same learning rate and same optimizer, but in principle a more general optimization scheme also suffices.

\subsection{Discovering Differential Equations}%
\label{sub:discovering_diffeq}

Instead of solving a \emph{known} differential equation this section describes
how to find an \emph{unknown} DE just from measurements and a library of plausible
basis functions.

As in section~\ref{sub:inferring_diffeq}, we start by introducing a parameterized differential equation which can be represented in a general form

\begin{equation}
  \label{eq:discoverdewithweights}
  F(u,x,t) = 0,
\end{equation}
where again the coefficients $\bm w$ must be inferred from some given data set. However, we parameterize $F$ such that it represents a large family of PDEs of which we expect a given instance to fit the observed data.
In particular, in the context of this work, we draw our attention to dynamical systems of the form
\begin{equation}
  \label{eq:regprob}
  \partial_t u = \bm w \cdot \bm\varphi(u),
\end{equation}
where $u$ is a scalar dependent variable and $\bm\varphi(u)$ is a library of
potentially important functions such as
\begin{equation}
  \label{eq:basis}
  \bm\varphi(u) = [1, x, t, u, u^2, u^3, \partial_x u, \partial^2_x u, u\partial_x u, \dots]^T.
\end{equation}
The library $\bm\varphi$ can be large if there is little prior knowledge about the
terms that are present in the data, or contain only few specific terms that are
likely to have generated the given dataset.
Finding a solution that is \emph{sparse} in $\bm w$ results in a simple expression
for the right-hand side (RHS) of equation~\ref{eq:regprob}. Hence, we can perform
equation discovery by finding a sparse coefficient vector $\bm w$.

Equation~\ref{eq:regprob} can be solved via STRidge regression (\cite{Brunton2016a}) if one
has access to measurements of both $u$ and $\partial_t u$ which is rare\footnote{If
measurements of $\partial_t u$ are not directly available they can be computed e.g.  via
the Savitzky-Golay filter (\cite{Schaferb}), but such methods typically
require very dense measurements.}.
\emph{Physics-informed neural networks} (PINN by \cite{Raissi2019}) which are normally used
for \emph{solving} DEs
can circumvent this problem by learning a surrogate $\hat u$ as described in Sec.~\ref{sub:solving_diffeq}.
The surrogate is fully differentiable and therefore gives easy access to the
derivatives needed to solve equation~\ref{eq:regprob}.
Combining PINNs and equation~\ref{eq:regprob} to DeepMoD (\cite{Both2021}) results in
a fully automated framework that can discover DEs purely from data and a
pre-specified library of basis functions $\bm\varphi(u)$.

On classical hardware neural networks are used to represent the approximate solution $\hat u$.
To perform equation discovery on NISQ devices we replace the classical neural
networks with quantum universal function approximators implemented on DQCs resulting
in our \emph{Quantum Model Discovery} (QMoD\xspace) approach \cite{patent}. This approach is schematically depicted in figure~\ref{fig:qmd_schematic}.

The QMoD\xspace loss is composed of three terms: the data-loss $\mathcal L_d$, the
DE (or physics) loss $\mathcal L_p$, and the regularization:
\begin{equation}
  \label{eq:discoverdeloss}
  \mathcal L(\bm\theta, \bm w) = \lambda_d \mathcal L_d + \lambda_p\mathcal L_p + \lambda_r\norm{\bm w}_1.
\end{equation}

Note, that like equation \ref{eq:inferdeloss}, we do not include a term to impose boundary conditions.
If boundary conditions are known, such a term might be introduced.
However, we focus on situations where the interior PDE dynamics are inferred from data alone.
Also, as in equation \ref{eq:inferdeloss}, we explicitly write the loss in terms of the variational parameters of the UFA ($\bm\theta$), and the coefficients of the learnt differential equation ($\bm w$).
We optimize over both these parameters on an equal footing, using the same optimizer and learning rate.
The term $\lambda_r\norm{\bm w}_1$ encourages sparsity in the solutions of $\bm w$ and is standard practice in statistical learning \cite{Hastie2009}.

The data-loss is the standard mean squared error (MSE)
\begin{equation}
  \label{eq:mse}
  \mathcal L_d = \frac{1}{N} \sum_{i=1}^N
                 \norm{\hat u(x_i,t_i,\bm\theta) - y_i}_2,
\end{equation}
running over $N$ data points with temporal, spatial and scalar responses given by $t_i$, $x_i$ and $y_i$ respectively.
The differential equation loss is computed via the linear combination of basis functions
\begin{equation}
  \mathcal L_p = \frac{1}{M} \sum_{i=1}^M
                 \norm{\partial_t\uh(x_i,t_i,\bm\theta) - \bm w \bm\varphi(\hat u)}_2.
\end{equation}
The points $(x_i,t_i)$ at which the DE loss is evaluated can be chosen
independently from the data grid of equation~\ref{eq:mse} via a fixed, random, or
adaptive grid which gives this method great flexibility.

Generalizing QMoD\xspace to higher dimensional $\bm u = [u_1,u_2,\,\dots]^T$ only
requires taking care of matching dimensions (e.g. the coefficient vector $\bm
w$ turns into a matrix). In the equation discovery setting it makes sense to
prevent the surrogates for the variables $u_i$ from interaction, because we
want the interaction of terms to occur only via the library $\bm\varphi$.
Therefore, we use one individual UFA per operator instead of one large UFA
with multiple outputs, that would allow interaction outside the library.

\section{Experiments}%
\label{sec:experiments}

In this section, we demonstrate quantum scientific machine learning to be a viable technique.
In addition to demonstrating parameter inference, we also compare our QMoD\xspace to its classical, neural network based
counter part (DeepMoD).
We generate data from several different differential equations and demonstrate applications of QMoD\xspace.

\subsection{Damped Harmonic Oscillator}%
\label{sub:damped_harmonic_oscillator}

We start by demonstrating the potential of quantum neural networks for parameter inference and equation discovery in second-order ODES.

We parameterize a damped harmonic oscillator as follows:
\begin{equation}
  \frac{d^2 u_d}{dt_d^2} = -\omega^2 u_d -\alpha\frac{du_d}{dt_d},
\end{equation}
where $u_d$ is the position, $\omega$ the frequency, and $\alpha$ the damping coefficient
of the oscillator. The index $d$ indicates that this is the ODE before non-dimensionalization.
We non-dimensionalize the equation such that the data lies in within the range
$(-1,1)$ by defining $u_d = u u_c$ and $t_d=t t_c$ which results in
\begin{equation}
  \label{eq:damped-harmonic-ode}
  \frac{d^2u}{dt^2} = -\omega^2 t_c^2 u -\alpha t_c^2\frac{du}{dt},
\end{equation}
With the constants $u_c=1\text{\,m}$, $t_c=2\pi\text{\,s}$, $\omega=1.5\frac{1}{\text{s}}$, and
$\alpha=0.5\frac{1}{\text{s}}$ we obtain training data as shown on the left plot of figure~\ref{fig:damped-harmonic-param-id}.

\subsubsection{Inferring Differential Equation Parameters}%
\label{ssub:parameter_identification}

As a first step we assume that the form the governing differential equation is
known and we must infer the values of the constants
$\omega$ and $\alpha$.
We compare the method outlined in section~\ref{sub:inferring_diffeq} with an implementation of DeepMoD without regularization coefficients.
The optimization problem in \eqref{eq:inferdeloss} with the specialized form of \eqref{eq:inferdewithweights} is defined as:
\begin{equation}\label{eq:varphi-param-id}
  F_{\bm w}(u,x,t) =  \frac{d^2u}{dt^2} + \bm w \left[t_c^2 u, t_c^2\frac{du}{dt}\right]^T.
\end{equation}
We include the known non-dimensionalization constants in the basis functions to scale the problem so that it is more amenable
to standard optimization hyperparameters (not doing this would result in the models having to identify coefficients $\gg 10$).

Figure~\ref{fig:damped-harmonic-param-id} shows the data and the obtained QMoD\xspace fit $\hat u$
on the left. We use a quantum neural network with 5 qubits as the UFA to approximate $\hat u$. We
reach an MSE loss value of $\mathcal L_d=10^{-7}$ after 200 BFGS optimization steps.  A
neural network with one hidden layer of size 32 and $\tanh$ activations reaches
similar MSE.
The central and right plot of figure~\ref{fig:damped-harmonic-param-id} show the
trajectories of the coefficients of the two basis functions
(equation~\ref{eq:varphi-param-id}), namely $\omega$ and $\alpha$, during training.
The target values for the two coefficients are $\omega=1.5$ and $\alpha=0.5$.

In this experiment we found good solutions using $\lambda_d=1$ (dataloss),
$\lambda_p=10^{-5}$ (physics loss) for both QMoD\xspace and DeepMoD within \eqref{eq:inferdeloss}.
The identified parameter values are shown in table~\ref{tab:damped_harmonic_param_id}.
\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{\input{damped-harmonic-param-id.tikz}}

  \caption{Damped harmonic parameter \textit{identification}. \figuretext
  Both QMoD\xspace and the classical DeepMoD infer the correct target parameters
  $\omega^2=2.25$ and $\alpha=0.5$.
  }%
  \label{fig:damped-harmonic-param-id}
\end{figure}
\begin{table}
  \label{tab:damped_harmonic_param_id}
  \centering
  \begin{tabular}{llll}
    \toprule
             &       & \multicolumn{2}{c}{Basis function} \\
    \cmidrule{3-4}
    Equation & Model & $u$  & $\frac{du}{dt}$ \\
    \midrule
    \multirow{3}{*}{$\frac{d^2u}{dt^2} = -\omega^2 t_c^2 u -\alpha t_c^2\frac{du}{dt}$}
    & QMoD\xspace      & $\omega=\sqrt{2.249} = 1.499$               & $\alpha=0.500$ \\
    & DeepMoD   & $\omega=\sqrt{2.247} = 1.498$               & $\alpha=0.500$ \\
    & Truth     & $\omega=1.5$                                & $\alpha=0.5$ \\
    \bottomrule
  \end{tabular}
  \caption{Final basis function parameters $\omega$ and $\alpha$ as inferred by QMoD\xspace
    and DeepMoD compared to the true coefficients that generated the data (called
    \emph{truth}).
    Non-dimensionalization coefficients in the basis functions are omitted for clarity.}
\end{table}

\subsubsection{Equation Discovery}%
\label{ssub:equation_discovery}

To go a step beyond the task of parameter inference, we next aim for full "equation discovery"; we
assume that only the LHS of equation~\ref{eq:damped-harmonic-ode} (i.e. $\partial u_t^2$) is known and choose
a set of plausible basis functions $\bm\varphi$ consisting of polynomials up to
degree 4 and the first derivative of $u$ for use in \eqref{eq:regprob}.
\begin{equation}\label{eq:varphioscillator}
  \bm\varphi(u) = \left[t_c^2\, u(t),\quad
                        u_c t_c^2\, u^2(t),\quad
                        u_c^2t_c^2\, u^3(t),\quad
                        u_c^3t_c^2\, u^4(t),\quad
                        t_c \frac{du}{dt}
                  \right]^T
\end{equation}
\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{\input{damped-harmonic.tikz}}

  \caption{Damped harmonic oscillator \emph{discovery}. \figuretext
    The difference to
    figure~\ref{fig:damped-harmonic-param-id} is the larger set of basis
    functions.  Both QMoD\xspace and DeepMoD correctly identify the relevant basis
    functions ($u$ and $\frac{du}{dt}$) whereas the irrelevant basis function
    coefficients are all close to zero.
  }%
  \label{fig:damped-harmonic}
\end{figure}
\newcommand{$<10^{-3}$}{$<10^{-3}$}
\begin{table}
  \centering
  \begin{tabular}{lllcccl}
    \toprule
             &       & \multicolumn{5}{c}{Basis function} \\
    \cmidrule{3-7}
    Equation & Model & $u$ & $u^2$ & $u^3$ & $u^4$ & $\frac{du}{dt}$ \\
    \midrule
    \multirow{3}{*}{$\frac{d^2u}{dt^2} = -\omega^2 t_c^2 u -\alpha t_c^2\frac{du}{dt}$}
    & QMoD\xspace    & $\omega=\sqrt{2.239} = 1.496$ &$<10^{-3}$&$<10^{-3}$&$<10^{-3}$& $\alpha=0.497$ \\
    & DeepMoD & $\omega=\sqrt{2.242} = 1.498$ &$<10^{-3}$&$<10^{-3}$&$<10^{-3}$& $\alpha=0.496$ \\
    & Truth   & $\omega= 1.5$ &$0$&$0$&$0$& $\alpha=0.5$ \\
    \bottomrule
  \end{tabular}
  \label{tab:damped_harmonic}
  \caption{Final basis function parameters resulting from QMoD\xspace/DeepMoD training with
    a larger library of basis functions (equation~\ref{eq:varphioscillator}).
    All irrelevant coefficients are very close to zero while the inferred $\omega$
    and $\alpha$ are close to their true values.
    Non-dimensionalization coefficients in the basis functions are omitted for clarity.
  }
\end{table}
We optimize \eqref{eq:discoverdeloss} and achieve good results with a regularization parameter $\lambda_r=10^{-4}$.
Otherwise we use the same training and model setup (5 qubits, 200 BFGS steps,
$\lambda_p=10^{-5}$) as in section~\ref{ssub:parameter_identification}.
Both our quantum equation discovery QMoD\xspace and DeepMoD successfully identify the
relevant basis functions and their correct coefficients. All other basis function
coefficients are close to zero after training.

\subsection{Lotka-Volterra System}%
\label{sub:lotka_volterra_system}
We now demonstrate that coupled non-linear differential equations can be discovered with quantum circuits.
Our model system is given by the Lotka-Volterra (LV) equations:
\begin{align}
  \label{eq:lotka-volterra}
  \frac{dx_d}{dt_d} &= \alpha x_d - \beta x_dy_d\\
  \frac{dy_d}{dt_d} &= \delta x_dy_d -  \gamma y_d
\end{align}
that describe a predator-prey system where $x_d$ is the number of prey and
$y_d$ the number of predators.  The parameters $\alpha$, $\beta$, $\gamma$, and
$\delta$ define the interaction of the two species which we chose to
$\alpha=1.5$, $\beta=1$, $\gamma=1$, $\delta=1$.
Again, the index $d$ indicates dimensional quantities. We non-dimensionalize
the LV system with critical time $t_c=8$ and a critical number of individuals
$x_c=2.5$ which results in
\begin{align}
  \label{eq:lotka-volterra-nd}
  \frac{dx}{dt} &= \alpha t_c x - t_c x_c \beta xy\\
  \frac{dy}{dt} &= t_c x_c \delta xy -  t_c\gamma y.
\end{align}
The non-dimensionalization coefficients are again included in the basis
functions $\bm\varphi$ which consist of polynomials up to degree 2
\begin{equation}
  \bm\varphi(x,y) = \left[t_c x,\, t_c y,\, t_c x_c x^2 ,\, t_c x_c y^2 ,\, t_c x_c xy\right]^T.
\end{equation}

QMoD\xspace is trained with a time-series of 30 datapoints for each variable $x$ and
$y$.  For this task we need a slightly larger featuremap with 7 qubits to fit the data.
Each variable of the LV system is modeled by a separate UFA driving its own set of (identical) library functions.
After 200 BFGS steps we discover the correct basis
function coefficients as shown in figure~\ref{fig:lotka-volterra} and
table~\ref{tab:lotka-volterra}. The QMoD\xspace loss coefficients are $\lambda_d=1$,
$\lambda_p=10^{-3}$, and $\lambda_r=10^{-5}$.

We compare QMoD\xspace to classical DeepMoD with two neural networks (one for each
variable $x$ and $y$) with one hidden layer of size 32.
Classical DeepMoD receives the same data. For this dataset we observe that the
classical neural networks suffer from stronger gradient pathologies which could
only be resolved by decreasing the physics loss coefficient to $\lambda_p=10^{-6}$
which weakens the training signal on the basis function coefficients.
This explains why DeepMoD needs more training steps to converge (see figure~\ref{fig:lotka-volterra}).
\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{\input{lotka-volterra.tikz}}

  \caption{Lotka-Volterra discovery. \figuretext
  Both QMoD\xspace and DeepMoD
  correctly identify the relevant basis functions for each variable.}%
  \label{fig:lotka-volterra}
\end{figure}
\begin{table}
  \label{tab:lotka-volterra}
  \centering
  \begin{tabular}{llccccc}
    \toprule
             &       & \multicolumn{5}{c}{Basis function} \\
    \cmidrule{3-7}
    Equation & Model & $x$  & $y$ & $x^2$ & $y^2$ & $xy$ \\
    \midrule
    \multirow{3}{*}{$\frac{dx}{dt} = \alpha t_c x - t_c x_c \beta xy$}
    & QMoD\xspace      & \caleft{$\alpha = 1.483$}    & $0.001$           & $-0.001$   & $-0.007$   & \caleft{$\beta=0.986$} \\
    & DeepMoD   & \caleft{$\alpha = 1.475$}    & $<10^{-3}$        & $<10^{-3}$ & $-0.003$   & \caleft{$\beta=0.979$} \\
    & Truth     & \caleft{$\alpha = 1.5$}      & $0$               & $0$        & $0$        & \caleft{$\beta=1$}\\
    \midrule
    \multirow{3}{*}{$\frac{dy}{dt} = t_c x_c \delta xy -  t_c\gamma y$}
    & QMoD\xspace      & $0.006$             & \caleft{$\delta = 0.993$}  & $0.008$    & $<10^{-3}$ & \caleft{$\gamma=0.991$} \\
    & DeepMoD   & $<10^{-3}$          & \caleft{$\delta = 0.987$}  & $0.011$    & $<10^{-3}$ & \caleft{$\gamma=0.979$} \\
    & Truth     & $0$                 & \caleft{$\delta = 1$}      & $0$        & $0$        & \caleft{$\gamma=1$}\\
    \bottomrule
  \end{tabular}
  \caption{Final basis function parameters of the LV system discovered by QMoD\xspace and DeepMoD
    compared to the true coefficients that generated the data (called
    \emph{truth}).  Non-dimensionalization coefficients in the basis functions
    are omitted for clarity. Both models identify the correct basis functions as relevant.}
\end{table}

\newpage
\subsection{Burgers' Equation}%
\label{sub:burgers_equation}

In this last section we present the discovery of Burgers' equation,
a non-linear second-order partial differential equation.
We chose the general (or \emph{viscous}) Burgers equation with one spatial
dimension $x$ and one time dimension $t$:
\begin{equation}
  \frac{\partial u}{\partial t} = \varepsilon\frac{\partial^2u}{\partial x^2} - \alpha u\frac{\partial u}{\partial x},
\end{equation}
where $u(x,t)$ is the velocity field, $\varepsilon$ the viscosity, and $\alpha$
the coefficient of the advection term. To obtain training data we solve
Burgers' equation on the domain $(x \times t)\in([0,1]\times[0,\frac{1}{2}])$
with parameters $\varepsilon=0.05$ and $\alpha=1$. The true solution with sinusoidal
initial conditions
\begin{equation}
  u(x,0) = \sin(2\pi x)
\end{equation}
is shown in the left plot of figure~\ref{fig:burgers-heatmap} as the gray-scale contour. We train on 200
uniformly sampled points\footnote{Note that no intricate grid definition is
needed for DQC/QMoD\xspace to solve PDEs} which are scattered on top of the left plot of figure~\ref{fig:burgers-heatmap}. For Burgers' equation we do not need to
perform non-dimensionalization as the domain is already within $(-1,1)$.
\begin{figure}
  \centering
  \includegraphics[width=\linewidth]{burgers-heatmap.pdf}
  \caption{\emph{Left}: Scatter plot of Burgers' equation dataset and the full solution as a gray-scale contour underneath
  (parameters $\varepsilon=0.05$, $\alpha=1$, and initial condition $u(x,0)=\sin(2\pi x)$).
  \emph{Center}: Learned QMoD\xspace solution. \emph{Right}: Error of QMoD\xspace solution.}%
  \label{fig:burgers-heatmap}
\end{figure}
We assume that the LHS $\partial u_t$ is known and want to discover $\varepsilon$
and $\alpha$ from the library
\begin{equation}
  \bm\varphi(u) = \left[u,\quad u^2,\quad \frac{\partial u}{\partial x},\quad
                        \frac{\partial^2u}{\partial x^2},\quad
                        u\frac{\partial u}{\partial x},\quad
                        u^2\frac{\partial u}{\partial x}
                  \right]^T.
\end{equation}
To solve this PDE we need a UFA with two inputs $(x,t)$ and one output.
Therefore, we employ two input featuremaps with 6 qubits each. The two featuremaps
are stacked in parallel onto one circuit, and connected to a variational circuit with a base on all 12 qubits, effectively entangling the two feature registers.
The physics loss coefficient is set to $\lambda_p=10^{-2}$ and the
regularization to $\lambda_r=10^{-3}$.

Figure~\ref{fig:burgers-discovery} and table~\ref{tab:burgers} show that
$\varepsilon$ is discovered very precisely but the advection term $\alpha$ is
slightly too low in both models.
However, we can still conclude that QMoD\xspace and DeepMoD correctly identify the two
relevant basis functions for advection and diffusion.
To get an even closer fit to the true basis function parameters it is possible to re-run
the model only with the two relevant basis functions an switch off the regularization (as we demonstrated
in section~\ref{ssub:parameter_identification}).
\begin{figure}
  \centering
  \resizebox{\textwidth}{!}{\input{burgers.tikz}}

  \caption{Burgers' equation discovery. \figuretext The trajectory of viscosity
    $\varepsilon$ (red line) is scaled by $\frac{1}{\varepsilon}$ for better
    visibility. Both QMoD\xspace and DeepMoD discover $\varepsilon$ with good
    precision. The advection term $\alpha$ (blue line) is clearly discovered as a relevant
    term but does not match the true value $\alpha=1$ perfectly.
  }%
  \label{fig:burgers-discovery}
\end{figure}
\begin{table}
  \centering
  \begin{tabular}{llcccllc}
    \toprule
             &       & \multicolumn{6}{c}{Basis function} \\
    \cmidrule{3-8}
    Equation & Model & $u$ & $u^2$ & $\frac{du}{dx}$ & $\frac{d^2u}{dx^2}$ & $u\frac{du}{dx}$ & $u^2\frac{du}{dx}$ \\
    \midrule
    \multirow{3}{*}{$\frac{\partial u}{\partial t} = \varepsilon\frac{\partial^2u}{\partial x^2} - \alpha u\frac{\partial u}{\partial x}$}
    & QMoD\xspace    & $<10^{-3}$ & $<10^{-3}$ & $<10^{-3}$ & $\varepsilon=0.050$ & $\alpha = 0.838$ & $<10^{-3}$ \\
    & DeepMoD & $<10^{-3}$ & $<10^{-3}$ & $<10^{-3}$ & $\varepsilon=0.049$ & $\alpha = 0.924$ & $<10^{-3}$ \\
    & Truth   & $0$     & $0$     & $0$     & $\varepsilon=0.05$  & $\alpha = 1$     & $<10^{-3}$ \\
    \bottomrule
  \end{tabular}
  \label{tab:burgers}
  \caption{Final basis function parameters for Burgers' equation discovered by QMoD\xspace
    and DeepMoD compared to the true coefficients that generated the data. Both QMoD\xspace and DeepMoD correctly identify that
    only two basis functions are needed to describe Burgers' equation. The parameters $\varepsilon$ is found very precisely
    by both models. DeepMoD is closer to the true parameter $\alpha$ but still slightly too low.}
\end{table}

\section{Discussion \& Conclusion}%
\label{sec:conclusion}
We demonstrated that techniques from scientific machine learning can be placed in a quantum computational setting.
In particular, our numerical experiments show that relatively small quantum circuits can be used for meaningful computations in scientific machine learning.
We have demonstrated how to implement an automated parameter inference and equation discovery
framework (QMoD\xspace) using a differentiable quantum circuit strategy, and reached results that are on
par with the classical machine learning method DeepMoD.
With this demonstration of what can be achieved with DQCs the fields of quantum
computing and scientific machine learning are moving closer together.

We should note that all experiments described here involve classical simulators of quantum computers with unrealistic assumptions such as zero sampling noise, no environment noise and no measurement noise. Although many error mitigation strategies for NISQ hardware exist, it should still be investigated how well the strategies described work in more realistic settings and importantly on real quantum hardware, and which quantum hardware is most compatible/natural for this strategy.
Furthermore, due to the intrinsic scaling issues of simulating quantum many-body systems, we did not consider more than 12 qubits in the largest simulations described here, and thus also compared with only small-sized classical neural networks. It will be important to investigate the performance for larger number of qubits on quantum hardware, and compare them to realistic classical PINN architectures.

The strategies described above work for systems governed by sets of Partial Differential Equations, but also for those governed by stochastic processes modelled as Stochastic (Partial) Differential Equations (SDEs), through the perspective of \textit{quantile mechanics} combined with DQCs as described in the recent paper on Quantum Quantile Mechanics (QQM) \cite{paine2021}. In such setting, an SDE may be learned (`discovered') based on just statistical data. This could prove valuable in the generative modelling case, where synthetic data generation of time-series is learned based on a combination of data and a parameterized SDE with a library of candidate model operators. In such real-world setting as financial markets, there is often not exact knowledge of the underlying model dynamics but there is access to historical data, where QMoD\xspace could infer the SDE from data and make better predictions compared to a purely data-only approach.

An interesting area of potential yet unexplored, is the identification of quantum computational speedup for scientific machine learning in other aspects than Neural Network enhancement, for example by combining QMoD\xspace with combinatorial optimization.

\textit{Ethics declaration.} A patent application for the method described in this manuscript has been submitted by Qu\&Co.

\printbibliography

\end{document}

redpajama_set_name: RedPajamaArXiv


text:
How the Most Vocal Legislators on Debt & Deficit Voted on the Budget
At 12:01 am on Friday morning the government shutdown. A vote to fund the government before midnight was delayed when Sen. Paul (R-KY) took to the floor to argue for stricter budget caps and the removal of the bill's debt-ceiling measure. An hour and 52 minutes later, the Senate would go on to pass the budget 71 to 28 and, after a brief standoff, the House followed suit at 5:30 am — 241 to 186. However, Sen. Paul is not the only lawmaker who's been raising flags on Congress' handling of the national debt and deficit over the past year. Here's a look at the top ten most vocal members on the debt and deficit during the 115th Congress alongside how they voted on yesterday's Bipartisan Budget Act of 2018:
Five of the top 10 most vocal senators on the debt and deficit voted in favor of the budget deal.
Combined, Minority Leader Schumer (D-NY) and Senators Perdue (R-GA), Brown (D-OH), Van Hollen (D-MD), and Durbin (D-IL) have mentioned the debt and deficit in 445 statements this Congress, and each voted in favor of the budget deal. Sen. Paul has been the second most vocal member in the chamber with 112 statements, succeeded by Sen. Sanders (I-VT) with 145.
Eight of the top 10 most vocal representatives on the debt and deficit are Democrats.
While coalitions of House Republicans, like the House Freedom Caucus, have historically been the most vocal on the debt and deficit, eight of the ten most vocal House members on the issue this Congress are Democrats. Freedom Caucus member Rep. Sanford (R-SC-1) cracked the top five with 115 statements. Reps. Kaptur (D-OH-9) and McCollum (D-MN-4) were the only members of the top 10 to vote in favor of the budget deal.

redpajama_set_name: RedPajamaCommonCrawl


text:
Q: Difference between declaration of function with extern and without it There is such code:
#include <iostream>

extern void fun();

int main(){
    fun();
    return 0;
}

void fun(){ std::cout << "Hello" << std::endl; }

Is there some difference between declarations:
extern void fun();
void fun();

? Code above behaves the same with extern and without extern keyword.

A: Function declarations do have external linkage by default, so adding the extern keyword to a function declaration makes no difference, it is redundant.

A: The difference between the two statements is this:
extern void fun();

tells the compiler and linker to look in another file when the code in this file refers to fun( ), perhaps by calling fun( ); This production is called a "declaration."
void fun ( ) {
  ...
}

Defines the function fun ( ) and, because it's defined in this file, obviates the need for the linker to look for the function elsewhere.
There's no harm in declaring the function extern: the linker does the right thing.

redpajama_set_name: RedPajamaStackExchange


An exciting range of original tableware and accessories. Ideal for complimenting your kitchen and updating it with fresh, modern elements.
Your kitchen is likely to be one of the most important spaces in your home. From kitchen cabinet choices to trays & cake stands, kitchen utensils and more besides – you're going to need to kit out your space with reliable, affordable fixtures and accessories you can buy online with ease.
Kitchens are both very sociable areas as well as where plenty of housework goes on – so you'll need access to specialist kitchen furniture design from leading local specialists to really get the most out of your space.
Need a new corner kitchen cabinet? Take a look at our quality kitchen table and chairs selection to welcome family and friends into the heart of your home – or why not consider a kitchen bar table for that modern touch and focal point you've been looking for?
Casanad are proud to be able to source leading local furniture design and manufacture for each and every room of the home. The kitchen, however, needs fixtures and fittings which emphasise comfort, style and functionality – often effortlessly.
If you're looking for cheap or discounted kitchen table sets or chinaware & glassware, let us show you the pinnacle of reliable and luxury standards – available for less than you'd expect to pay elsewhere on the high street.
We deliver low cost rates to our customers wherever possible, meaning that savings for you are always high on our list of priorities.
Our online store and furniture shop is stocked with a variety of luxury, reliable kitchen standards. From traditional folding kitchen table designs to round kitchen table units to really help make things sociable, simply take a close look at our fantastic range of tables, chairs and cabinets to start putting together your perfect kitchen suite.
Looked at Ikea kitchen table designs but are unsure about the look or the prices? Allow us to offer you a fantastic range of kitchen furniture choices for sale – at a low cost to suit all household budgets.
We're leading local suppliers and specialists in furniture throughout the home – and we're proud to be working with some of the most talented designers and manufacturers in kitchen fittings and accessories across London and beyond.
Our trusted outlet will be able to help you find affordable small kitchen tables, wine racks, jars, kitchen utensils and more besides – even if you're looking for a specialist kitchen table with bench or similar.
We have a wide selection of stunning items for you to choose from. Simply let us know if there is anything you're looking for in our catalogue – we're always happy to help!
When it comes to kitchen furniture and accessories, look for a local London company near me with a fantastic range – Casanad are here to help you redesign and restock your kitchen your way.
Looking to capture a certain look or to create your own unique space? Take a look at our kitchen catalogue and see what inspires you! Call us on 0203 795 1021 to learn more.


text:
Chicago Reporter (https://www.chicagoreporter.com/chicago-schools-policy-luncheon-march-20-2007/)
Chicago Schools Policy Luncheon: March 20, 2007
By Catalyst Chicago | April 1, 2007
Leslie Jacobs and Shenita Johnson Garrard describe how New Orleans is striving to convert a failed, corrupt school system into a decentralized system of schools. With some 60 percent of its public school students enrolled in charter schools, New Orleans has the most school autonomy and competition of any district in the country.