# Solving physics-based initial value problems with unsupervised machine
learning

Jack Griffiths  0000-0001-7794-6687<EMAIL_ADDRESS>Joint Quantum Centre (JQC)
Durham–Newcastle, Department of Physics, Durham University, Durham DH1 3LE,
United Kingdom  Steven A. Wrathmall  0000-0003-1770-9721
<EMAIL_ADDRESS>Joint Quantum Centre (JQC) Durham–Newcastle,
Department of Physics, Durham University, Durham DH1 3LE, United Kingdom
Simon A. Gardiner  0000-0001-5939-4612<EMAIL_ADDRESS>Joint Quantum
Centre (JQC) Durham–Newcastle, Department of Physics, Durham University,
Durham DH1 3LE, United Kingdom

###### Abstract

Initial value problems — a system of ordinary differential equations and
corresponding initial conditions — can be used to describe many physical
phenomena including those arise in classical mechanics. We have developed a
novel approach to solve physics-based initial value problems using
unsupervised machine learning. We propose a deep learning framework that
models the dynamics of a variety of mechanical systems through neural
networks. Our framework is flexible, allowing us to solve non-linear, coupled,
and chaotic dynamical systems. We demonstrate the effectiveness of our
approach on systems including a free particle, a particle in a gravitational
field, a classical pendulum, and the Hénon–Heiles system (a pair of coupled
harmonic oscillators with a non-linear perturbation, used in celestial
mechanics). Our results show that deep neural networks can successfully
approximate solutions to these problems, producing trajectories which conserve
physical properties such as energy and those with stationary action. We note
that probabilistic activation functions, as defined in this paper, are
required to learn any solutions of initial value problems in their strictest
sense, and we introduce coupled neural networks to learn solutions of coupled
systems.

machine learning; differential equations

††preprint: AIP/123-QED

##  I Introduction

Machine learning models embedded with physical constraints, equations, data,
or other physical principles (including variational principles such as the
Rayleigh–Ritz method or Hamilton’s principle, which is what we present in the
current work) offer a powerful framework for ensuring physically consistent
solutions for differential equations [1, 2, 3, 4, 5, 6, 7] or discovering new
physics [8, 9, 10, 11, 12, 13]. Incorporating informed priors in machine
learning models allows effective handling of sparse data, improved
generalisation, and offer deeper insights into the underlying mechanisms of
the systems they are designed to analyse.

Machine learning is an optimisation procedure. Machine learning is _physics-
informed_ if physical data or principles are present in i) the training data,
ii) the underlying neural network architecture, iii) the cost function and iv)
the optimisation algorithms. Training data may include data from physical
experiments or simulations [14, 15]. Training data may be augmented [16] by
applying transformations (e.g., translations, scalings or rotations) in line
with known symmetries of the system, without having to generate new data. The
network architecture may include structures which obey physical laws. Ling,
Kurzawski and Templeton [17] predict turbulent fluid flows using a deep neural
network with an auxiliary network which ensures invariance to Galilean
transformations (which one would expect from such classical fluid problems).

Regularisation terms may be added to the cost function to constrain the neural
network output based on, for example, the underlying differential equation,
whether the solution should be divergence free [18], or any other mathematical
description of the data [19]. The cost function therefore may consist of a
weighted combination of: i) the differential equation (or other equations
which constrain the solution), ii) the initial conditions, iii) the boundary
conditions (where appropriate), and iv) sparse data sampled from throughout
the problem domain (the training data). Other problem-specific conditions may
also be included in the cost function. Xu, _et al._ [6], solve Fokker–Planck
equations with machine learning methods, which additionally _requires_
constraining the normalisation of their solutions (which are probability
distributions). The training data may be sparse, as this approach effectively
interpolates between data points by enforcing the constraint given by the
differential equation.

We stress through our own experiments and through earlier research of Vapnik
[20, 21] that one should choose the simplest cost function possible to avoid
overfitting neural networks or — in some cases — complete failure of their
training. One may know the differential equation, initial conditions,
Hamiltonian [14], and Lagrangian [15], but one should not, in general, include
all of these terms in the cost function.

In the present work, we describe a framework to solve initial value problems
in the strictest sense, that is without any constraints which are not needed
mathematically. We do not include any extraneous regularisation terms as is
common in typical physics-informed neural networks [1, 2, 3, 4, 5, 6]. Our
framework is an implementation of Hamilton’s principle through machine
learning methods; for all dynamical problems we consider, the only constraints
are the underlying differential equation and its associated initial
conditions. There is no training data — this is _unsupervised_ learning. Our
discussion starts with the most minimally conceivable neural network to
describe free particle dynamics — this allows us to set notation and introduce
our neural network representation of initial value problems. We consider
linear approximations of a particle in a gravitational field and in a harmonic
oscillator. We briefly discuss observed deviations from the theoretical
results of the harmonic oscillator. Guided by the need to increase the
representational capacity of our neural networks, we introduce deep learning
techniques to solve the classical, arbitrary-angle pendulum. We introduce
coupled neural networks and optimisation routines, using the Hénon–Heiles
system [22] (akin to three-body dynamics [23]) as an example. We finally make
comments relating the machine learning trajectories back to Hamilton’s
principle — our solutions are eventually consistent with the principle of
stationary action, despite some peculiarities in how the action evolves during
training.

##  II Statement of the problem

###  II.1 Statement of the physical problem

Consider a mechanical system characterised by $m$ generalised coordinates
$\\{q_{i}(t)\\}$, $i\in\\{1,\cdots,m\\}$, and described by the Lagrangian $L$.
Application of Hamilton’s principle leads to the Euler–Lagrange equations

$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial
L}{\partial\dot{q}_{i}}-\frac{\partial L}{\partial{q}_{i}}=0.$ (1)

We chose to consider only Lagrangians of the form $L=T-V$, where

$T=\frac{1}{2}\sum_{i=1}^{m}\mu_{i}\dot{q}_{i}(t)^{2}$ (2)

is the kinetic energy and $V\equiv V(q_{1}(t),q_{2}(t),\cdots,q_{m}(t))$ is a
potential energy. We consider only potential energies that are explicit
functions of the generalised coordinates, which excludes time-dependent
potentials and the Lorentz force. The Euler–Lagrange equations lead to
equations of motion of the form

$\mu_{i}\ddot{q}_{i}(t)-F_{i}(q_{1}(t),q_{2}(t),\cdots,q_{m}(t))=0,$ (3)

where $\mu_{i}$ are generalised coefficients — these may be associated with a
mass or moment of inertia, depending upon the dynamical system — and $q_{i}$
are the generalised coordinates — which may be associated with positions or
angles or other suitable coordinates that uniquely define the configuration of
the system. Their dimensions will depend on the nature of each coordinate
(e.g., an angle coordinate will be dimensionless, whereas a length coordinate
will have dimensions of length). We can associate

$F_{i}(q_{1}(t),q_{2}(t),\cdots,q_{m}(t))=-\frac{\partial
V(q_{1}(t),q_{2}(t),\cdots,q_{m}(t))}{\partial q_{i}(t)},$ (4)

with a generalised force. In order to eliminate the explicit dependence upon
$\mu_{i}$, let the rescaled generalised force (which may have units of
acceleration or angular acceleration, as appropriate) be
$f_{i}=F_{i}/\mu_{i}$. This leads to the equations of motion

$\ddot{q}_{i}(t)-f_{i}(q_{1}(t),q_{2}(t),\cdots,q_{m}(t))=0,$ (5)

where we choose to cast the equations of motion in such a way that they equal
zero. This simplifies the subsequent development of the machine learning
problem in section II.2.

The domain of the problem is $t\in[0,T]$ (the time origin is set to zero
without loss of generality). The problem is subject to $2m$ initial
conditions, where each generalised coordinate is associated with an initial
position $s_{i,0}$ and velocity $v_{i,0}$.

###  II.2 Statement of the machine learning problem

We will denote the neural network’s representation of the solution by
$\mathbf{q}_{i}(\mathbf{t})=\left(q_{i,0},q_{i,1},\cdots,q_{i,n},\cdots,q_{i,N_{T}}\right)$.
Each element of the neural representation, $q_{i,n}$, is associated with the
continuous solution at a particular time $t_{n}=n\Delta t$, i.e.,
$q_{i,n}\approx q_{i}(n\Delta t)$, where $\Delta t=T/N_{T}$ is the time step.
We are seeking discretised vector representations
$\\{\mathbf{q}_{i}(\mathbf{t})\\}$ of the continuous functions
$\\{q_{i}(t)\\}$ over the temporal domain
$\mathbf{t}=(0,t_{1},\cdots,t_{n},\cdots,t_{N_{T}})$.

The machine learning problem is to find an optimised neural network
representation which approximates the solution of the initial value problem by
minimising an appropriate cost metric. Our cost metric will be a dimensionless
function (to not be dependent upon any particular system of units). We will
associate all timescales with a characteristic time $\tau$ and all length
scales with a characteristic length $\lambda_{i}$. This implies a
dimensionless generalised position coordinate
$\tilde{q}_{i}=q_{i}/\lambda_{i}$ and a dimensionless rescaled generalised
force $\tilde{f}_{i}=f_{i}\tau^{2}/\lambda_{i}$. Furthermore, the
dimensionless parameters $\alpha_{i},\beta_{i}$ and $\gamma_{i}$ may be
introduced to adjust the relative weighting of each term in the cost function.
Dropping all tildes, the dimensionless cost function we wish to minimise is

$\displaystyle\mathcal{C}$
$\displaystyle=\sum_{i=1}^{m}\frac{1}{\lambda_{i}^{2}}\Big{[}\alpha_{i}\tau^{4}\left\|\mathbf{\ddot{q}}_{i}(\mathbf{t})-\mathbf{f}_{i}\left(\mathbf{q}_{1}(\mathbf{t}),\mathbf{q}_{2}(\mathbf{t}),\cdots,\mathbf{q}_{m}(\mathbf{t})\right)\right\|^{2}+$
(6)
$\displaystyle\quad\quad\quad\quad+\beta_{i}\tau^{2}\left(\dot{q}_{i,0}-v_{i,0}\right)^{2}+\gamma_{i}\left(q_{i,0}-s_{i,0}\right)^{2}\Big{]}$

where $\mathcal{C}$ produces, in general, a positive, real and scalar output.
We define, in general, the notation

$\|\mathbf{g}(\mathbf{t})\|^{2}\equiv\frac{\Delta
t}{\tau}\sum_{n=0}^{N_{T}}g(t_{n})^{2}=\frac{T}{\tau
N_{T}}\sum_{n=0}^{N_{T}}g(t_{n})^{2},$ (7)

such that

$\lim_{\Delta t\to
0}\|\mathbf{g}(\mathbf{t})\|^{2}=\frac{1}{\tau}\int_{0}^{T}\mathrm{d}t\,g(t)^{2}.$
(8)

Scaling the generalised coordinates, $q_{i}$, and time, $t$, to be expressed
in terms of $\lambda_{i}$ and $\tau$, respectively, ensures a dimensionless
form; this corresponds to setting $\lambda_{i}=\tau=1$, which we do in every
example in this paper. Furthermore, we always set
$\alpha_{i}=\beta_{i}=\gamma_{i}=1$. Under these assumptions, the cost
function takes the form

$\displaystyle\mathcal{C}$
$\displaystyle=\sum_{i=1}^{m}\Big{[}\left\|\mathbf{\ddot{q}}_{i}(\mathbf{t})-\mathbf{f}_{i}\left(\mathbf{q}_{1}(\mathbf{t}),\mathbf{q}_{2}(\mathbf{t}),\cdots,\mathbf{q}_{m}(\mathbf{t})\right)\right\|^{2}+$
(9)
$\displaystyle\quad\quad\quad+\left(\dot{q}_{i,0}-v_{i,0}\right)^{2}+\left(q_{i,0}-s_{i,0}\right)^{2}\Big{]}$

The first line is the residual of the differential equation (i.e., the
differential equation should equal zero). The second line constrains
$\dot{q}_{i}(0)$ and $q_{i}(0)$ to not deviate significantly from the initial
velocity, $v_{i,0}$, and the initial position, $s_{i,0}$. The cost function,
in principle, fully specifies the initial value problem, and when the
differential equation and initial conditions are satisfied, the neural network
solution should tend towards the expected solution of the initial value
problem.

##  III Linear neural network representation of initial value problems

###  III.1 Overview

In order to set notation, we first explore free particle dynamics by
considering the simplest conceivable neural network for modelling linear
solutions to initial value problems, as shown in Fig. (1). This neural network
can be used to learn free particle dynamics and can behave as a linear
approximant for a particle in a gravitational field. We briefly explore using
the linear neural network to model the harmonic oscillator.

###  III.2 Linear neural network for a free particle

####  III.2.1 Statement of the equation of motion and cost function

Consider a system with the free particle Lagrangian

$L=\frac{1}{2}\mu\dot{q}(t)^{2}.$ (10)

Solving the Euler–Lagrange equations leads to the equation of motion

$\ddot{q}(t)=0,$ (11)

subject to the initial position $s_{0}$ and initial velocity $v_{0}$. Equation
(11) has solutions $q(t)=v_{0}t+s_{0}.$

Figure 1: A neural network to solve the differential equation (11), as
represented by Eq. (12). The network contains two neurons: i) an input neuron,
which represents the discretised temporal domain, and ii) an output neuron,
which represents the discretised solution to the initial value problem. The
network consists of two parameters: a weight between the input and output
vector and a bias associated with the output parameter.

For this dynamical system, we seek neural solutions

$\mathbf{q}(\mathbf{t})=w\mathbf{t}+b\mathbbm{1},$ (12)

where $\mathbbm{1}\in\mathbb{R}^{N_{T}}$ is a vector of ones.

Figure 1 represents the discretised solution in Eq. (12) as a neural network.
Neural networks are composed of nodes connected by edges, forming a graph
where each node is referred to as a neuron. In this paper, neurons are always
vectors — this choice means that there are fewer weights and biases to learn
and that a direct equivalency of Eq. (12) with linear regression may be
observed. Every neuron-neuron connection is associated with a weight, and
every neuron (in all layers except the input) is associated with a bias. The
free particle system is described entirely by the simplest conceivable neural
network: it contains no hidden layers or activation functions. We define the
zeroth layer neuron value as $\mathbf{t}$ and the output layer neuron as
$\mathbf{q}(\mathbf{t})$. The identification of $v_{0}$ with the weight $w$
and $s_{0}$ with the bias is trivial, i.e., $q(0)=b$, and
$\dot{q}(t)=\dot{q}(0)=w$. We can, therefore, describe the cost function as a
function of the weight and bias,

$\displaystyle\mathcal{C}(w,b)$
$\displaystyle=\|f(w\mathbf{t}-b\mathbbm{1})\|^{2}+(w-v_{0})^{2}+(b-s_{0})^{2}$
(13)
$\displaystyle=\frac{T}{N_{T}}\sum_{n=0}^{N_{T}}f(wt_{n}+b)^{2}+(w-v_{0})^{2}+(b-s_{0})^{2}.$

Given that $\ddot{q}=0\implies f(\cdot)=0$, the cost function reduces to

$\mathcal{C}(w,b)=(w-v_{0})^{2}+(b-s_{0})^{2},$ (14)

which is minimised by $w=v_{0}$ and $b=s_{0}$.

####  III.2.2 Initialisation

All parameters in any machine learning model are initialised randomly. This
initialisation has no physical interpretation; it is part of the computational
methodology, as explored in appendix C. The weight and bias for the neural
network are sampled from the uniform distribution [24]

$w,b\sim U\left(-1,1\right).$ (15)

The range of initial values of the weight and bias and their convergence to
the expected parameters, $v_{i,0}$ and $s_{i,0}$, for 250 iterations of the
training procedure is shown in Fig. (2).

Figure 2: _Free particle:_ The convergence of the weight and bias for an
initial value problem describing a free particle with $s_{0}=2.0$ and
$v_{0}=0.5$ for 250 instances of the initialisation procedure. The weight and
bias are initialised randomly according to Eq. (15) and tend towards their
expected values after sufficient iterations of an optimisation algorithm have
passed. After around 2,500 epochs, the training may be stopped as an
approximate solution has been obtained. The shaded region indicates the
maximum and minimum value of both parameters at a given epoch, and the line
indicates the average value of the parameters throughout (which is an
arithmetic mean for the underlying uniform distribution).

####  III.2.3 Learning the optimal parameters using gradient descent

Optimisation of a two parameter model can be achieved by gradient descent
(usually attributed to Cauchy). Let the model parameters be denoted by the
tuple $(w,b)$. We wish to adjust the weight and bias by some $\Delta w$ and
$\Delta b$, respectively, where the adjustment is assumed to be small to avoid
overshooting the global minima. Let the adjustments be described by the tuple
$(\Delta w,\Delta b)$. Ensuring a descent in the direction of the negative
gradient, the result of the first order expansion of the cost function about
$(\Delta w,\Delta b)$ is

$(\Delta w,\Delta b)=-\eta\left(\frac{\partial\mathcal{C}(w,b)}{\partial
w},\frac{\partial\mathcal{C}(w,b)}{\partial b}\right),$ (16)

where $\eta>0$ is a fixed global hyperparameter (a parameter which describes a
neural network) known as the step size or learning rate. Equation (16) is the
epoch-to-epoch update formula for the weight and the bias.

From the original definition of the cost function in Eq. (14),

$\frac{\partial C(w,b)}{\partial
w}=2(w-v_{0})\quad\text{and}\quad\frac{\partial C(w,b)}{\partial
b}=2(b-s_{0}).$ (17)

Equation (16) describes the gradient descent throughout the parameter space.
One set of parameter updates is associated with one epoch (one complete
forward and backwards pass through the neural network). After a sufficient
number of epochs have been completed, the weight and bias should tend toward
the dimensionless initial velocity and initial position, respectively. Figure
3 shows the gradient descent of the weight and biases throughout the parameter
space, updating according to Eq. (16).

Figure 3: The gradient descent/update procedure for the weights and biases
towards the target parameters $v_{0}$ and $s_{0}$ are shown, according to Eq.
(16). The contours represent the values of the cost function in Eq. (14).
Figure 4: _Free particle:_ The expected solution, the neural solution, and the
difference (residual) of the two are shown.

Figure 4 demonstrates the neural solution given by Eq. (12) for the learnt
parameters $w$ and $b$. Given that the values of $s_{0}$ and $v_{0}$ obtained
by the neural network are approximations, the solution they describe will
always diverge from the expected solution after sufficient time, as shown by
the residual plot.

###  III.3 Linear neural network as a linear approximant for a particle in a
gravitational field

(a)

(b)

Figure 5: (a) The linear approximant to Eq. (12) for a particle in a
gravitational field obtained by machine learning, alongside the expected
particle trajectory. Since we are trying to fit a linear function to a
quadratic function, the residuals show an underlying quadratic pattern as time
evolves. (b) The average convergence of the bias and weight for 250 iterations
of the training algorithm is again shown. The shaded region indicates the
maximum and minimum value of both parameters at a given epoch, and the line
indicates the average value of the parameters throughout (which is an
arithmetic mean for the underlying uniform distribution).

This framework additionally serves as a linear approximant for functions which
exhibit non-linear solutions. Consider a particle in a uniform gravitational
field given by the dimensionless equation of motion

$\ddot{q}(t)+1=0,$ (18)

which — in the linear neural network framework — continues to have solutions
of the form in Eq. (12). The cost function is again a function of the weight
and bias, with a non-zero forcing term

$\displaystyle\mathcal{C}(w,b)$
$\displaystyle=\frac{T}{N_{T}}\sum_{n=0}^{N_{T}}1^{2}+(w-v_{0})^{2}+(b-s_{0})^{2}$
(19) $\displaystyle=\frac{T}{N_{T}}(N_{T}+1)+(w-v_{0})^{2}+(b-s_{0})^{2}$
$\displaystyle\approx T+(w-v_{0})^{2}+(b-s_{0})^{2}$

as $N_{T}\to\infty$. The cost function is again minimised by $w$ and $b$ but
increases without bound as $T\to\infty$. Figure 5 demonstrates the linear
neural network’s ability to perform as a linear approximant for non-linear
solutions.

###  III.4 Linear neural network as a linear approximant for a harmonic
oscillator

Consider now the harmonic oscillator given by the equation of motion

$\ddot{q}(t)+q(t)=0.$ (20)

Letting $f(t)=-q(t)=-wt-b$, the cost function is

$\displaystyle\mathcal{C}(w,b)$
$\displaystyle=\frac{T}{N_{T}}\sum_{n=0}^{N_{T}}\left[w(t_{n})+b\right]^{2}+(w-v_{0})^{2}+(b-s_{0})^{2}$
(21)
$\displaystyle\approx\int_{0}^{T}\mathrm{d}t\,\\{(wt+b)^{2}\\}+(w-v_{0})^{2}+(b-s_{0})^{2}$
$\displaystyle=\frac{w^{2}}{3}(T-t_{0})^{2}+wb(T-t_{0})^{2}+b^{2}(T-t_{0})+$
$\displaystyle\quad+(w-v_{0})^{2}+(b-s_{0})^{2},$

which is minimised by

$\frac{\partial\mathcal{C}}{\partial
w}=\frac{2}{3}w(T-t_{0})^{3}+b(T-t_{0})^{2}+2(w-v_{0})=0,$ (22)

and

$\frac{\partial\mathcal{C}}{\partial
b}=w(T-t_{0})^{2}+2b(T-t_{0})+2(b-s_{0})=0.$ (23)

Solving Eqs. (22) and (23) simultaneously gives the weight and bias in terms
of the initial conditions $v_{0}$ and $s_{0}$ only,

$w=\frac{-6T^{2}s_{0}+12Tv_{0}+12v_{0}}{T^{4}+12T+4T^{3}+12},$ (24)

and

$b=\frac{4T^{3}s_{0}-6T^{2}v_{0}+12s_{0}}{T^{4}+12T+4T^{3}+12}.$ (25)

It follows that, as $T\to 0$, $b=s_{0}$ and $w=v_{0}$. Furthermore, as
$T\to\infty$, $b\to 0$ and $w\to 0$. That is to say, when probing long
dynamical regimes as $T$ grows larger, the linear approximant suggests that
the learnt weight and bias should always tend towards zero regardless of the
true values of the initial conditions. We did not observe this, and we did not
explore this further. This may suggest the need for increasing the
representational capacity of the neural network.

Figure 6 shows the values of the learnt conditions as functions of $\tau$ for
$(s_{0},v_{0})=(0.5,1.0)$.

Figure 6: The theoretical values of the weight (Eq. (24)) and bias (Eq. (25))
describing the linear approximant of the harmonic oscillator are shown. The
weight and bias are expected to tend towards zero as the time $T$ increases,
regardless of the true values of the initial conditions.

Learning linear solutions or linear approximations to solutions is both
computationally inexpensive and easy to follow mathematically. The cost
function in Eq. (14) is entirely analogous to linear regression in two
dimensions (since $\ddot{q}(t)=0$, we are only finding two optimal
parameters).

##  IV Non-linear neural network representation of initial value problems

###  IV.1 Overview

We will now explore two initial value problems with non-linear solutions: the
arbitrary angle pendulum in section IV.2 and the Hénon–Heiles system in
section IV.3. To model non-linear behaviour, we will extend the already
introduced notation to consider a deep learning framework, as shown in Fig.
(7). Deep learning increases the representational capacity of the neural
network by adding more layers (and, in turn, more network parameters) and
introduces non-linear (activation) functions. Our work on the Hénon–Heiles
system necessitates the introduction of a coupled neural network and
optimisation scheme for its solutions.

###  IV.2 Classical pendulum

####  IV.2.1 Statement of the equation of motion and cost function

Consider a one-dimensional pendulum with kinetic energy

$T(\dot{\varphi}(t))=\frac{1}{2}\mu\ell^{2}\dot{\varphi}(t)^{2},$ (26)

where $\mu$ is the mass of the pendulum bob, $\ell$ is the length of the
pendulum, and $\dot{\varphi}(t)$ is the angular velocity of the pendulum,
associated with the angle coordinate $\varphi(t)$. The potential energy is
$V(\varphi(t))=\mu g\ell[1-\cos\varphi(t)]$, where $g$ is the acceleration due
to gravity. The Lagrangian is therefore

$L(\varphi(t))=\frac{1}{2}\mu\ell^{2}\dot{\varphi}(t)^{2}-\mu
g\ell[1-\cos\varphi(t)].$ (27)

Introducing the characteristic time $\tau=\sqrt{\ell/g}$, and solving the
Euler–Lagrange equation, the dimensionless equation of motion is
$\ddot{\varphi}(t)+\sin\varphi(t)=0,$ subject to the pendulum’s initial
angular velocity $v_{0}$ and initial angle $s_{0}$.

The cost function for this problem is

$\displaystyle\mathcal{C}$
$\displaystyle=\left\|\boldsymbol{\ddot{\varphi}}(\mathbf{t})-\sin(\boldsymbol{\varphi}(\mathbf{t}))\right\|^{2}+({\dot{\varphi}}_{0}-v_{0})^{2}+({\varphi}_{0}-s_{0})^{2},$
(28)

where $\dot{\varphi}_{0}$ and $\varphi_{0}$ are the neural network’s current
values of the initial angular velocity and initial position. The cost function
is minimised by an appropriate set of weights and biases, which we will
determine soon.

####  IV.2.2 Towards deep learning

Unlike the linear neural network, there is no direct physical interpretation
of the weights and biases other than in the small region where the solution
behaves linearly. For initial value problems where their solutions exhibit
non-linear behaviour, we observe that neural networks require of the order of
$10^{4}$ or $10^{5}$ weights and biases to approximate their solution. We will
use a similar architecture to Fig. (1) where the input neuron represents the
discretised temporal domain, $\mathbf{t}$, and the output neuron represents
the discretised solution, $\boldsymbol{\varphi}(\mathbf{t})$. To increase the
representational capacity of the neural network, we construct a multi-layer,
multi-neuron network as shown in Fig. (7). Every neuron-neuron connection is
associated with a weight parameter and every neuron is associated with a bias
parameter. Appendix A provides further details of the output of a scalar-
valued neuron.

Figure 7: A fully connected and feedforward deep neural network. We annotate a
weights matrix between the $\ell$th and $(\ell-1)$th layer, assuming that both
layers contain $n$ neurons. Every connection between every neuron is
associated with a weight. We further annotate a bias vector on the $\ell$th
layer — every neuron is associated with a scalar bias. The number of layers
and neurons in each layer are hyperparameters which may need to be modified
for different dynamical systems.

The layers between the input and output layers are called hidden layers.
Computing successive hidden linear layers would be equivalent to just one
linear operation—everything collapses into a linear model so that no non-
linear problem can be learnt. To avoid this, we apply a non-linear function
$\mathcal{A}$ to every neuron; this is referred to as an activation function,
and we exploit the probabilistic interpretation of certain activation
functions (see appendix B.2) to introduce heuristically a form of model
averaging (see appendix B.1) to our solutions. For every neural network
describing an initial value problem that we explore in this section,
probabilistic activation functions are a fundamental component of the machine
learning architecture—non-probabilistic activation functions are not conducive
to learning in this context.

We assume that there are $n$ neurons in each hidden layer (in general, this
need not be the case), where each neuron has an index $j\in\\{1,n\\}$. We
further assume that the network has a total of $L$ layers (excluding the input
layer, which may be considered the $\ell=0$th layer). Then, the value of a
hidden neuron in the first layer is

$\displaystyle\mathbf{h}_{k}^{1}$
$\displaystyle=\mathcal{A}(\mathbf{z}_{k}^{1}),$ (29)
$\displaystyle\mathrm{where}\quad\mathbf{z}_{k}^{1}$
$\displaystyle=w_{k1}^{1}\mathbf{t}+b_{k}^{1}\mathbbm{1},$

and $\mathbbm{1}\in\mathbb{R}^{N_{T}}$ is a vector of ones. In the subsequent
layer(s), $\ell\in\\{2,\cdots,L-1\\}$, we similarly have

$\displaystyle\mathbf{h}_{j}^{\ell}$
$\displaystyle=\mathcal{A}(\mathbf{z}_{j}^{\ell}),$ (30)
$\displaystyle\mathrm{where}\quad\mathbf{z}_{j}^{\ell}$
$\displaystyle=\sum_{k=1}^{n}w_{jk}^{\ell}\mathbf{h}_{k}^{\ell-1}+b_{j}^{\ell}\mathbbm{1}.$

The final (output) layer, $\ell=L$, of the neural network contains a single
neuron with value

$\displaystyle\boldsymbol{\varphi}(\mathbf{t})$
$\displaystyle=\mathcal{A}(\mathbf{z}_{1}^{L}),$ (31)
$\displaystyle\mathrm{where}\quad\mathbf{z}_{1}^{L}$
$\displaystyle=\sum_{k=1}^{n}w_{1j}^{L}\mathbf{h}_{j}^{L-1}+b^{L}_{1}\mathbbm{1}.$

We find that $L=3$ layer neural networks appear to be sufficient to learn the
solutions to initial value problems, including the classical pendulum. Figure
8 shows the complete neural network representation of the classical pendulum
in a three-layer architecture (with two hidden layers and one output layer),
alongside Eqs. (29–31).

Figure 8: The neural network in Fig. (7) is annotated to show which layers
compute Eqs. (29–31). It demonstrates that the neural network is essentially a
composition of many functions. A discretised temporal domain $\mathbf{t}$ with
$N_{T}$ time steps is input to the network. The discretised solution
$\boldsymbol{\varphi}(\mathbf{t})$ is output. Information flows from bottom to
top. Every neuron is associated with a tensor of dimensions $[1,N_{T},1]$.
Every hidden layer is a tensor of shape $[1,N_{T},n]$.

####  IV.2.3 Initialisation

We initialise the weights and biases in the hidden layers from the uniform
distribution

$w_{jk}^{\ell},b_{j}^{\ell}\sim
U\left(-\sqrt{\frac{3}{n}},+\sqrt{\frac{3}{n}}\right).$ (32)

See appendix (C) for more details.

####  IV.2.4 Learning the optimal parameters using adaptive moment estimation

We wish to learn the cost function minimising mapping between the input layer
and the output layer. Gradient descent will eventually reach a global minimum
deterministically, although it may take a very long time to traverse flat
regions of the parameter space where parts of the gradient of the cost
function are approximately zero. Instead, in all examples in this section, we
use the adaptive moment estimation (Adam) algorithm, which we derive in
appendix D.4. The gradients of the weights and biases are determined by an
algorithm called _backpropagation_. The standard equations of backpropagation
for scalar-valued neurons are presented in appendix D.5. We will now derive
the equations of backpropagation for vector-valued neurons.

####  IV.2.5 Backwards propagation of errors with vector-valued neurons

Forward propagation computes an estimate of the vector
$\boldsymbol{\varphi}(\mathbf{t})$ based on the existing weights and biases.
The cost function, $\mathcal{C}$, is then calculated for that particular
estimate. We want to determine the _sensitivity_ of the cost function with
respect to the weights and biases throughout the entire network. The
sensitivities of the cost function are given by

$\displaystyle\frac{\partial\mathcal{C}}{\partial w_{jk}^{\ell}}$
$\displaystyle=\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}\cdot\nabla_{w_{jk}^{\ell}}\mathbf{z}_{j}^{\ell}$
(33a) $\displaystyle\frac{\partial\mathcal{C}}{\partial b_{j}^{\ell}}$
$\displaystyle=\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}\cdot\nabla_{b_{j}^{\ell}}\mathbf{z}_{j}^{\ell}=\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}\cdot\mathbbm{1}=\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}.$
(33b)

which we obtain by the chain rule. We use $\cdot$ to denote the scalar product
for notational simplicity. Using the definition of the output of a neuron in
the first hidden layer in Eq. (29) and in an arbitrary layer, $\ell$, in Eq.
(30), we can see that the sensitivity of the cost function with respect to the
network weights is

$\frac{\partial\mathcal{C}}{\partial
w_{jk}^{\ell}}=\begin{cases}\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}\cdot\mathbf{h}_{j}^{\ell},&\text{if
$\ell\geq 2$,}\\\
\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}\cdot\mathbf{t},&\text{if
$\ell=1$.}\end{cases}$ (34)

We seek an equation for $\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}$, which is
calculated backwards through the network for computational efficiency [25].
This is the sensitivity of the cost function with respect to the weighted
input to the neuron. In the final layer, this sensitivity is given by

$\nabla_{\mathbf{z}_{1}^{L}}\mathcal{C}=\mathop{\nabla_{\boldsymbol{\varphi}}\mathcal{C}}\mathop{\nabla_{\mathbf{z}_{1}^{L}}{\boldsymbol{\varphi}}}.$
(35)

The sensitivity is, in a linear algebra sense, a Jacobian matrix (of the
partial derivatives of the elements of $\boldsymbol{\varphi}$ with respect to
the elements of $\mathbf{z}_{1}^{L}$) multiplied by a row vector (of the
partial derivatives of $\mathcal{C}$ with respect to each element of the
solution vector $\boldsymbol{\varphi}$). Since the activation function is
applied element-wise, the Jacobian
$\nabla_{\mathbf{z}_{1}^{L}}\boldsymbol{\varphi}$ is diagonal, so the
sensitivity reduces to an element-wise multiplication, as indicated by the
Hadamard product, $\odot$. It therefore follows that

$\nabla_{\mathbf{z}_{1}^{L}}\mathcal{C}=\nabla_{\boldsymbol{\varphi}}\mathcal{C}\odot(\nabla_{\mathbf{z}_{1}^{L}}\odot{\boldsymbol{\varphi}}),$
(36)

where the substitution $\boldsymbol{\varphi}=\mathcal{A}(\mathbf{z}_{1}^{L})$
can then be made — this directly introduces the derivative of the activation
function into the backpropagation algorithm.

We then work top to bottom through the neural network (hence, backpropagation)
to determine the optimality of each neuron. Let us consider two arbitrary
layers in the network, $\ell$ and $\ell-1$, with neuron indices $j$ and $k$,
respectively. This preserves the previously defined ordering of indices on,
e.g., the weights $w_{jk}^{\ell}$. For the layer $\ell-1$, the sensitivity of
the cost function with respect to the $k$th neuron can be written in terms of
all of the sensitivities in the layer after it (this neuron influences all
neurons in the subsequent layers $\ell$, $\ell+1$, and so on, in a fully
connected neural network). The sensitivity of the output of a neuron in layer
$(\ell-1)$ is therefore also determined by the chain rule

$\nabla_{\mathbf{z}_{k}^{\ell-1}}\mathcal{C}=\sum_{j=1}^{n_{\ell}}\mathop{\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}}\mathop{\nabla_{\mathbf{z}_{k}^{\ell-1}}\mathbf{z}_{j}^{\ell}},$
(37)

where each term in the summation is a Jacobian matrix (the matrix of partial
derivatives of pre-activations from one layer to the next) multiplied by a
vector (the sensitivity in the cost function with respect to a change in the
pre-activations in the subsequent layers).

Let the Jacobian matrix
$\mathsf{J}^{jk}=\mathop{\nabla_{\mathbf{z}_{k}^{\ell-1}}\mathbf{z}_{j}^{\ell}}$
with elements

$\mathsf{J}_{mn}^{jk}=\frac{\partial(z_{j}^{\ell})_{m}}{\partial(z_{k}^{\ell-1})_{n}}.$
(38)

By the definition of the pre-activation in Eq. (30), it follows that each
element $\mathsf{J}_{mn}^{jk}$ is

$\mathsf{J}_{mn}^{jk}=\sum_{k^{\prime}=1}^{n_{\ell}}w_{jk^{\prime}}^{\ell}\frac{\partial(\mathcal{A}(({z}_{k^{\prime}}^{\ell-1})_{m}))}{\partial(({z}_{k}^{\ell-1})_{n})},$
(39)

which can only be nonzero if $k^{\prime}=k$ and $m=n$, such that the function
being differentiated is a function of the variable it is being differentiated
with respect to. Therefore, the matrices $\mathsf{J}^{jk}$ are again diagonal,
with

$\mathsf{J}_{nn}^{jk}=w_{jk}^{\ell}\frac{\partial(\mathcal{A}(({z}_{k}^{\ell-1})_{n}))}{\partial(({z}_{k}^{\ell-1})_{n})}$
(40)

and $\mathsf{J}_{m\neq n}^{jk}=0$. Substituting Eq. (39) into equation Eq.
(37) gives the error in the $(\ell-1)$th layer in terms of the errors in the
$\ell$th layer,

$\nabla_{\mathbf{z}_{k}^{\ell-1}}\mathcal{C}=\sum_{j=1}^{n_{\ell}}w_{jk}^{\ell}\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}\odot\left(\nabla_{\mathbf{z}_{k}^{\ell-1}}\odot\mathcal{A}(\mathbf{z}_{k}^{\ell-1})\right).$
(41)

It immediately follows that between the output and the penultimate layer,

$\nabla_{\mathbf{z}_{j}^{\ell}}\mathcal{C}=w_{1j}^{L}\nabla_{\mathbf{z}_{1}^{L}}\mathcal{C}\odot\left(\nabla_{\mathbf{z}_{j}^{\ell}}\odot\mathcal{A}(\mathbf{z}_{j}^{\ell})\right).$
(42)

Figure 9: The phase portrait for four open and four closed trajectories. Small
artefacts can be observed (for example, on the closed trajectory with the
highest initial angle) since the error in the neural solution increases as
time increases. Note that the unstable trajectory at $\varphi(0)=\pi$ and
$\dot{\varphi}(0)=0$ does not appear on the phase space portrait — it only
exists at a single point. Figure 10: Cost metrics for all trajectories in the
classical pendulum that were probed in Fig. (9). The closed trajectories have
initial angles $\varphi_{1,0}=\pi/8$, $\varphi_{2,0}=\pi/3$,
$\varphi_{3,0}=13\pi/24$ and $\varphi_{4,0}=3\pi/4$ (these are four equally
spaced angles covering a broad range of dynamical behaviour) and initial
angular velocities $\dot{\varphi}_{i=1,\cdots,4}=0$. The open, positive
trajectories have initial angles $\varphi_{5,0}=-\pi$ and $\varphi_{6,0}=-\pi$
and initial angular velocities $\dot{\varphi}_{5,0}=\pi/2$ and
$\dot{\varphi}_{6,0}=\pi$. The open, negative trajectories have initial angles
$\varphi_{7,0}=\pi$ and $\varphi_{8,0}=\pi$ and initial angular velocities
$\dot{\varphi}_{7,0}=-\pi/2$ and $\dot{\varphi}_{8,0}=-\pi$. All cost metrics
use a logarithmic scale.

####  IV.2.6 Results and discussion

The exact solution for the pendulum dynamics, for arbitrary initial angle and
angular velocity, may be written in terms of Jacobi elliptic functions, but
the standard method of solution is via numerical integration schemes. There is
a numerical instability when $\varphi(0)=\pi$ and $\dot{\varphi}(0)=0$; the
pendulum should remain at $\varphi(t)=\pi$ for all time, but machine precision
means that the pendulum will always fall, as the pendulum is never at exactly
$\varphi(0)=\pi$.

We sample both open and closed trajectories of the classical pendulum. The
closed trajectories have initial positions from $s_{0}=\pi/8$ to
$s_{0}=3\pi/4$ — equally spaced and chosen to cover a broad range of dynamical
behaviour — and zero velocity, $v_{0}=0$. The open trajectories have initial
positions $s_{0}=\pm\pi$ and velocities $v_{0}=\pm\pi/2$ or $v_{0}=\pm\pi$
(the open and positive trajectories have negative initial angle and positive
initial angular velocity, and the open and negative trajectories have positive
initial angle and negative initial angular velocity). All trajectories are
shown in Fig. 9. Figure 10 shows the various cost metrics for each trajectory.
The machine learning modelling estimates the solution of the initial value
problem for a given set of weights and biases, and determines the mean square
error of the trajectory across the problem domain. In order to plot a quantity
with physical dimensions, we take the square root of this estimate; we define
the quantity
$\sqrt{\left\|\mathbf{\ddot{q}}_{i}(\mathbf{t})-\mathbf{f}_{i}\left(\mathbf{q}_{1}(\mathbf{t}),\mathbf{q}_{2}(\mathbf{t}),\cdots,\mathbf{q}_{m}(\mathbf{t})\right)\right\|^{2}}$
to be the root-mean-square of the differential equation (RMS of the ODE).
Remarkably, we note that the machine learning modelling successfully captured
the numerical instability at $\pi$ radians, although this took around ten
times more epochs than other trajectories.

###  IV.3 Hénon–Heiles system

####  IV.3.1 Statement of the equation of motion and cost function

The Hénon–Heiles system [22] describes the non-linear motion of a star around
its galactic centre. The motion can be viewed as a pair of harmonic
oscillators with a perturbation (associated with a galactic potential) which
couples the oscillators together. The motion is chaotic and non-dissipative.
This system can, again, be described in terms of a Lagrangian of the form
$L=T-V$, where the kinetic energy is

$T(\dot{x}(t),\dot{y}(t))=\frac{1}{2}m\left[\dot{x}^{2}(t)+\dot{y}^{2}(t)\right],$
(43)

and the potential is

$V(x(t),y(t))=\frac{1}{2}m\omega^{2}\left[x^{2}(t)+y^{2}(t)\right]+k\left[x^{2}(t)y(t)-\frac{y^{3}(t)}{3}\right].$
(44)

Ensuring a dimensionless system of equations by setting $m=\omega=k=1$, and by
solving the Euler–Lagrange equations, we obtain the system of coupled
equations of motion for each coordinate

$\displaystyle\ddot{x}(t)+x(t)-2x(t)y(t)$ $\displaystyle=0,$ (45)
$\displaystyle\ddot{y}(t)+y(t)+y^{2}(t)-x^{2}(t)$ $\displaystyle=0.$

The system is time-independent, and so the total energy $E=T+V$ is a conserved
quantity.

The equations of motion and initial conditions lead to the cost function

$\displaystyle\mathcal{C}$
$\displaystyle=\|\mathbf{\ddot{x}}(\mathbf{t})+\mathbf{x}(\mathbf{t})-2\mathbf{x}(\mathbf{t})\mathbf{y}(\mathbf{t})\|^{2}$
(46)
$\displaystyle\quad+\|\mathbf{\ddot{y}}(\mathbf{t})+\mathbf{y}(\mathbf{t})+\mathbf{y}^{2}(\mathbf{t})-\mathbf{x}^{2}(\mathbf{t})\|^{2}$
$\displaystyle\quad+(x_{0}-s_{x,0})^{2}+(y_{0}-s_{y,0})^{2}$
$\displaystyle\quad+(\dot{x}_{0}-v_{x,0})^{2}+(\dot{y}_{0}-v_{y,0})^{2}.$

####  IV.3.2 Coupled neural networks

We have two dynamical degrees of freedom, which we handle by employing two
neural networks — one for each coordinate — with separate parameterisations
that we simultaneously optimise. The cost function is shared between the two
networks, although note that there are two separate backwards propagations of
errors needed to update the parameters of each network. The gradients computed
during both backward passes are based on the cost function in Eq. (46), which
includes terms for both coordinates. Because these coordinates are coupled in
the equations of motion, the gradient of the cost function with respect to one
set of parameters will naturally include contributions from both coordinates.
The presence of coupling in the cost function means that we must ensure that
the updates to the parameters for one network reflect the influence of the
other network. The first backward propagation computes the gradients in the
network representing $x$; the gradient information is retained after the
errors associated with the first coordinate are computed. The second backward
propagation computes the gradients in the network representing $y$; all
gradient information is then discarded since all backwards passes through all
coordinates are complete.

Once we calculate the gradients, the optimiser takes a step to update the
parameters in the direction that minimises the cost function. The process of
zeroing (or resetting) the gradients, calculating new gradients, and then
stepping the optimiser is repeated many times in a training loop until the
networks converge to a solution that minimises the cost function. In
principle, the coupled optimiser could include as many optimisers as computer
memory allows. We refer the reader to our code [26] for further implementation
details.

Figure 11 shows a schematic for the coupled neural network for the
Hénon–Heiles system.

Figure 11: A coupled neural network for the Hénon–Heiles system. It propagates
independent parameters for each generalised coordinate (given a common
temporal domain) and couples only in the evaluation of the cost functions.

####  IV.3.3 Results and discussion

(a)  (b)

Figure 12: Two Hénon–Heiles trajectories obtained by machine learning with
different energies: (a) corresponds to $E=1/6$ (the chaotic trajectory) and
(b) corresponds to $E=1/12$ (the quasi-periodic trajectory). We show
comparative metrics for the same trajectories obtained by Runge–Kutta methods
in Fig. (13) and Fig. (14). Top: the time evolution of the individual
coordinates. Bottom: phase profiles, showing $(x,y)$ pairs as parametric
functions of $t$. Figure 13: _Hénon–Heiles system for energy $E=1/6$_. Row 1:
the evolution of the action as training progresses. The dashed, grey
horizontal lines indicate the action of the system obtained by using a
Runge–Kutta 4th-order integration scheme. The dashed, yellow trajectories are
associated with $N_{T}=512$ timesteps. The solid, orange trajectories are
associated with $N_{T}=256$ timesteps. As the training evolves, an extremum of
the action is also found, taking similar values to the Runge–Kutta integration
schemes. Row 2: the root mean square of the residual of the differential
equations for the respective coordinates. Row 3: the residuals in the initial
positions for both coordinates. Row 4: the residuals in the initial velocities
for both coordinates. The initial conditions are easier to match since the
solution begins to diverge as the time coordinate increases. Figure 14:
_Hénon–Heiles system for energy $E=1/12$_. Row 1: the evolution of the action
as training progresses. The dashed, grey horizontal lines indicate the action
of the system obtained by using a Runge–Kutta 4th-order integration scheme.
The dashed, pink trajectories are associated with $N_{T}=512$ timesteps. The
solid, purple trajectories are associated with $N_{T}=256$ timesteps. As the
training evolves, an extremum of the action is also found, taking similar
values to the Runge–Kutta integration schemes. Row 2: the root mean square of
the residual of the differential equations for the respective coordinates. Row
3: the residuals in the initial positions for both coordinates. Row 4: the
residuals in the initial velocities for both coordinates. The initial
conditions are easier to match since the solution begins to diverge as the
time coordinate increases.

We will only consider confined (closed) trajectories. Hénon and Heiles give an
upper limit for the energy to produce a confined trajectory as $E=1/6$. This
trajectory is entirely chaotic. We also study the confined trajectory with
energy $E=1/12$, which is quasi-periodic. We choose always to set the initial
position of the first coordinate such that $s_{x,0}=0$. For the chaotic
trajectory, where $E=1/6$, we set the initial velocity for the $x$-coordinate
to $v_{x,0}=0.10$, and the initial position and velocity for the
$y$-coordinate to $s_{y,0}=0.57$ and $v_{y,0}=0.057$. For the quasi-periodic
trajectory, where $E=1/12$, we set the velocity for the $x$-coordinate to
$v_{x,0}=0.05$, and the initial position and velocity for the $y$-coordinate
to $s_{y,0}=0.40$ and $v_{y,0}=0.041$.111The initial velocity of the first
coordinate is constrained by
$\dot{x}(0)=\sqrt{2E-v_{y,0}^{2}-s_{y,0}^{2}+2s_{y,0}^{3}/3}$, where we have
set either $E=1/6$ or $E=1/12$. Given these constraints, the initial position
and initial velocity for the second coordinate can be arbitrarily chosen,
providing that $\dot{x}(0)$ remains real-valued. To ensure an appropriate
choice of the initial conditions for the second coordinate, we find $y(0)$ by
using root finding algorithms on the constraint $V<E$, and setting the result
as 10% of this value. We find $\dot{y}(0)$ by evaluating $T<E-V$, and setting
the result as 10% of this value.

Figures 12a and 12b show these two Hénon–Heiles trajectories. We used the
SechLU activation function, given in Eq. (53b). GELU, given in Eq. (52b), was
not conducive to learning and all trajectories failed to obtain correct
dynamics; we are unsure why this is the case, but it is our empirical
observation that SechLU is the most reliable of the probabilistic activation
functions that we have tested. Figures 13 and 14 show the individual cost
metrics for a Hénon–Heiles system with energies $E=1/6$ and $E=1/12$,
respectively.

The neural solutions successfully minimise the action, $S$, as training
progresses (we stress that the action is not an explicit term in the cost
function, although it generates the dynamical equations through Hamilton’s
principle). We approximate the action at every epoch using $S\approx\Delta
t\sum_{k=0}^{N_{T}}L(x_{k},\dot{x}_{k},y_{k},\dot{y}_{k})$. We note that the
action of the machine learning trajectories does not necessarily increase or
decrease monotonically with epoch. We note that one might expect that forming
a cost function around the action only and then minimising it would be
sufficient, as an almost literal realisation of Hamilton’s principle of
stationary action (the fact that the action must be extremised rather than
strictly minimised notwithstanding). We have not found this to be the case,
neither in conjunction with a trajectory’s initial values, nor with coordinate
values at the path’s two endpoints. The Euler–Lagrange equations, i.e., the
dynamical consequence of Hamilton’s principle, in conjunction with the initial
conditions, appear to be necessary, which also have the useful property of the
minimising value of the cost function always being equal to zero. We note that
when the root mean square of the differential equation remains constant over
many epochs, it corresponds to unchanging values of the action. This is due to
the equivalence of the principle of stationary action and the underlying
equations of motion.

The trajectories obtained by the Runge–Kutta integration routines may not give
the ‘true’ value of the action (especially in systems with chaotic dynamics).
However, since multiple methods are approaching similar values of the action,
this value could be interpreted as the most stationary.

We finally note that a finer temporal discretisation appears to improve the
learning capabilities of the neural network, especially with the chaotic
trajectory, without an increase in network parameters (since the neurons are
vector-valued). The stability of probing chaotic dynamics using numerical
integration schemes such as the Runge–Kutta methods can be improved (up to a
point) by using a finer temporal resolution, which is what we observe in our
neural solutions.

##  V Conclusions

We have presented the neural initial value problem: a machine learning
framework which approximates the solutions to a range of physics-based initial
value problems in the absence of training data (i.e., an unsupervised learning
approach). We have formulated a general cost function that captures the
dynamics of a wide range of mechanical systems. Through extensive
experimentation, we demonstrated that our approach can successfully model
systems with non-linear, coupled, and chaotic behaviours. For the free
particle and particle in a gravitational field, our linear neural network
provided accurate solutions and served as a linear approximant for problems
with non-linear solutions. In the case of the pendulum (for arbitrary angles,
i.e., outside the small angle approximation), the deep neural network
framework proved effective in capturing non-linear dynamics. For the
Hénon–Heiles system, our coupled neural network was successful in capturing
the dynamics of the two coordinates, including in chaotic regimes. We have
also introduced and evaluated several machine learning techniques to enhance
our framework: probabilistic activation functions, which appear to be
necessary to learn the solutions of initial value problems (they may be
interpreted as a general class of model averaging functions), and the coupled
optimisation routines allow for the simultaneous optimisation of many neural
networks representing individual generalised coordinates.

Whilst not a direct replacement for standard numerical integration routines
such as the Runge–Kutta methods, we have demonstrated that neural initial
value problems are a complementary approach which may be used to verify the
validity of other numerical methods, for example in highly chaotic regimes.
Whilst not presented in the current work, we have also found that coupled
neural networks — alongside the neural initial value problem framework — can
solve complex-valued initial value problems by treating the underlying
differential equations as a system of ordinary differential equations for the
real and imaginary parts (i.e., considering the complex domain to be an
ordered pair of real numbers). The methods presented in this paper could be
augmented with the presence of training data, where the probabilistic
activation functions could provide greater accuracy and/or generalisation
capabilities for a broader range of dynamical systems.

Future work may also be directed towards understanding the non-monotonicity of
the action in the optimisation procedure: for example whether the neural
network settles temporarily on a different solution (which one might consider
to be an ODE equivalent of ‘hallucination’ phenomena observed in AI solutions
more generally) obeying Hamilton’s principle, possibly one which satisfies the
contemporaneous approximation for the initial conditions. Another question of
principle is whether a strict action minimisation procedure (choosing two
endpoints and then trying to find the minimal action between the points) can
be implemented using the framework in this paper.

##  VI Acknowledgements

We thank the EPSRC, grant EP/T015241/1, for funding this project.

##  VII Code availability

The code for this project is available on GitHub [26].

##  Appendix A Output of a scalar-valued neuron

In the main text we consider vector neurons exclusively. Here we briefly
present the notation for scalar neurons, with the intention of making the
parallels between scalar and vector neurons readily apparent.

Let the $k$th neuron from the $(\ell-1)$th layer have an output value of
$y_{k}^{\ell-1}$ and let the $j$th neuron from the $\ell$th layer have an
output value of $y_{j}^{\ell}$. Let $w_{jk}^{\ell}$ be the weight from the
$k$th neuron in the $(\ell-1)$th layer to the $j$th neuron in the $\ell$th
layer. We assume that the neural network is fully connected (all neurons in
one layer are connected to all neurons in the subsequent layer) and
feedforward (there are no loops or recurrent connections). Figure 15 shows
this notation between two layers of a neural network.

Figure 15: Two layers of a feedforward neural network are shown. A connection
between the first neuron of the fourth layer and the third neuron of the fifth
layer is highlighted. This connection has a weight $w_{jk}^{\ell}=w_{31}^{5}$.

Using this notation, the output, $y_{j}^{\ell}$, of the $j$th neuron in the
$\ell$th layer then evaluates as

$\displaystyle z_{j}^{\ell}$
$\displaystyle=\sum_{k=1}^{n_{\ell-1}}w_{jk}^{\ell}y_{k}^{\ell-1}+b_{j}^{\ell},$
(47a) $\displaystyle y_{j}^{\ell}$
$\displaystyle=\mathcal{A}\left(z_{j}^{\ell}\right),$ (47b)

where Eq. (47a) is akin to a linear regression and Eq. (47b) is referred to as
the activation or output of a neuron in the deep neural network.

Activation functions are introduced between layers to invoke the _universal
approximation theorem_ [28, 29, 30] — a statement that, with the addition of
nonlinear activation functions, a possibly infinitely wide neural network
trained for possibly infinite time can represent any continuous function222The
universal approximation theorems do not tell us _how_ to find the optimal
weights and biases, only that they exist.. In general, it is not obvious _a
priori_ which activation function to choose; an activation function suitable
for one problem is not necessarily suitable for a different problem, as a
consequence of the _No Free Lunch_ theorem [32].

##  Appendix B Probabilistic activation functions

###  B.1 Model averaging and dropout

Neural networks exhibit run-to-run variation — different initialisations of
the network parameters and stochastic descent through the (possibly vast)
parameter space mean that each neural network will return a different output
compared to another instance of a neural network with the same architecture,
cost function, inputs and hyperparameters. Different instances of the neural
network explore different local minima, which in turn produce different sets
of optimal parameters. It may be that some runs fail to learn the mapping
between input and output, whilst others learn a more accurate mapping. It has
been noted extensively [33, 34, 35, 36, 37] that performing an average
(arithmetic or geometric mean) of many neural networks produces a more
accurate model than any individual neural network.

Model averaging may be heuristically achieved using an algorithm called
dropout. Consider the output of the $j$th neuron in the $\ell$th layer, given
in Eqs. (47a) and (47b). Dropout, introduced by Srivastava _et al._ [38],
removes neurons with probability $p$ and preserves neurons with probability
$1-p$ at the start of each epoch. We can represent the effect of dropout
through a set of Bernoulli trial outcomes $r_{k}^{\ell-1}$ ($=0$ with
probability $p$ and otherwise $=1$). We then define
$\tilde{y}_{k}^{\ell-1}=r_{k}^{\ell-1}y_{k}^{\ell-1}$, such that, instead of
Eq. (47a), we have

$z_{j}^{\ell}=\sum_{k=1}^{n_{\ell-1}}w_{jk}^{\ell}\tilde{y}_{k}^{\ell-1}+b_{j}^{\ell},$
(48)

with $y_{j}^{\ell}=\mathcal{A}(z_{j}^{\ell})$, as in Eq. (47b). Figure 16
shows a feedforward neural network, with and without dropout. Configurations
with dropout are called subnetworks, as they are a subset of the larger
network. Neurons in all layers except the output layer can be subject
todropout. It is in principle possible that there are no connections between
the input and output layer, but this becomes vanishingly unlikely for larger
networks.

(a)  (b)

Figure 16: (a) A fully connected, feedforward neural network without dropped
neurons and (b) the same neural network with dropped neurons (dark grey) and
the remaining connections. The neurons which are dropped are stochastically
chosen at each epoch.

Dropout prevents overfitting (referred to as regularisation), and also
approximately averages the output of multiple neural networks. Dropout is
carried out at the start of each epoch. For a network with $N$ neurons in the
hidden layers, there are $2^{N}$ possible configurations of the neurons in the
network. Dropout is thus akin to an ensemble of $2^{N}$ instances of a neural
network, and can be described as an average of many models produced by the
subnetworks.

A process somewhat analogous to dropout can be carried out via a class of
activation functions we call probabilistic activation functions, extending
work by Hendrycks and Gimpel [39]. It is our observation that probabilistic
activation functions are required for learning solutions of initial value
problems.

###  B.2 Connection to probabilistic activation functions

We define the _probabilistic activation function_

$\mathcal{A}(x)=x\phi(x)=x\int_{-\infty}^{x}\mathrm{d}x^{\prime}\,p(x^{\prime}),$
(49)

where $p(x)$ is a probability distribution supported on the domain
$x\in(-\infty,\infty)$ and the integral defines its cumulative distribution
function $\phi(x)$ 333Probability distributions with support only on the
domain $x\in(0,\infty)$ are not conducive to learning, possibly due to the
lack of symmetry. It is our empirical observation that some distributions
which are not perfectly symmetric about the origin — such as that used in Mish
(see Eq. (56b)) — are still conducive to learning..

Although not smooth or differentiable, it is insightful to consider the Dirac
delta distribution, $p(x)=\delta(x)$. The associated cumulative distribution
function is the Heaviside step function,

$\phi(x)=\int^{x}_{-\infty}\mathrm{d}x^{\prime}\,\delta(x^{\prime})=\Theta(x),$
(50)

which leads to the standard (and also non-differentiable) ReLU activation
function

$\mathrm{ReLU}(x)=x\mathop{\Theta(x)}=\max(0,x).$ (51)

The effect is therefore comparable to a form of dropout; if $x$ is viewed as a
random variable, then neurons are either preserved or dropped, with
probability $0.5$. More generally, the cumulative distribution function
$\phi(x)$ for a continuous, differentiable distribution function $p(x)$ will
interpolate smoothly between being close to $0$ over much of the initial part
of its range, and close to $1$ over the latter part of its range. Note also
that if $\phi(x)$ is the probability for the outcome of a Bernoulli trial to
be $=1$, it is also the expectation value for the Bernoulli trial. Hence, a
probabilistic activation function $\mathcal{A}(x)=x\mathop{\phi(x)}$ scales
all neuron outputs by the expectation value of a Bernoulli trial, effectively
providing something resembling a smooth implementation of dropout and a form
of model averaging.

In this way one can derive the Gaussian error linear unit (GELU) activation
function [39]

$\displaystyle p(x)$
$\displaystyle=\frac{1}{\rho\sqrt{2\pi}}\mathrm{exp}\left(-\frac{x^{2}}{2\rho^{2}}\right),$
(52a) $\displaystyle\mathcal{A}(x)$
$\displaystyle=\mathrm{GELU}(x)=\frac{x}{2}\left[1+\mathrm{erf}\,\left(\frac{x}{\rho\sqrt{2}}\right)\right].$
(52b)

An activation function we have found in this work to be very useful we call
SechLU:

$\displaystyle p(x)$
$\displaystyle=\frac{1}{2\rho}\mathrm{sech}^{2}\left(\frac{x}{\rho}\right),$
(53a) $\displaystyle\mathcal{A}(x)$
$\displaystyle=\mathrm{SechLU}(x)=\frac{x}{1+\mathrm{exp}(-2x/\rho)}.$ (53b)

Starting with a Cauchy (Lorentzian) distribution, we can similarly obtain an
activation function we call CauchyLU:

$\displaystyle p(x)$
$\displaystyle=\frac{1}{\pi\rho}\left(1+\frac{x^{2}}{\rho^{2}}\right)^{-1},$
(54a) $\displaystyle\mathcal{A}(x)$
$\displaystyle=\mathrm{CauchyLU}(x)=\frac{x}{2}\left[1+\frac{2}{\pi}\arctan\left(\frac{x}{\rho}\right)\right].$
(54b)

Noting it is not everywhere differentiable, we can obtain an activation
function which we call LaplaceLU:

$\displaystyle p(x)$
$\displaystyle=\frac{1}{2\rho}\mathrm{exp}\left(-\frac{|x|}{\rho}\right)$
(55a) $\displaystyle\mathcal{A}(x)$
$\displaystyle=\mathrm{LaplaceLU}(x)=\frac{x}{2}\left\\{1+\operatorname{sgn}(x)\left[1-\exp\left(-\frac{|x|}{\rho}\right)\right]\right\\}.$
(55b)

Finally, the Mish activation function, introduced by Misra in 2019 [41], can
also be derived from a (somewhat unusual) probability distribution.

$\displaystyle p(x)$
$\displaystyle=\frac{1}{\rho}\mathop{\mathrm{sech}^{2}\left(\ln\left|1+\mathrm{e}^{x/\rho}\right|\right)}\frac{\mathrm{e}^{x/\rho}}{1+\mathrm{e}^{x/\rho}},$
(56a) $\displaystyle\mathcal{A}(x)$
$\displaystyle=\mathrm{Mish}(x)=x\mathop{\tanh\left(\ln\left|1+\mathrm{e}^{x/\rho}\right|\right)}.$
(56b)

The probability distributions, cumulative distributions, activation functions
and their derivatives are shown in Fig. 17.

In all cases, the width parameter $\rho$ may be learnt as an additional
hyperparameter by following the same initialisation and optimisation routines
as for the weights and biases. We observe faster and more reliable convergence
to appropriate physical solutions if we treat $\rho$ as a free parameter.

##  Appendix C Weight and bias initialisation

The weights and the biases in the neural network must be initialised before
any training or learning can take place, where it is important to note that
this is not a “best guess” at the solution. The values for this parameter
initialisation should be neither too large (to avoid diverging values of the
gradients) nor too small (to avoid slow training).

In what follows we outline an initialisation scheme general to any
differentiable activation function based on that described by Glorot and
Bengio [42], which assumed the consistent use throughout the network of the
hyperbolic tangent. Without loss of generality we may ignore the bias terms to
simplify the mathematics (the biases are, in any case, drawn from the same
distribution as the weights). We assume the weights distribution has mean $=0$
and variance $=\sigma^{2}$, and that the inputs $y_{k}^{\ell-1}$ are
distributed with mean $=0$ and variance $=\rho^{2}$ (in general
$\sigma^{2}\neq\rho^{2}$). The parameters are initialised with zero mean to
avoid erratic and undesirable zig-zagging through the parameter space during
optimisation, which may cause overshooting of the minima. The requirement for
constant variance, i.e.,
$\mathrm{Var}(y_{k}^{\ell-1})\approx\mathrm{Var}(y_{j}^{\ell})$, prevents the
layers’ output values from exponentially vanishing or exploding as one moves
from bottom to top through the network, which would make the network difficult
or impossible to train.

Early on in the training, the values of $z_{j}^{\ell}$ can be considered to be
small [42], such that the action of a differentiable activation function can
be approximated by
$y_{j}^{\ell}=\mathcal{A}(z_{j}^{\ell})\approx\mathcal{A}(0)+\mathcal{A}^{\prime}(0)z_{j}^{\ell}$.
We can therefore consistently approximate the expectation value by

$\mathbb{E}(y_{j}^{\ell})\approx\mathcal{A}(0)+\mathop{\mathcal{A}^{\prime}(0)}\mathop{\mathbb{E}(z_{j}^{\ell})},$
(57)

and the variance
$\mathrm{Var}(y_{j}^{\ell})=\mathbb{E}[(y_{j}^{\ell})^{2}]-[\mathbb{E}(y_{j}^{\ell})]^{2}$
by

$\mathrm{Var}(y_{j}^{\ell})\approx\mathcal{A}^{\prime}(0)^{2}\mathop{\mathrm{Var}}(z_{j}^{\ell}).$
(58)

Hence, assuming $w_{jk}^{\ell}$ and $y_{k}^{\ell-1}$ to be independent random
variables,
$\mathbb{E}(w_{jk}^{\ell}y_{k}^{\ell-1})=\mathop{\mathbb{E}}(w_{jk}^{\ell})\mathop{\mathbb{E}(y_{k}^{\ell-1})}$,
and

$\mathbb{E}(y_{j}^{\ell})\approx\mathcal{A}(0)+\mathcal{A}^{\prime}(0)\sum_{k=1}^{n_{\ell-1}}\mathop{\mathbb{E}}(w_{jk}^{\ell})\mathop{\mathbb{E}(y_{k}^{\ell-1})}.$
(59)

As we have already assumed the distribution from which the weights are drawn,
and the distribution of the neuron inputs, both have zero mean, it follows
that $\mathbb{E}(y_{j}^{\ell})\approx\mathcal{A}(0)$. Similarly,

$\begin{split}\mathrm{Var}(y_{j}^{\ell})\approx&\mathcal{A}^{\prime}(0)^{2}\sum_{k,k^{\prime}=1}^{n_{\ell-1}}\mathbb{E}(w_{jk}^{\ell}y_{k}^{\ell-1}w_{jk^{\prime}}^{\ell}y_{k^{\prime}}^{\ell-1})\\\
=&\mathcal{A}^{\prime}(0)^{2}\sum_{k=1}^{n_{\ell-1}}\mathop{\mathbb{E}[(w_{jk}^{\ell})^{2}]}\mathop{\mathbb{E}[(y_{k}^{\ell-1})^{2}]},\end{split}$
(60)

as
$\mathbb{E}(w_{jk}^{\ell}y_{k}^{\ell-1}w_{jk^{\prime}}^{\ell}y_{k^{\prime}}^{\ell-1})=\mathop{\mathbb{E}(w_{jk}^{\ell})}\mathop{\mathbb{E}(y_{k}^{\ell-1})}\mathop{\mathbb{E}(w_{jk^{\prime}}^{\ell})}\mathop{\mathbb{E}(y_{k^{\prime}}^{\ell-1})}=0$
if $k\neq k^{\prime}$ and the variables are therefore all independent. Hence,
$\mathrm{Var}(y_{j}^{\ell})\approx\mathcal{A}^{\prime}(0)^{2}n_{\ell-1}\sigma^{2}\rho^{2}.$

For the variance to be the same from layer to layer, i.e.,
$\mathrm{Var}(y_{j}^{\ell})=\mathrm{Var}(y_{k}^{\ell-1})=\rho^{2}$, it follows
that
$\mathcal{A}^{\prime}(0)^{2}n_{\ell-1}\sigma^{2}=\mathcal{A}^{\prime}(0)^{2}n_{\ell}\sigma^{2}=1$
(noting that $n_{\ell-1}\neq n_{\ell}$ in general). We therefore take the
average $\mathcal{A}^{\prime}(0)^{2}(n_{\ell-1}+n_{\ell})\sigma^{2}/2=1$, from
which it follows that

$\sigma^{2}=\frac{2}{\mathcal{A}^{\prime}(0)^{2}(n_{\ell-1}+n_{\ell})}.$ (61)

Generalised Xavier initialisation uses a uniform (rectangular) distribution
over the interval $[-a,a]$, which has mean $=0$ and variance
$\sigma^{2}=a^{2}/3$. Equating this with Eq. (61) gives the values of the
support of the distribution in terms of the number of neurons in each layer

$w_{jk}^{\ell}\sim
U\left(-\sqrt{\frac{6}{\mathcal{A}^{\prime}(0)^{2}(n_{\ell-1}+n_{\ell})}},\sqrt{\frac{6}{\mathcal{A}^{\prime}(0)^{2}(n_{\ell-1}+n_{\ell})}}\right).$
(62)

For $\tanh(x)$, the scaling factor $\mathcal{A}^{\prime}(0)=1$. For any
probabilistic activation function, the scaling factor
$\mathcal{A}^{\prime}(0)=\phi(0)$. If the underlying probability distribution
is symmetric, then ${\mathcal{A}^{\prime}(0)=1/2}$.

Note that the initialisation we employ for linear neural networks in section
III.2.2 is the Kaiming He scheme, where one samples from an unscaled uniform
distribution [24]:

$w_{jk}^{\ell}\sim
U\left(-\sqrt{\frac{1}{n_{\ell}}},\sqrt{\frac{1}{n_{\ell}}}\right).$ (63)

If $n_{\ell}=1$, then this reduces to Eq. (15) in the main text.

##  Appendix D Learning algorithms

###  D.1 Statement of the problem

The machine learning problem is to find the set of weights and biases which
approximate our solution, and we use our cost function to check if we are
converging to the expected solution. In all examples in this paper, the
weights and biases are adjusted by the adaptive moment estimation or _Adam_
optimisation algorithm [43]. Adam is the descendent of other optimisation
algorithms which have emerged over the last few decades. We will follow Adam’s
lineage, starting with gradient descent, before arriving at the equations for
its implementation.

###  D.2 Gradient descent

One can extend the gradient descent algorithm introduced in section III.2.3 to
consider any number of weights and biases. All cost functions are
parameterised by their weights and biases $\boldsymbol{\theta}$. Let us again
ignore the possibility of biases in the network and focus on the adjustments
required for the weights. For an $L$-layer neural network,

$\mathcal{C}(\boldsymbol{\theta})\equiv\mathcal{C}(w_{11}^{1},\cdots,w_{n_{\ell}n_{\ell-1}}^{\ell},\cdots,w_{n_{L}n_{L-1}}^{L}).$
(64)

We wish to adjust the weights by some amount $\Delta w_{jk}^{\ell}$ (and, in
general, the biases by some assumed small amount $\Delta b_{j}^{\ell}$). These
adjustments are described collectively by the vector
$\Delta\boldsymbol{\theta}$. Taylor expanding,

$\mathcal{C}(\boldsymbol{\theta}+\Delta\boldsymbol{\theta})\approx\mathcal{C}+(\Delta\boldsymbol{\theta})\cdot\nabla_{\boldsymbol{\theta}}\mathcal{C}(\boldsymbol{\theta}),$
(65)

where by definition

$\nabla_{\boldsymbol{\theta}}\mathcal{C}(\boldsymbol{\theta})=\left(\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{11}^{1}},\cdots,\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{n_{\ell}n_{\ell-1}}^{\ell}},\cdots,\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{n_{L}n_{L-1}}^{L}}\right).$ (66)

For $\mathcal{C}\to 0$, we require in general that
$\mathcal{C}(\boldsymbol{\theta}+\Delta\boldsymbol{\theta})<\mathcal{C}(\boldsymbol{\theta})$,
and that the adjustments point in the opposite direction to the gradient,
i.e.,

$\Delta\boldsymbol{\theta}=-\eta\nabla_{\boldsymbol{\theta}}\mathcal{C}(\boldsymbol{\theta}),$
(67)

where $\eta>0$ is (so far) a fixed global hyperparameter (a parameter which
describes a neural network) known as the step size or learning rate. One,
therefore, has the epoch-to-epoch update formulas for the weights and, by
corollary, for the biases

$\Delta
w_{jk}^{\ell}=-\eta\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\quad\text{and}\quad\Delta
b_{j}^{\ell}=-\eta\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
b_{j}^{\ell}}.$ (68)

###  D.3 Adaptive gradient methods (AdaGrad and RMSProp)

The cost function may be highly sensitive to certain directions in the
parameter space $\boldsymbol{\theta}$ and insensitive to others. A class of
learning algorithms (AdaGrad [44], RMSProp [45] and Adam [43]) instead
introduce adjustable hyperparameters $\eta_{ij}^{\ell}$ for every parameter in
the neural network to adjust the rate of descent across the parameter space.

AdaGrad adapts the learning rates for each parameter by scaling them by a
rolling history of previous gradients $r_{ij}^{\ell}$,

$\displaystyle r_{jk}^{\ell}$
$\displaystyle=r_{jk}^{\ell}+\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right]^{*}\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ $\displaystyle\Delta w_{jk}^{\ell}$
$\displaystyle=-\eta\frac{1}{\delta+\sqrt{r_{jk}^{\ell}}}\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ (69)

where ∗ denotes the complex conjugate in the case of complex-valued neural
networks, $\eta$ is the _global_ learning rate and $\delta\sim 10^{-8}$ is a
small parameter for numerical stability. The initial values of $r_{jk}^{\ell}$
are all zero.

RMSProp introduces a decay rate for the accumulated gradients, reducing the
dependence on earlier gradients and prioritising the more recent history of
the descent through the parameter space (e.g., when RMSProp finds a highly
convex part of the parameter space, it will very rapidly converge on the
minimum, whereas AdaGrad may never make it into such a region of the parameter
space). Let $\rho_{1}\in[0,1)$ be a decay parameter. The RMSProp algorithm is

$\displaystyle r_{jk}^{\ell}$
$\displaystyle=\rho_{1}r_{jk}^{\ell}+(1-\rho_{1})\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right]^{*}\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ $\displaystyle\Delta w_{jk}^{\ell}$
$\displaystyle=-\eta\frac{1}{\delta+\sqrt{r_{jk}^{\ell}}}\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right].$ (70)

###  D.4 Adaptive moment estimation (Adam)

Adaptive moment estimation (Adam) is an extension of RMSProp and is considered
to be the state of the art optimisation algorithm due to its fast and reliable
convergence. In addition to the accumulation of the square of the gradients
(or square of the modulus of the gradients), Adam also computes an
accumulation of the gradients, decayed by some parameter $\rho_{2}$

$\displaystyle r_{jk}^{\ell}$
$\displaystyle=\rho_{1}r_{jk}^{\ell}+(1-\rho_{1})\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right]^{*}\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ $\displaystyle s_{jk}^{\ell}$
$\displaystyle=\rho_{2}s_{jk}^{\ell}+(1-\rho_{2})\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ $\displaystyle\Delta w_{jk}^{\ell}$
$\displaystyle=-\eta\frac{s_{jk}^{\ell}}{\delta+\sqrt{r_{jk}^{\ell}}},$ (71)

where $r_{jk}^{\ell}$ and $s_{jk}^{\ell}$ are both initialised to be zero at
epoch 0. Since the accumulated gradient and square of the gradients are
initially zero, there is a tendency for training to be biased towards zero
values, which is undesirable. Adam introduces correction terms to
$r_{jk}^{\ell}$ and $s_{jk}^{\ell}$ to avoid such bias,

$\displaystyle r_{jk}^{\ell}$
$\displaystyle=\rho_{1}r_{jk}^{\ell}+(1-\rho_{1})\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right]^{*}\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ $\displaystyle s_{jk}^{\ell}$
$\displaystyle=\rho_{2}s_{jk}^{\ell}+(1-\rho_{2})\left[\frac{\partial\mathcal{C}(\boldsymbol{\theta})}{\partial
w_{jk}^{\ell}}\right],$ $\displaystyle\hat{r}_{jk}^{\ell}$
$\displaystyle=\frac{r_{jk}^{\ell}}{1-\rho_{1}^{\mathcal{E}}}$
$\displaystyle\hat{s}_{jk}^{\ell}$
$\displaystyle=\frac{s_{jk}^{\ell}}{1-\rho_{2}^{\mathcal{E}}}$
$\displaystyle\Delta w_{jk}^{\ell}$
$\displaystyle=-\eta\frac{\hat{s}_{jk}^{\ell}}{\delta+\sqrt{\hat{r}_{jk}^{\ell}}},$
(72)

where lines 3 and 4 are the corrected accumulated gradients (and square of
gradients) and the superscript denotes the power of the current epoch number
$\mathcal{E}$. As with gradient descent, the biases are updated in a similar
way.

###  D.5 Backwards propagation of errors with scalar-valued neurons

In all learning algorithms — from gradient descent to Adam — we need knowledge
of the quantities

$\frac{\partial\mathcal{C}}{\partial
w_{jk}^{\ell}}\quad\text{and}\quad\frac{\partial\mathcal{C}}{\partial
b_{j}^{\ell}},$ (73)

which we refer to as the _sensitivity_ of the cost function with respect to a
given weight or bias in the network. When we have completed a forward pass
through the network, we have a number of outputs which need evaluating against
our already defined cost function. The output layer neuron(s) are dependent
upon all neurons which come before it. Since a neural network is a composition
of many functions, we note that the quantities in Eq. (73) can be determined
from the chain rule and simplified,

$\frac{\partial\mathcal{C}}{\partial
w_{jk}^{\ell}}=\frac{\partial\mathcal{C}}{\partial z_{j}^{\ell}}\frac{\partial
z_{j}^{\ell}}{\partial w_{jk}^{\ell}}=\frac{\partial\mathcal{C}}{\partial
z_{j}^{\ell}}y_{k}^{\ell-1},$ (74)

and

$\frac{\partial\mathcal{C}}{\partial
b_{j}^{\ell}}=\frac{\partial\mathcal{C}}{\partial z_{j}^{\ell}}\frac{\partial
z_{j}^{\ell}}{\partial b_{j}^{\ell}}=\frac{\partial\mathcal{C}}{\partial
z_{j}^{\ell}}.$ (75)

We now seek an equation for ${\partial\mathcal{C}}/{\partial z_{j}^{\ell}}$.
This quantity compares a neuron’s pre-activation against the cost function and
is used to make adjustments to the weights and biases in the network. We will
continue with the convention that neurons with index $k$ feed into neurons
with index $j$.

We can express the sensitivity in the output layer, $L$, by looking at how
small or large the rate of change of the cost function is with respect to the
weighted input $z_{j}^{L}$ (we could, in principle, look at the rate of change
with respect to the full output of the neuron, but this involves an extra
computational step for the activation function, and the weighted input
contains all the information we need). Thus,

$\frac{\partial\mathcal{C}}{\partial
z_{j}^{L}}=\frac{\partial\mathcal{C}}{\partial y_{j}^{L}}\frac{\partial
y_{j}^{L}}{\partial z_{j}^{L}}=\frac{\partial\mathcal{C}}{\partial
y_{j}^{L}}\mathcal{A}^{\prime}(z_{j}^{L}).$ (76)

If Eq. (76) is small, then the neuron is nearly optimal.

We then work top to bottom through the neural network (hence, backpropagation)
in order to determine the optimality of each neuron. For the next layer,
$L-1$, the sensitivity of the cost function with respect to the neuron pre-
activations is

$\frac{\partial\mathcal{C}}{\partial
z_{k}^{L-1}}=\sum_{j=1}^{n_{L}}\frac{\partial\mathcal{C}}{\partial
z_{j}^{L}}\frac{\partial z_{j}^{L}}{\partial z_{k}^{L-1}},$ (77)

where

$z_{j}^{L}=\sum_{k=1}^{n_{L-1}}w_{jk}^{L}\mathcal{A}(z_{k}^{L-1})+b_{j}^{L}.$
(78)

Using Eq. (78) in Eq. (77) gives the sensitivity of the cost function with
respect to a neuron’s pre-activation in the penultimate layer

$\frac{\partial\mathcal{C}}{\partial
z_{k}^{L-1}}=\sum_{j=1}^{n_{L}}\frac{\partial\mathcal{C}}{\partial
z_{j}^{L}}w_{jk}^{L}\mathcal{A}^{\prime}(z_{k}^{L-1}).$ (79)

Then, in any arbitrary layer $\ell-1$, $\ell\in\\{1,\cdots,L-2\\}$, the
sensitivity of the cost function with respect to any arbitrary neuron in that
layer is

$\frac{\partial\mathcal{C}}{\partial
z_{k}^{\ell-1}}=\sum_{j=1}^{n_{\ell}}\frac{\partial\mathcal{C}}{\partial
z_{j}^{\ell}}w_{jk}^{\ell}\mathcal{A}^{\prime}(z_{k}^{\ell-1}).$ (80)

The parameter updates in Eqs. (75) and (74) are obtained by a possibly very
lengthy chain rule. In machine learning, we use _automatic differentiation_ to
compute these quantities 444Automatic differentiation is explored in
literature such as Baydin, _et al._ [47] and Griewank and Walther [25].

Figure 17: A selection of probabilistic activation functions are shown: the
ReLU (Eq. (51)), GELU (Eq. (52b)), SechLU (Eq. (53b)), CauchyLU (Eq. (54b)),
LaplaceLU (Eq. (55b)) and Mish (Eq. (56b)). All figures use either the
variance or scaling parameter $\rho=1$ (even if this is not the most optimal
value). The first column is the probability distribution $p(x)$. The second
column is the cumulative distribution
$\phi(x)=\int_{-\infty}^{x}\mathrm{d}tp(t)$. The third column is the
probabilistic activation function $\mathcal{A}(x)=x\phi(x)$. The fourth column
is the derivative of the activation function $f^{\prime}(x)$. The functions,
with the exception of ReLU, are plotted on common axes to see clearly the
relative magnitude and curvature of the functions.

## References

  * Lagaris  _et al._ [1998] I. Lagaris, A. Likas, and D. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE Trans. Neural Netw. Learn. Syst. 9, 987 (1998).
  * Karniadakis  _et al._ [2021] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nat. Rev. Phys. 3, 422 (2021).
  * Han  _et al._ [2018] J. Han, A. Jentzen, and W. E, Solving high-dimensional partial differential equations using deep learning, Proc. Natl. Acad. Sci. U.S.A. 115, 8505 (2018).
  * Sirignano and Spiliopoulos [2018] J. Sirignano and K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, J. Comput. Phys. 375, 1339 (2018).
  * Lu  _et al._ [2021] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, DeepXDE: A deep learning library for solving differential equations, SIAM Rev. 63, 208 (2021).
  * Xu  _et al._ [2020] Y. Xu, H. Zhang, Y. Li, K. Zhou, Q. Liu, and J. Kurths, Solving Fokker–Planck equation using deep learning, Chaos 30, 013133 (2020).
  * Zhang and Yuen [2022] Y. Zhang and K.-V. Yuen, Physically guided deep learning solver for time-dependent Fokker–Planck equation, Int. J. Non-Linear Mech. 147, 104202 (2022).
  * Bongard and Lipson [2007] J. Bongard and H. Lipson, Automated reverse engineering of nonlinear dynamical systems, Proc. Natl. Acad. Sci. U.S.A. 104, 9943 (2007).
  * Schmidt and Lipson [2009] M. Schmidt and H. Lipson, Distilling free-form natural laws from experimental data, Science 324, 81 (2009).
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proc. Natl. Acad. Sci. U.S.A. 113, 3932 (2016).
  * Raissi  _et al._ [2020] M. Raissi, A. Yazdani, and G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367, 1026 (2020).
  * Meng and Karniadakis [2020] X. Meng and G. E. Karniadakis, A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems, J. Comput. Phys. 401, 109020 (2020).
  * Chen  _et al._ [2019] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, Neural ordinary differential equations (2019).
  * Greydanus  _et al._ [2019] S. Greydanus, M. Dzamba, and J. Yosinski, Hamiltonian neural networks (2019).
  * Cranmer  _et al._ [2020] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho, Lagrangian neural networks (2020).
  * Lecun  _et al._ [1995] Y. Lecun, L. Jackel, L. Bottou, C. Cortes, J. Denker, H. Drucker, I. Guyon, U. Muller, E. Sackinger, P. Simard, and V. Vapnik, Learning algorithms for classification: A comparison on handwritten digit recognition (World Scientific, 1995) pp. 261–276.
  * Ling  _et al._ [2016] J. Ling, A. Kurzawski, and J. Templeton, Reynolds averaged turbulence modelling using deep neural networks with embedded invariance, J. Fluid Mech. 807, 155 (2016).
  * Wang and Yu [2023] R. Wang and R. Yu, Physics-guided deep learning for dynamical systems: A survey (2023).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378, 686 (2019).
  * Vapnik [1995] V. N. Vapnik,  _The Nature of Statistical Learning Theory_ (Springer New York, 1995).
  * Vapnik [1982] V. N. Vapnik,  _Estimation of Dependences Based on Empirical Data_ (Springer-Verlag, Berlin, 1982).
  * Hénon and Heiles [1964] M. Hénon and C. Heiles, The applicability of the third integral of motion: Some numerical experiments, Astron. J. 69, 73 (1964).
  * Bastos de Figueiredo  _et al._ [1998] J. Bastos de Figueiredo, C. Grotta Ragazzo, and C. Malta, Two important numbers in the hénon-heiles dynamics, Phys. Lett. A 241, 35 (1998).
  * He  _et al._ [2015] K. He, X. Zhang, S. Ren, and J. Sun, Delving deep into rectifiers: Surpassing human-level performance on imagenet classification (2015).
  * Griewank and Walther [2008] A. Griewank and A. Walther,  _Evaluating Derivatives_ , 2nd ed. (Society for Industrial and Applied Mathematics, 2008).
  * Griffiths [2024] J. Griffiths, Neural initial value problems codebase (2024).
  * Note [1] The initial velocity of the first coordinate is constrained by $\dot{x}(0)=\sqrt{2E-v_{y,0}^{2}-s_{y,0}^{2}+2s_{y,0}^{3}/3}$, where we have set either $E=1/6$ or $E=1/12$. Given these constraints, the initial position and initial velocity for the second coordinate can be arbitrarily chosen, providing that $\dot{x}(0)$ remains real-valued. To ensure an appropriate choice of the initial conditions for the second coordinate, we find $y(0)$ by using root finding algorithms on the constraint $V<E$, and setting the result as 10% of this value. We find $\dot{y}(0)$ by evaluating $T<E-V$, and setting the result as 10% of this value.
  * Hornik  _et al._ [1989] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators, Neural Networks 2, 359 (1989).
  * Pinkus [1999] A. Pinkus, Approximation theory of the mlp model in neural networks, Acta Numerica 8, 143 (1999).
  * Leshno  _et al._ [1993] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, Multilayer feedforward networks with a nonpolynomial activation function can approximate any function, Neural Networks 6, 861 (1993).
  * Note [2] The universal approximation theorems do not tell us _how_ to find the optimal weights and biases, only that they exist.
  * Wolpert and Macready [1997] D. Wolpert and W. Macready, No free lunch theorems for optimization, IEEE Trans. Evol. Comput. 1, 67 (1997).
  * Breiman [1996] L. Breiman, Bagging predictors, Mach. Learn. 24, 123 (1996).
  * Dietterich [2000] T. G. Dietterich, Ensemble methods in machine learning, in  _Multiple Classifier Systems_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2000) pp. 1–15.
  * Hansen and Salamon [1990] L. Hansen and P. Salamon, Neural network ensembles, IEEE Trans. Pattern Anal. Mach. 12, 993 (1990).
  * Krogh and Vedelsby [1994] A. Krogh and J. Vedelsby, Neural network ensembles, cross validation, and active learning, in  _Adv. Neural Inf. Process. Syst. 32_ , Vol. 7, edited by G. Tesauro, D. Touretzky, and T. Leen (MIT Press, 1994).
  * Opitz and Maclin [1999] D. Opitz and R. Maclin, Popular ensemble methods: An empirical study, J. Artif. Int. Res. 11, 169 (1999).
  * Srivastava  _et al._ [2014] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, Dropout: A simple way to prevent neural networks from overfitting, J. Mach. Learn. Res. 15, 1929 (2014).
  * Hendrycks and Gimpel [2023] D. Hendrycks and K. Gimpel, Gaussian error linear units (GELUs) (2023).
  * Note [3] Probability distributions with support only on the domain $x\in(0,\infty)$ are not conducive to learning, possibly due to the lack of symmetry. It is our empirical observation that some distributions which are not perfectly symmetric about the origin — such as that used in Mish (see Eq. (56b)) — are still conducive to learning.
  * Misra [2020] D. Misra, Mish: A self regularized non-monotonic activation function (2020).
  * Glorot and Bengio [2010] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks (Proceedings of the 13th International Conference on Artificial Intelligence and Statistics (AISTATS), Chia Laguna Resort, Sardinia, Italy, 2010).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014).
  * Duchi  _et al._ [2011] J. Duchi, E. Hazan, and Y. Singer, Adaptive subgradient methods for online learning and stochastic optimization, J. Mach. Learn. Res. 12, 2121 (2011).
  * Hinton [2018] G. Hinton, Coursera: Neural networks for machine learning (lecture 6) (2018).
  * Note [4] Automatic differentiation is explored in literature such as Baydin, _et al._ [47] and Griewank and Walther [25].
  * Baydin  _et al._ [2017] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: A survey, J. Mach. Learn. Res. 18, 5595 (2017).