# Deep Learning: Hydrodynamics, and Lie–Poisson Hamilton–Jacobi Theory

Nader T. Ganaba E-mail address: mailto:nganaba<EMAIL_ADDRESS>
###### Abstract

The interpretation of deep learning as a dynamical system has gained a
considerable attention in recent years as it provides a promising framework.
It allows for the use of existing ideas from established fields of mathematics
for studying deep neural networks. In this article we present deep learning as
an equivalent hydrodynamics formulation which in fact possess a Lie-Poisson
structure and this further allows for using Poisson geometry to study deep
learning. This is possible by considering the training of a deep neural
network as a stochastic optimal control problem, which is solved using
mean–field type control. The optimality conditions for the stochastic optimal
control problem yield a system of partial differential equations, which we
reduce to a system of equations of quantum hydrodynamics. We further take the
hydrodynamics equivalence to show that, with conditions imposed on the
network, that it is equivalent to the nonlinear Schrödinger equation.

As the such systems have a particular geometric structure, we also present a
structure–preserving numerical scheme based on the Hamilton–Jacobi theory and
its Lie–Poisson version known as Lie–Poisson Hamilton–Jacobi theory. The
choice of this scheme was decided by the fact it can be applied to both
deterministic and mean–field type control Hamiltonians. Not only this method
eliminates the possibility of the emergence of fictitious solutions, it also
eliminates the need for additional constraints when constructing high order
methods.

##  1 Introduction

The kinship between deep learning and dynamical systems, discussed thoroughly
and lucidly in [102], has gained a considerable attention in recent years as
it was studied in a number of papers such as [67, 68, 10]. This comes as a
consequence of the increased complexity of deep learning and their success in
various real world applications, such as performing outstandingly in
recognition tasks for images [30, 99, 62], and speech [42, 93], automated
language translation [98], drug discovery [69, 34, 22] and predicting splicing
patterns of RNA and their role in genetic diseases [8, 66, 104]. Yet, the
inner workings are not well understood and the lack of mathematical framework
in which one can analyse the systems. Further, basic tools of linear algebra
are not sufficient to explain and analyse the behaviour of complex systems.

This article is part of a series that view deep learning as a hydrodynamics
system within the geometric mechanics framework to better understand the
behaviour of deep networks. And here, we introduce a Hamiltonian hydrodynamic
system on Poisson manifold which is equivalent to training of deep neural
networks. The result is a system of partial differential equations similar to
quantum Euler equation, or Madelung equation and this equation can be further
transformed into the nonlinear Schrödinger equation. As these equations have
special properties, we will also present geometric numerical schemes to be
used as training algorithms.

The passage from deep learning to hydrodynamics is not direct and at its core
is optimal control interpretation of deep learning. The use of optimal control
method for deep learning is not new, in fact, it is used to train residual
neural networks [67, 103, 68, 10]. In addition to that, it is a candidate for
a framework for studying deep neural networks [40, 10]. Besides training this
point of view is important as there is an association between optimal control
and geometry [97, 13, 35] and this offers the opportunity to use methods and
theories beyond linear and nonlinear algebra to study deep learning. This
association is due in part to Pontryagin’s Maximum Principle [89] for solving
the optimal control problem, where it gives a Hamiltonian function whose
associated vector field generates a flow that solves the training problem. In
Hamiltonian mechanics, the solutions can be considered as coordinate
transformations between two symplectic manifolds known as _symplectic
transformations_ , or _symplectomorphisms_. The underlying geometry allows one
to make the transition into continuum mechanics [14, 43] and such transition
is reminiscent of [9] use of fluid dynamics for solving optimal transport
problems. Here, we consider a general case when the deep neural network is
stochastic as a way to quantify uncertainties. That said, the training problem
becomes a stochastic optimal control problem and there are different ways to
solve this problem. The approach chosen here is mean–field type control [11],
which casts the problem as a deterministic problem. What is advantageous is
that the conditions for optimality can be reduced to an equivalent
hydrodynamics system, which is indeed a _Lie–Poisson system_. The resultant
system of partial differential equations has a structure as the one for ideal
compressible fluid.

Such formulation for deep learning is novel, and it is presented in this
article. Not only it helps in applying Hamilton–Jacobi for Poisson
integrators, but also provides a framework for studying and analysing deep
learning problems. The said Lie–Poisson system for deep learning presented
here has the same structure as the ideal compressible fluids; it allows us to
use theories of hydrodynamics to study neural networks. The question of
stability can help choosing hyperparameters heuristically and methods such as
energy–momentum [95, 94], and energy–Casimir both deterministic [45] and
stochastic [4] can be used to study the stability of Lie–Poisson systems. That
said, the study of stability is reserved for a future article. Here the focus
is on employing Hamilton–Jacobi theory and its variants for systems with
symmetry for deriving structure–preserving integrators for Hamiltonian
systems.

Since these systems are Hamiltonian, one needs to use numerical schemes that
are symplectic or Poisson. There are a number of ways one can derive such
integrators, such as the partitioned Runge–Kutta method [1], or variational
integrators [78] and here we focus on a type that uses generating functions.
It relies on the _Hamilton–Jacobi theory_ , which is an alternative
formulation of Hamiltonian mechanics in terms of wavefronts and instead of
solving for a single path, it generates a family of solutions. At the heart of
the Hamilton–Jacobi theory, there is the _Hamilton–Jacobi equation_ , as it
appears in geometrical optics [19, 16], and the equation for modelling the
wave front propagation and this is considered a passage from quantum mechanics
to classical mechanics [53]. The solution of the Hamilton–Jacobi equation,
known as _generating function_ , which from the point of view of symplectic
geometry, allows to heuristically construct coordinate transformations that
preserve a geometric structure known as _symplectic structure_. This forms a
basis for the numerical scheme, as obtaining an approximation of the
generating function we can then use it to generate symplectic transformations.
In other words, the numerical scheme is a symplectic transformation.

Here, we limit the application of the Hamilton–Jacobi theory to the
construction of a structure–preserving numerical scheme, based on the
generating function, by projecting a continuous neural network model onto a
finite dimensional space while conserves a number of the network’s properties.
The second application is more elaborate as it takes into account
uncertainties in the network and uses mean–field games interpretation of deep
learning and it uses a variant of Hamilton–Jacobi theory known as _Lie–Poisson
Hamilton–Jacobi theory_ [107]. As mentioned, the resultant system of partial
differential equations has the same Poisson structure as the ideal
compressible fluid [76], and because it is defined on semidirect product
spaces, it necessitates the derivation of Lie–Poisson Hamilton–Jacobi theory
of semidirect products spaces.

Most numerical methods for the class of Lie–Poisson equations for infinite
dimensional systems use a discretisation of the vector field on the group of
diffeomorphisms and not the group’s elements [79, 29] based on the works of
[106]. However, it is difficult to discretise the group of diffeomorphisms
directly, as it is an infinite dimensional group, and apply Lie–Poisson
Hamilton–Jacobi theory. One option is to consider the discretisation of [87]
where the group of diffeomorphisms is approximated in finite dimensional space
as the Lie group $GL(N)$. Such discretisation comes with its own challenges,
such as the introduction of flat operators. For this reason we introduce a
discretisation that uses the group of pseudodifferential symbols instead of
the group of diffeomorphisms, which is a submanifold of the group of
pseudodifferential symbols, and it makes systems on the infinite dimensional
space compatible with the Lie–Poisson Hamilton–Jacobi theory. The reward of
going to these lengths for numerical integrators is that the scheme preserves
important properties such as the coadjoint orbit, meaning that the numerical
solution would always belong to the solution set. Figure 1 demonstrates
visually the solution of the angular momentum of Hamiltonian system of the
$\mathrm{SO}(3)$ rigid body obtained by three different numerical schemes. For
this rigid body, the coadjoint orbit is a sphere. It is evident the explicit
Euler and $4^{th}$ order Runge–Kutta method yield solutions that do not belong
to the set of solutions, whereas the Lie–Poisson Hamilton–Jacobi scheme yields
a solution that remains on the surface of the sphere.

The three main reasons for using generating functions to derive numerical
schemes are high order integrators, possible lack of hyperregularity and its
extension to the stochastic case. The first reason concerns the construction
of high-order schemes and using generating functions one can construct high-
order methods without the need to resort to carefully select coefficients,
such as symplectic Runge–Kutta method [1], or compose a number low order
integrators to form a high order one [105]. The second reason is due to the
nature of the optimal control Hamiltonian. In Hamiltonian mechanics, the move
into the Lagrangian formulation is done by using Legendre transform, however,
this transform is not bijective when the Hamiltonian is not hyperregular. One
approach to solving Hamiltonian systems numerically is to move into the
Lagrangian side and apply variational integrators [78]. Moreover, the
application of structure–preserving integrators for optimal control was first
done in the Lagrangian picture [48] and further studied in [63, 64, 15, 58].
The lack of hyperregularity impedes the move to the Lagrangian side and for
this reason a Hamiltonian structure–preserving integrator is needed [65]. The
third reason is of practicality as our approach can be applied to both
deterministic and stochastic optimal control problems, albeit with additional
steps.

(a)

(b)

(c)

Figure 1: The solution of $SO(3)$ Rigid body, which is a Lie–Poisson system,
computed using three different numerical methods to highlight the importance
of structure–preserving integrators. The equation of motion is
$\dot{\Pi}=\Pi\times\Omega$, with $\Pi=\mathbb{I}\Omega$ and
$\mathbb{I}=\mathrm{diag}(1,2,3)$ is the inertia tensor. In figure 1(a), the
Euler scheme is used to approximate the solution and it is evident that it
does not stay on the sphere. The second scheme used is the $4^{th}$ order
Runge–Kutta method, shown in figure 1(b), although it fares better than Euler
method, it too does not preserve the structure. The last scheme is
structure–preserving scheme that uses Lie–Poisson Hamilton–Jacobi theory of
[107], the method is discussed in section 3.1, and the results that it yields,
shown in figure 1(c), are satisfactory.

##  2 Preliminaries

###  2.1 Hamilton–Jacobi Theory

As it plays central role for the derivation of structure–preserving schemes,
here, we give a brief summary of the Hamilton–Jacobi theory and the important
concepts used throughout this article. One source that we highly recommend is
[77]. We consider the the space of all possible states of the dynamical system
and it is called _configuration manifold_ $Q$ with coordinates
$q=(q^{1},\dots,q^{n})$, where $n$ is the dimension of $Q$. The tangent space
of the configuration bundle is $TQ$ and its dual is the cotangent bundle
$T^{\ast}Q$, whose generalised coordinates are
$(q,p)=(q^{1},\dots,q^{n},p^{1},\dots,p^{n})$. In mechanics, there are a
number of ways to describe the motion of the system and one of these is
Newtonian mechanics. Despite the fact that the intuition behind Newtonian
mechanics is easier to grasp, when there are physical and geometric
constraints, it becomes difficult to model. Lagrangian mechanics was
introduced to tackle this and it is a formulation that is widely used for
constraint particle motion. Such formulation of mechanics requires
$(q,\dot{q})\in TQ$ to describe the motion of the system. An alternative
formulation, and the one used here, is the Hamiltonian formulation which uses
position and conjugate momenta $(q,p)\in T^{\ast}Q$ to describe evolution of
the system. Another important object in Hamiltonian mechanics is the
Hamiltonian function $H:T^{\ast}Q\to\mathbb{R}$, which carries all the
information of the dynamics of the system. The Hamiltonian, $H$, is related to
the energy of the system and, in some cases, it coincides with the total
energy. Once we have the Hamiltonian, $H$, and with fixed boundary terms
$q(0)$ and $q(T)$, the variational principle states that

$\displaystyle\delta\int_{0}^{T}\left\langle
p,\dot{q}\right\rangle-H(q,p)\,dt=0,$ (1)

and it yields Hamilton’s equation

$\displaystyle\frac{d}{dt}q=\frac{\partial H}{\partial
p}(q,p),\quad\text{and}\quad\frac{d}{dt}p=-\frac{\partial H}{\partial
q}(q,p).$ (2)

Often case, Hamilton’s equation are described as what is known as _a
Hamiltonian vector field_ $X_{H}=(\frac{\partial H}{\partial
p},-\frac{\partial H}{\partial q})$ whose flow is a symplectomorhism,
$\Phi_{X_{H}}:(q_{0},p_{0})\to(q,p)$ , and it maps the initial points
$(q_{0},p_{0})\in T^{\ast}Q$ to the final points $(q(t),p(t))\in T^{\ast}Q$ at
time $t$. This mapping leaves the symplectic structure, $\omega=dq\wedge dp$
invariant, i.e.

$\displaystyle\left(\Phi_{X_{H}}\right)^{\ast}\omega-\omega=0,$ (3)

and this allows for constructing coordinate transformation that considerably
simplifies the equation of motion. The mapping
$\left(\Phi_{X_{H}}\right)^{\ast}$ is the pull-back of of $\omega$ by flow of
$X_{H}$ and (3) can be easily proven by taking the time derivative and using
Cartan’s formula for Lie derivatives [77].

In fact, one choice of coordinates $(w,J)$, known as action-angle coordinates,
chooses the momentum $J$ to be a constant of motion and this coordinate system
is central to the celebrated _Kolmogorov-Arnold-Moser theorem_. Meanwhile, the
symplectic form can be written as an exact form

$\displaystyle\omega=\mathbf{d}\theta=dq\wedge
dp,\quad\text{where}\quad\theta=p\mathbf{d}q$ (4)

where $\mathbf{d}$ is a map from $k$-forms on a manifold to $(k+1)$-forms on
the same manifold known as _exterior derivative_ and it is linear and
satisfies the product rule [2]. In coordinates, if a $k$-form is given by

$\displaystyle\alpha=\alpha_{i_{1},\dots,i_{k}}dx^{i_{1}}\wedge\dots\wedge
dx^{i_{k}},$

the exterior derivative is given by

$\displaystyle\mathbf{d}\alpha=\sum_{j=1}^{k}\frac{\partial\alpha_{i_{1},\dots,i_{k}}}{\partial
x^{j}}dx^{j}\wedge dx^{i_{1}}\wedge\dots\wedge dx^{i_{k}}$

It has the properties that $\mathbf{d}^{2}=0$ and when $f$ is a scalar
function, $\mathbf{d}f$ is equal to the differential of $f$. The relation (4)
holds is due to the fact that $\omega$ is a closed form, i.e.
$\mathbf{d}\omega=0$. This is important, as it allows us to write

$\displaystyle\left(\Phi_{X_{H}}\right)^{\ast}\mathbf{d}\theta-\mathbf{d}\theta=\mathbf{d}\left(\left(\Phi_{X_{H}}\right)^{\ast}\theta-\theta\right)=0=\mathbf{d}^{2}S$
(5)

and it implies there the difference between the one-form and its pullback by
the flow of the Hamiltonian system is a closed one-form. This closed one-form
is defined as an exterior derivative of the function, $S:Q\times
Q\times\mathbb{R}\to\mathbb{R}$ , and it is known as _generating function_.
The generation of symplectic transformation from $S$ comes from the relation

$\displaystyle\mathbf{d}S=p\mathbf{d}q-p_{0}\mathbf{d}q_{0},$ (6)

which yields

$\displaystyle p(t)=\frac{\partial S}{\partial
q}(q_{0},q(t),t),\quad\text{and}\quad p_{0}=\frac{\partial S}{\partial
q_{0}}(q_{0},q(t),t).$ (7)

Luckily, the generating function satisfies a partial differential equation

$\displaystyle\frac{\partial S}{\partial t}+H\left(q,\frac{\partial
S}{\partial q}\right)=0.$ (8)

known as _Hamilton–Jacobi equation_. Solving Hamilton’s equations and solving
Hamilton–Jacobi equation is one and the same thing. The main difference
between the two approaches is that the latter deals with functions defined on
the configuration manifold, $Q$, while the former the dynamics belong to the
cotangent bundle $T^{\ast}Q$.

For certain systems, it is taxing to work with generating function $S$ that
generates a canonical transformation from the old to the new coordinates,
$(q_{0},p_{0})$ and $(q(t),p(t))$ respectively, and for this it is not unwise
to consider a combination of old and new coordinates. The said combinations
can be constructed using Legendre transformations and one transformation that
is widely used is the type-II generating function

$\displaystyle S_{2}(q_{0},p(t),t)$
$\displaystyle=S(q_{0},q(t),t)+\left\langle q(t),p(t)\right\rangle.$ (9)

Substituting $S$ in terms of $S_{2}$ we obtain a Hamilton–Jacobi equation (8)
for each type of generating function.

######  Remark 2.1.

For numerical schemes based on generating functions, it is practical to use
type-II generating function as it allows for initial condition to be
$S_{2}(p_{0},p(0),0)=0$ and, in the same time, satisfy the identity
transformation required by the Hamilton–Jacobi theory.

###  2.2 Optimal Control Theory

In this section, we give a brief overview of optimal control theory and
formulating deep learning as an optimal control problem. The reader is invited
to read of the following references [3, 49, 17] for more detailed exposition
on the subject. Let $Q$ be a configuration manifold and a system’s state, on
the time interval $[0,T]$, is described by $q(t)\in Q$. The goal of optimal
control is to move the system from the state $q(0)=q_{0}$ to the state
$q(T)=q_{T}$, where $q_{0},q_{T}\in Q$ are predetermined states. This is
actualised by choosing $\theta=(u,b)\in\mathcal{U}$, where $\mathcal{U}$ is
called _the space of admissible controls_ , such that the _cost functional_

$\displaystyle S(q,\theta)=\int_{0}^{T}\ell(q(t),\dot{q(t)},\theta(t))\,dt$
(10)

is maximised/minimised and at the same time $q(t)\in Q$ satisfies

$\displaystyle\frac{d}{dt}q=\sigma\left(u\cdot q+b\right),$ (11)

where $q\in Q$, $\frac{d}{dt}q\in TQ$ and $\theta=(u,b)\in\mathcal{U}$ and
$\sigma:Q\times\mathcal{U}\to\mathbb{R}$. The term
$\ell:TQ\times\mathcal{U}\to\mathbb{R}$, in (10), is the _loss function_ or
_cost Lagrangian_. The choice of loss function depends on the application, and
it is usually chosen as a measure of how far is the state of the system,
$q(t)$, from a predetermined target. There are two methods that are widely
used to compute $\theta$: dynamic programming and Pontryagin’s Maximum
Principle. The former breaks the problem into a smaller problem on the
interval $[t,t+\Delta t]$ and uses the solution of the smaller problem to
construct a solution on $[0,T]$ and it yields a solution which satisfies the
necessary and sufficient condition for optimality. However, it is not
considered in this article and instead we use Pontryagin’s Maximum Principle.

The method in essence is finding a curve on the Pontryagin’s bundle
$T^{\ast}Q\times\mathcal{U}$, denoted by $(q(t),p(t),\theta(t))$, that
minimises (10) and it also satisfies the following Hamilton’s equation:

$\displaystyle\frac{d}{dt}q(t)=\frac{\partial H}{\partial
p},\quad\frac{d}{dt}p(t)=-\frac{\partial H}{\partial
q},\quad\mbox{and}\quad\frac{\partial H}{\partial\theta}=0,$ (12)

with the Hamiltonian function $H:T^{\ast}Q\times\mathcal{U}\to\mathbb{R}$,
defined by $H(q,p,\theta)=\left\langle
p,\sigma\left(q(t),\theta(t)\right)\right\rangle-\ell\left(q(t),\theta(t)\right)$.
These equations, (12), are the stationary variations of Hamilton’s principle

$\displaystyle 0=\delta\mathbb{S}=\delta\int_{0}^{T}\left\langle
p,\dot{q}\right\rangle-H(q,p,\theta)\,dt.$ (13)

The equation $\dot{p}=-\frac{\partial H}{\partial q}$ is the costate equation
and also known as the adjoint equation and it has gained attention in deep
learning research as it is the core of methods such as neural ordinary
differential equation [23, 96, 25]. Solving such Hamiltonian system is an
arduous task, thus numerical methods are required to approximate the solution.
The application of numerical integrators for optimal control problems has been
explored in [48, 63, 64, 15, 58]and the key idea is to use a discrete version
of the variational principle. On the contrary, in this article, the idea is,
instead of applying discrete variational principle, we use the
structure–preserving numerical scheme derived using Hamilton–Jacobi theory,
given by (82), to approximate the solution of Pontryagin’s Maximum Principle
(12). Concerning the initial conditions for $(q(t),p(t))$, we are only given
$q_{0}$ and $q_{T}$ and for this reason an additional step is needed to find
the correct $p^{0}$ that would make $(q^{k},p^{k},\theta^{k})$, where
$k=0,\dots N_{t}$ satisfy $q^{0}=q_{0}$ and $q^{N_{t}}=q_{T}$. There are a
number of choices for the extra step, such as the shooting method where
$p^{0}$ is chosen at random and then adjusted. The other approach is to use
nonlinear programming, and the idea is that it makes the solution of the
algebraic expressions (82) and the boundary conditions as the solution of a
nonlinear system of equations.

###  2.3 Deep Neural Network Training As Optimal Control Problem

Now considering deep learning, we first need to a representation of the neural
networks as a dynamical system. Artificial neural networks can be viewed as a
group of simple sub-systems referred to as synapses connected in a specific
way. Mathematically, a synapsis is modelled as the following transformation
$Q_{k}\to Q_{k+1}$, without the loss of generality, in this article we solely
consider the case where this transformation is an endomorphism, i.e. a
transformation $Q\to Q$, given by

$\displaystyle q_{i}^{k+1}$
$\displaystyle=F^{n}\left(q_{i}^{k},\theta^{k}\right),\quad\mbox{for}\quad
k=0,\ldots,N_{t}-1$ (14)

where $q^{k+1}_{i}$ is the $i^{th}$ synapsis of the $k+1$ layer, $q^{k}\in Q$,
$\theta^{k}=(u^{k},b^{k})\in\mathcal{U}$ is the network’s parameters and the
map $F^{k}:Q\times\mathcal{U}\to Q$ is the state transition map. The vector
$q^{k+1}=(q^{k+1}_{1},\ldots,q^{k+1}_{N_{l}})\in Q\subset\mathbb{R}^{N_{l}}$
constitute a single layer in the network. For general cases, the configuration
space $Q$ need not be vector spaces, they can be Hilbert spaces or
differentiable manifold, and in that case (14) can be regarded as a pullback
operation from $Q_{k}$ to $Q_{k+1}$. The choice of the state transition map
determines the type of network and here we consider a transition of the type

$\displaystyle q_{i}^{k+1}$
$\displaystyle=F\left(q^{k}_{i},\theta^{k}\right)=q^{k}_{i}+\alpha\sigma\left(\sum\limits_{j=0}^{N_{l}}u^{k}_{i,j}q^{k}_{j}+b^{k}_{i}\right),\quad\alpha\in\mathbb{R}_{>0},$
(15)

and this renders (14) as residual neural network, ResNet henceforth. The
function $\sigma:Q\times\mathcal{U}\to Q$ is the activation function and it is
applied element-wise and few of the most commonly used ones are: linear
rectifier, denoted by ReLU, hyperbolic tan, sigmoid function, and softmax.
Here, $u^{k}_{i,j}$ is the element of the $i^{th}$ row and $j^{th}$ column of
weight matrix and $b^{k}_{i}$ is the $i^{th}$ element of bias vector for the
$k^{th}$ layer. We do not specify $\theta^{k}=(u^{k},b^{k})$ for
$k=0,\ldots,N_{t}$, before hand, at first, random weights and biases are
sampled form a distribution function. Finding the values for $\theta^{k}$ to
make the neural network model a specific system is the focus of training
methods studied in the field of machine learning.

Equation (15) is viewed as Euler discretisation of a continuous dynamical
system [26, 27, 23, 91]. The move to the continuous domain can be interpreted
as taking the limit when the number of hidden layers approaches infinity,
$N_{l}\to\infty$. The result is a dynamical system defined on the manifold $Q$
of the form

$\displaystyle\frac{d}{dt}q(t)=\sigma(q(t),\theta(t)),\quad
t\in[0,T],\quad\mbox{and}\quad q(0)=q_{0},$ (16)

where $q_{0}\in Q$.

######  Remark 2.2.

In some literature [23, 25], the limit of ResNet, (16), is referred to as
_neural ordinary differential equation_. In this article, we substitute this
naming with continuous ResNet.

Training the neural network can be treated as an optimal control problem with
(16) as a constraint, and the dataset
$\\{(q_{0}^{(i)},c^{(i)})\\}_{i=0,\dots,N}$ specifies the initial conditions
and the final conditions [67]. We define $\mathcal{C}$ as the space of target
data. Thus, we state the training of continuous ResNet as follows

######  Problem 2.3.

Given the training set $\\{(q_{0}^{(i)},c^{(i)})\\}_{i=0,\dots,N}$, compute
$\theta=(u,b)$ such that it minimises

$\displaystyle
S=\int_{0}^{T}\ell_{d}\left(q(t),\theta(t)\right)\,dt+\sum_{i=0}^{N}L(\pi((q(T)^{(i)}),c^{(i)}),$
(17)

such that it satisfies

$\displaystyle\frac{d}{dt}q(t)=\sigma\left(u(q(t))+b(t)\right),$ (18)

and $q(0)=q_{0}^{(i)}\,\text{and}\,\pi(q(T))=c^{(i)}$ for $i=0,\dots,N$.

Here, where $\ell:Q\times\mathcal{U}\to\mathbb{R}$ is the cost Lagrangian and
in machine learning it acts as a regularisation term. The last term of the
functional $S$ is the terminal condition and it is a metric for measuring how
far the output of (57), $q(T)^{(i)}$, whose starting point is
$q(0)=q_{0}^{(i)}$, is from the target data $c^{(i)}$. A suitable candidate is
$\sum\limits_{i=0}^{N}L(\pi(q(T)^{(i)}),c^{(i)})$, where
$L:\mathcal{C}\times\mathcal{C}\to\mathbb{R}$ and the projection
$\pi:Q\to\mathcal{C}$ and in this article it is chosen as
$\pi(x)=\frac{1}{1+e^{-x}}$. As for $L$, in many classification problems, a
common and a suitable choice is the categorical cross entropy loss function.

##  3 Stochastic Optimal Control for Deep Learning

In this section, we consider a continuous residual neural network that takes
into consideration the network’s uncertainties. These uncertainties are
accounted for by modelling them as a stochastic process and for the time-
being, we only consider Wiener processes. That means the dynamical model for
continuous ResNet, given by (16), is now a stochastic differential equation.
The addition of stochastic terms not only render the model robust for
potential variations, but also it is very rich mathematically as shall we
explore in this section.

Similar to the previous sections, we present an algorithm for training deep
learning using a generating function for stochastic continuous ResNet. As the
dynamical system is now stochastic, the application of the Pontryagin’s
Maximum Principle used in this article cannot be done directly. For this we
require stochastic optimal control and especially mean–field type optimal
control [11], and, in this paradigm, the solution to the optimal control
problem satisfies a system of partial differential equations instead of
ordinary differential equations. The idea of using stochastic optimal control
for deep learning was studied in [27, 88], however, the main difference is we
use elimination theorem of [24] to reduce the system of partial differential
equations to hydrodynamic equation with advected quantity.

One major consequence of mean–field optimal control is that we can no longer
utilise the Hamilton–Jacobi theory used in section
LABEL:sec:hamilton_jacobi_intro. However, there are two possible ways to
construct structure–preserving schemes for optimal control: using
multisymplectic integrators and using integrators based on generating
functions on Poisson manifolds. In this article, we only focus on the latter,
as the former is studied in [32]. The key here is to use Poisson map analogous
to the symplectic transformation and we, first, reduce the system of partial
differential equations that we obtain by mean–field type optimal control to a
hydrodynamics equation. Then we restate the said hydrodynamics equation as a
Lie–Poisson equation, which is the same equation but from the Hamiltonian
viewpoint. The Poisson geometry counterpart for the Hamilton–Jacobi theory was
introduced in [107] and it is known as the Lie–Poisson Hamilton–Jacobi theory.

We give a brief overview of the theory in section 3.1 and we introduce a
variant of the theory that concerns semidirect products. In section 3.1.1 we
introduce the algebras of pseudodifferential symbols whose dual space is where
the Lie–Poisson equation evolves. In section 4.1, we use algebras of
pseudodifferential symbols on discrete meshes to derive numerical scheme based
on generating function on Poisson space is presented in proposition 4.11.

###  3.1 Hamilton–Jacobi for Systems with Symmetry

Having demonstrated how Hamilton–Jacobi theory is used as a formulation for
Hamiltonian systems on $Q$, here we consider a special case when the
Hamiltonian is $G$-invariant. What we mean by that, we consider a
configuration manifold $Q$ and when $g\in G$ is a Lie group action acting on
the Hamiltonian $H:T^{\ast}Q\to\mathbb{R}$, the following holds

$\displaystyle H(g\cdot q,(g)^{\ast}p)=H(q,p),\quad(q,p)\in T^{\ast}Q.$

This allows us to perform the Lie–Poisson reduction [77]. We will not cover
reduction theory here, the reader is advised to refer to [77, 76, 72]. In
essence, the invariance allows us to choose action such that $T^{\ast}Q/G$ and
the resultant quotient space is isomorphic to dual of the Lie algebra
$\mathfrak{g}^{\ast}$. For example, when $Q=SO(3)$ and $G=SO(3)$, the quotient
space $T^{\ast}SO(3)/SO(3)$ is isomorphic to $\mathfrak{so}^{\ast}(3)$, the
dual of the Lie-algebra associated to $SO(3)$. In general, the dual of Lie
algebra $\mathfrak{g}^{\ast}$ is an example of what is known as a Poisson
manifold [76]. A manifold $(P,\\{\cdot,\cdot\\})$ is said to be a _Poisson
manifold_ when a manifold $P$ is equipped with a bilinear form
$\\{\cdot,\cdot\\}$ that has the properties that is for any function on
$C^{\infty}(P)$, the form is anti-symmetric, satisfies both the Jacobi
identity and Leibniz’s rule. The bilinear form $\\{\cdot,\cdot\\}$ is known as
_the Poisson structure_. Let $(P_{1},\\{\cdot,\cdot\\}_{1})$ and
$(P_{2},\\{\cdot,\cdot\\}_{2})$ be two Poisson manifolds, then a map
$\phi:P_{1}\to P_{2}$ that satisfies
$\\{f\circ\phi,h\circ\phi\\}_{1}=\\{f,h\\}_{2}\circ\phi$ for any $f,h\in
C^{\infty}(P_{2})$ is called a _Poisson map_. A detailed exposition of
reduction is beyond the scope of this article, however, it is only summarised.
For the dynamics, whose Hamiltonian is $H:T^{\ast}Q\to\mathbb{R}$, its
integral curves when projected to a Poisson manifold $P$ using the map
$\pi_{P}:T^{\ast}Q\to P$, coincide with the integral curves of dynamics, whose
Hamiltonian is $H_{P}:P\to\mathbb{R}$. We can say that the flow $F$ of the
Hamiltonian $H$ and the flow $F_{P}$ of $H_{P}$ are related by $\pi_{P}\circ
F=F_{P}\circ\pi_{P}$.

${T^{\ast}Q}$${T^{\ast}Q}$${P}$${P}$$\scriptstyle{F}$$\scriptstyle{F_{P}}$

The Hamiltonian dynamics on $T^{\ast}Q$ induce Lie–Poisson equation on
$\mathfrak{g}^{\ast}$

$\displaystyle\frac{d}{dt}F=\\{F,H_{P}\\},$ (19)

which is equivalent to

$\displaystyle\frac{d}{dt}\mu=\mathrm{ad}^{\ast}_{\frac{\delta
H_{P}}{\delta\mu}}\mu,$ (20)

where
$\mathrm{ad}^{\ast}:\mathfrak{g}\times\mathfrak{g}^{\ast}\to\mathfrak{g}^{\ast}$
is the coadjoint map.

The Lie–Poisson reduction can also be used for the Hamilton–Jacobi theory
[107]. Starting with the Hamilton–Jacobi theory

$\displaystyle\frac{\partial S}{\partial t}(q,q_{0})+H\left(q,\frac{\partial
S}{\partial q}\right)=0,\quad p_{0}=-\frac{\partial S}{\partial
q_{0}},\,p(t)=\frac{\partial S}{\partial q},$ (21)

where $S:Q\times Q\to\mathbb{R}$ is the generating function, the Hamiltonian
$H:T^{\ast}Q\to\mathbb{R}$ is a $G$-invariant Hamiltonian. The generating
function yields a canonical transformation from $(q_{0},p_{0})\to(q(t),p(t))$.
Let $S_{L}(g)=S(q^{-1}q,q^{-1}q_{0})=S(e,g)$, where $g=q^{-1}q_{0}$, the
Lie–Poisson reduction reduces (21) to

$\displaystyle\frac{\partial S_{L}}{\partial
t}(g)+H\left(-\mathbf{J}_{L}\circ\mathbf{d}S_{L}\right)=0,\quad\mu_{0}=TL_{g}^{\ast}\mathbf{d}_{g}S_{L},\quad\mu=TR_{g}^{\ast}\mathbf{d}_{g}S_{L},$
(22)

with $\mu_{0},\mu\in\mathfrak{g}^{\ast}$. Here we have $S_{L}:Q\to\mathbb{R}$,
$g\in G$ is the Lie group action, $TL_{g}^{\ast}$ and $TR_{g}^{\ast}$ are the
left and right cotangent lift of Lie group action $g$. The map
$\mathbf{J}_{L}:T^{\ast}Q\to\mathfrak{g}^{\ast}$ is the _the momentum map_ for
the cotangent lift and it is equal to

$\displaystyle\mathbf{J}_{L}(\alpha)=T_{e}R_{g}^{\ast}\alpha,$ (23)

while $\mathbf{J}_{R}:T^{\ast}Q\to\mathfrak{g}^{\ast}$ is given by

$\displaystyle\mathbf{J}_{R}(\alpha)=T_{e}L_{g}^{\ast}\alpha.$ (24)

The partial differential equation and two expressions in (22) are known as
_Lie–Poisson Hamilton–Jacobi theory_ introduced in [107]. In fact, the
integrator that was used in figure 1, to showcase one of the advantages of
structure–preserving integrators, uses the discrete version (22) on $SO(3)$.

Now, the generating function $S_{L}$, generates a Poisson map
$\phi:\mu_{0}\to\mu$, which is the analogue of symplectic map on the cotangent
bundle $T^{\ast}Q$. This reduction is not only valid for semisimple Lie groups
such as $SO(3)$, it can be applied to semidirect products. First, let us
recall some basic facts about semidirect products. Let $S=G\oplus V$ be a
semidirect product space with Lie group $G$ and vector space $V$. Let
$\zeta:G\to\mathrm{Aut}(V)$ be the representation of the Lie group $G$ in $V$,
and $\zeta^{\prime}:\mathfrak{g}\to\mathrm{End}(V)$ is the Lie algebra
representation in $V$. The multiplication of $S$ is defined as

$\displaystyle(g_{1},v_{1})(g_{2},v_{2})=\left(g_{1}g_{2},v_{1}+\zeta(g_{1})v_{2}\right),$

where $(g_{1},v_{1}),(g_{2},v_{2})\in S$. The corresponding Lie algebra
$\mathfrak{s}=\mathfrak{g}\times V$ has the Lie bracket

$\displaystyle\left[(\xi_{1},v_{1}),(\xi_{2},v_{2})\right]=\left(\left[\xi_{1},\xi_{2}\right],\zeta^{\prime}(\xi_{1})v_{2}-\zeta^{\prime}(\xi_{2})v_{1}\right),$

where $(\xi_{1},v_{1}),(\xi_{2},v_{2})\in\mathfrak{s}$. The left action of $S$
on to itself is defined as

$\displaystyle L_{(g,v)}(h,u)=\left(gh,v+\zeta(g)u\right),$ (25)

and tangent lift of the left action is

$\displaystyle
T_{(h,u)}L_{(g,v)}(v_{h},u,w)=\left(T_{h}L_{g}(v_{h},v+\zeta(g)u,\zeta(g)w\right),$
(26)

where $(v_{h},u,w)\in T_{(h,u)}\left(G\oplus V\right)$. The tangent lift at
the element $(g,v)(h,u)\in S$ and we are introducing it, so we can define left
action momentum map $\mathbf{J}_{L}$ later.

$\displaystyle
T_{(g,v)(h,u)}L_{(g,v)^{-1}}(v_{gh},v+\zeta(g)u,w)=\left(T_{gh}L_{g^{-1}}(v_{gh},u,\zeta(g^{-1})w\right),$
(27)

for any $(v_{gh},v+\zeta(g)u,w)\in T_{(g,v)(h,u)}\left(G\oplus V\right)$. The
cotangent lift is

$\displaystyle\left(T_{(g,v)(h,u)}L_{(g,v)^{-1}}\right)^{\ast}\left(\alpha_{h},v,a\right)=\left(\left(T_{gh}L_{g^{-1}}\right)^{\ast}\alpha_{h},u+\zeta(g)v,\zeta_{\ast}(g)a\right),$
(28)

where $(\alpha_{h},v,a)\in T^{\ast}_{(h,u)}\left(G\oplus V\right)$. On the
other hand, the right action of $S$ on itself is

$\displaystyle R_{(g,v)}(h,u)=\left(hg,u+\zeta(h)v\right),$ (29)

and the tangent lift is

$\displaystyle
T_{(h,u)}R_{(g,v)}(v_{h},u,w)=\left(T_{h}L_{g}(v_{h},u+\zeta(h)v,\zeta(h)w\right),$
(30)

and the lift over the point $(g,v)(h,u)\in S$,

$\displaystyle T_{(h,u)(g,v)}R_{(g,v)^{-1}}$
$\displaystyle(v_{hg},u+\zeta(h)v,w)=$ (31)
$\displaystyle\left(T_{hg}R_{g^{-1}}(v_{hg},v,w-T_{hg}\zeta(v_{hg})\cdot\zeta(g^{-1})v\right),$

and the cotangent lift is

$\displaystyle\left(T_{(h,u)(g,v)}R_{(g,v)^{-1}}\right)^{\ast}$
$\displaystyle\left(\alpha_{h},v,a\right)=$ (32)
$\displaystyle\left(\left(T_{hg}R_{g^{-1}}\right)^{\ast}\alpha_{h}-\mathbf{d}f_{\zeta(g^{-1})v}(hg),v+\zeta(h)v,a\right)$

With that we can write the left momentum map,
$\mathbf{J}_{L}:T^{\ast}\left(G\oplus V\right)\to\mathfrak{s}^{\ast}_{+}$ and
the right momentum map, $\mathbf{J}_{R}:T^{\ast}\left(G\oplus
V\right)\to\mathfrak{s}^{\ast}_{-}$ by using (28) and (32) with
$(g,v)=(g,v)^{-1}$ and $(h,u)=(g,v)$ and that gives

$\displaystyle\mathbf{J}_{L}$
$\displaystyle=\left(T_{(e,0)}R_{(g,v)}\right)^{\ast}\left(\alpha_{g},u,a\right)=\left(\left(T_{e}R_{g}\right)^{\ast}\alpha_{g}+(\zeta^{\prime}_{u})^{\ast}a,a\right),$
(33) $\displaystyle\mathbf{J}_{R}$
$\displaystyle=\left(T_{(e,0)}L_{(g,v)}\right)^{\ast}\left(\alpha_{g},u,a\right)=\left(\left(T_{e}L_{g}\right)^{\ast}\alpha_{g},\zeta^{\ast}(g)a\right).$
(34)

######  Theorem 3.1.

Given the left $G$-invariant Hamiltonian, $H$, then the left reduced
Lie–Poisson Hamilton–Jacobi equation for semidirect products for a function
$S_{L}:G\oplus V\times\mathbb{R}\to\mathbb{R}$ is the following:

$\displaystyle\partial_{t}S_{L}((g,v),t)$
$\displaystyle+H_{r}\left(TR^{\ast}_{g}\mathbf{d}_{g}S_{L}+(\zeta^{\prime}_{v})^{\ast}\mathbf{d}_{v}S_{L},\mathbf{d}_{v}S_{L}\right)=0$
(35)

The solution of equation (35), known as Lie–Poisson Hamilton–Jacobi equation,
generates a Poisson map that defines the flow of the reduced Hamiltonian
$H_{r}$. The Poisson map transforms the initial momentum map
$\mu_{0}\to\mu(t)$, where

$\displaystyle(\mu_{0},\Gamma_{0})=$
$\displaystyle\left(-TL_{g}^{\ast}\mathbf{d}_{g}S_{L},\zeta^{\ast}(g)\mathbf{d}_{v}S_{L}\right),$
(36) $\displaystyle\text{and}\quad(\mu(t),\Gamma(t))=$
$\displaystyle\left(TR^{\ast}_{g}\mathbf{d}_{g}S_{L}+(\zeta^{\prime}_{v})^{\ast}\mathbf{d}_{v}S_{L},\mathbf{d}_{v}S_{L}\right),$
(37)

where $g\in G$ and $v\in V$.

For construction of a numerical scheme, the steps followed are similar to the
ones taken in section 4. First, the equation is approximated numerically in
time, using a numerical scheme, such as Euler method or Runge–Kutta method,
while maintaining that is it is parameterised in $(g,v)$. For example, a first
order in $t$ approximation is given by

$\displaystyle S_{L}((g,v))\approx S_{0}((g,v))+\Delta
tH\left(TR^{\ast}_{g}\mathbf{d}_{g}S_{0}+(\zeta^{\prime}_{v})^{\ast}\mathbf{d}_{v}S_{0},\mathbf{d}_{v}S_{0}\right),$
(38)

where $S_{0}$ is the initial condition of $S_{L}$ and it is chosen such that
it satisfies the identity mapping. The next step is to then use this
approximation of $S_{L}$ and substitute it into equation (36) and it allows to
obtain the group action $(g,v)$ that would result in $(\mu_{0},\Gamma_{0})$.
Once the action $(g,v)$ is found, it is then used in (37) to compute
$(\mu(t),\Gamma(t))$.

####  3.1.1 Lie Algebra of Pseudodifferential Symbols

Here we give briefly go over the basics of Lie group whose Lie algebra are
pseudodifferential symbols. A profound study of such groups is not required at
the moment as we only require the basics. For detailed exposition on this
topic, the reader can refer to [12, 57, 56]. This Lie algebra is required as
it allows for the application of Lie–Poisson Hamilton–Jacobi to Lie–Poisson
partial differential equations.

We begin by defining what we mean by an operator of order $m$,

$\displaystyle P=\sum\limits_{|\alpha|\leq m}a_{\alpha}(x)D^{\alpha},$ (39)

where $D$ is a differential operator and $\alpha$ is a multiindex and
$a_{\alpha}\in H^{s}(\mathbb{R})$. The _symbol_ of the operator $P$ is a
polynomial in $\xi$ variable

$\displaystyle p(x,\xi)=\sum\limits_{|\alpha|\leq
m}a_{\alpha}(x)\xi^{\alpha}.$ (40)

In this section we denote the algebra of pseudodifferential symbols by
$\Psi\mathrm{DS}_{m}$, where the $m$ subscript indicates the order. The union
of such algebras $\Psi\mathrm{DS}_{m}$ forms an infinite dimensional algebra
$\Psi\mathrm{DS}=\bigcup\limits_{k}\Psi\mathrm{DS}_{k}$ and it is closed under
multiplication defined by

$\displaystyle U(x,\xi)\circ M(x,\xi)=\sum\limits_{n\geq
0}\frac{1}{n!}\left(\frac{d^{n}}{d\xi^{n}}U(x,\xi)\right)\left(\frac{d^{n}}{dx^{n}}M(x,\xi)\right),$
(41)

where $U\in\Psi\mathrm{DS}_{m}$ and $M\in\Psi\mathrm{DS}_{n}$. The
multiplication, $U(x,\xi)\circ M(x,\xi)$, belongs to the group
$\Psi\mathrm{DS}_{n+m}\subset\Psi\mathrm{DS}$, i.e. multiplication is a closed
group operation. The Lie bracket of the Lie algebra $\Psi\mathrm{DS}$ is given
by

$\displaystyle\left[U,M\right]$ $\displaystyle=U\circ M-M\circ U$ (42)
$\displaystyle=\sum\limits_{n\geq
0}\frac{1}{n!}\left\\{\left(\frac{d^{n}}{d\xi^{n}}U(x,\xi)\right)\left(\frac{d^{n}}{dx^{n}}M(x,\xi)\right)-\left(\frac{d^{n}}{d\xi^{n}}M(x,\xi)\right)\left(\frac{d^{n}}{dx^{n}}U(x,\xi)\right)\right\\}.$

The algebra $\Psi\mathrm{DS}$ can be written as the direct sum of two other
subalgebras, i.e.

$\displaystyle\Psi\mathrm{DS}=\Psi\mathrm{DS}_{\mathrm{DO}}\oplus\Psi\mathrm{DS}_{\mathrm{INT}},$
(43)

where

$\displaystyle\Psi\mathrm{DS}_{\mathrm{DO}_{m}}:=\left\\{\sum\limits_{\alpha=0}^{m}a_{\alpha}(x)\partial^{\alpha}\right\\},\quad\text{and}\
\quad\Psi\mathrm{DS}_{\mathrm{INT}_{m}}:=\left\\{\sum\limits_{\alpha=-m}^{-1}u_{\alpha}(x)\partial^{\alpha}\right\\},$

and $\Psi\mathrm{DS}_{\mathrm{DO}}=\bigcup\limits_{m\geq
0}\Psi\mathrm{DS}_{\mathrm{DO}_{m}}$ and
$\Psi\mathrm{DS}_{\mathrm{INT}}=\bigcup\limits_{m\geq
0}\Psi\mathrm{DS}_{\mathrm{INT}_{m}}$.

Usually for Lie algebras, there’s an associated Lie group and to map the Lie
algebra to the Lie group an exponential map is used. The map needs to be
bijective and here is defined as the solution of the following differential
equation

$\displaystyle\frac{d}{ds}G(x,s)=U\circ G(s),$ (44)

where $G$ is the group element and $U$ belongs to the Lie algebra, and the
initial condition is chosen as $G(0)=\mathrm{I_{d}}$. The differential
equation (44) actually yields a system of differential equations, in most
cases it is a triangular system, for each coefficients of symbols of $G$ and
solving this system is required to obtain the coefficients of the group
element $G$ [57]. Here we use the $\Psi G$ as the group whose vector fields at
the identity are the $\Psi\mathrm{DS}$ symbols. Basically, $\Psi G$ is a group
whose elements are solutions of (44). The reasoning behind not considering it
a full Lie group is that the full algebra of symbols $\Psi\mathrm{DS}$ does
not have a corresponding Lie group nor its subalgebra
$\Psi\mathrm{DS}_{\mathrm{DO}}$. However, for the subalgebra
$\Psi\mathrm{DS}_{\mathrm{INT}}$ does have an associated Lie group $G_{INT}$
and the exponential map is bijective
$\mathrm{exp}:\mathfrak{g}_{\mathrm{INT}}\to G_{\mathrm{INT}}$. We denote the
algebra $\Psi\mathrm{DS}_{\mathrm{INT}}$ by $\mathfrak{g}_{\mathrm{INT}}$
henceforth. The pairing between the Lie algebra $\mathfrak{g}_{\mathrm{INT}}$
and its dual $\mathfrak{g}_{\mathrm{INT}}^{\ast}$,
$\left\langle\cdot,\cdot\right\rangle:\mathfrak{g}_{\mathrm{INT}}\times\mathfrak{g}_{\mathrm{INT}}^{\ast}\to\mathbb{R}$
is defined by

$\displaystyle\left\langle M,U\right\rangle=\mathrm{Tr}\left(M\circ
U\right)=\int\mathrm{res}\left(M\circ U\right),\quad
U\in\mathfrak{g}_{\mathrm{INT}},\,M\in\mathfrak{g}_{\mathrm{INT}}^{\ast}.$
(45)

and we actually find that
$\mathfrak{g}_{\mathrm{INT}}^{\ast}\simeq\mathfrak{g}_{\mathrm{DO}}$, where
$\mathfrak{g}_{\mathrm{DO}}$ denotes the Lie algebra of differential operators
$\Psi\mathrm{DS}_{\mathrm{DO}}$. For Hamiltonian systems, we need to compute
functional derivatives. Let $M$ be an element of
$\mathfrak{g}^{\ast}_{\mathrm{INT}}$, and consider the functional
$F:\mathfrak{g}_{\mathrm{INT}}\to\mathbb{R}$, its functional derivative is
obtained using the following pairing

$\displaystyle DF(M)\cdot M=\left\langle\frac{\delta F}{\delta
M},M\right\rangle=\sum\limits_{i}\int\frac{\delta F}{\delta m_{i}}m_{i}\,dx,$
(46)

thus we obtain,

$\displaystyle\frac{\delta F}{\delta M}=\sum\limits_{i=1}^{m}\frac{\delta
F}{\delta m_{i}}\xi^{-i}.$ (47)

Using the functional derivatives, one can define the Lie–Poisson bracket by

$\displaystyle\left\\{F,H\right\\}\left(M\right)=\left\langle
M,\left[\frac{\delta F}{\delta M},\frac{\delta H}{\delta
M}\right]\right\rangle=\int\mathrm{res}\left(M\circ\left[\frac{\delta
F}{\delta M},\frac{\delta H}{\delta M}\right]\right),$ (48)

where $H$ and $F$ are functionals that depend on the coefficients of the
pseudodifferential symbols.

To apply theorem 3.1, for deep learning, we first need to compute the
cotangent lifts of the action on the dual of the Lie algebra. Let $A$ be an
element of the group $G_{INT}$. The multiplication operator is defined by (41)
and to compute the tangent lift, let $B\in\Psi G_{INT}$ such that
$\left.\frac{d}{ds}\right|_{s=0}B(s)=X$, where
$X\in\mathfrak{g}_{\mathrm{INT}}$. Multiplying $A$ with $B$ and taking the
derivative with respect to $s$ and evaluate at $s=0$, we have

$\displaystyle TL_{A}X=\left.\frac{d}{ds}\right|_{s=0}A(t)\circ B(s)=A(t)\circ
X,$ (49)

and

$\displaystyle TR_{A}X=\left.\frac{d}{ds}\right|_{s=0}B(s)\circ A(t)=X\circ
A(t).$ (50)

To compute the cotangent lifts, $TL^{\ast}_{A}$ and $TR^{\ast}_{A}$, are
computed by pairing the tangent lifts to an element
$\alpha\in\mathfrak{g}_{\mathrm{INT}}^{\ast}$, the pairing is done using the
trace operator, i.e

$\displaystyle\mathrm{Tr}\left(TL^{\ast}_{A}\alpha\circ X\right)$
$\displaystyle=\mathrm{Tr}\left(\alpha\circ
TL_{A}X\right)=\int\mathrm{res}\left(\alpha\circ A(t)\circ X\right),$ (51)

and

$\displaystyle\mathrm{Tr}\left(TR^{\ast}_{A}\alpha\circ X\right)$
$\displaystyle=\mathrm{Tr}\left(\alpha\circ
TR_{A}X\right)=\int\mathrm{res}\left(\alpha\circ X\circ A(t)\right).$ (52)

These yield the following cotangent lift of the left and right action

$\displaystyle TL^{\ast}_{A}\alpha$ $\displaystyle=\alpha\circ A(t),$ (53)
$\displaystyle TR^{\ast}_{A}\alpha$ $\displaystyle=A(t)\circ\alpha.$ (54)

###  3.2 Hydrodynamics Formulation of Deep Learning

From the basic iterative relation of deep neural networks, we can interpret
them, in the language of geometric mechanics, as applying a pullback to a
curve on a manifold. It is very reminiscent of the particle relabelling
symmetry of fluids labels, i.e the action of the group of diffeomorphism,
denoted by $\mathrm{Diff}(Q)$, on the Lagrangian labels [7, 28]. The vector
field associated with $\mathrm{Diff}(Q)$ is denoted by $\mathfrak{X}(Q)$. The
neurons of the hidden layers can be viewed as Lagrangian labels

$\displaystyle\frac{d}{dt}q(x,t)=\sigma(u\circ q(x,t)),$ (55)

where $q\in\mathrm{Diff}(Q)$ and $u\in\mathfrak{X}(Q)$, and each time they
propagate through to the next layer, a transformation is applied to the
labels. Let $(\Omega,\mathcal{F},\mathbb{P})$ be the probability space, with
$\Omega$ the sample space, and $\mathcal{F}$ is $\sigma$-algebra of measurable
events of $\Omega$ and $\mathbb{P}$ is the probability measure on
$(\Omega,\mathcal{F})$. We consider the stochastic differential equation for
continuous ResNet

$\displaystyle d_{t}q(t)=(u\circ q)\,dt+\nu dW_{t},\quad
q(0)=q_{0},\quad\text{and}\quad q(T)=q_{T},$ (56)

where $u\in\mathfrak{X}(Q)$ is a vector field that plays the role of network
parameters, and $\xi_{k}\in\mathfrak{X}(Q)$ are vector fields that couple the
Brownian motion $W_{t}^{k}$ to $\mathfrak{X}(Q)$ and $\nu\in\mathbb{R}$.
Notice in (56) the activation function was removed because it is redundant in
this form. Here, the network parameters are specified to be a vector field on
the group of diffeomorphisms on $Q$ and can be viewed as a combination of the
weight matrix and activation function. Since the system is stochastic, the
cost functional (57) cannot be used in its current form, instead, we use its
expected value. Thus, we can state the training problem of stochastic
continuous ResNet as follows:

######  Problem 3.2.

Let $(\Omega,\mathcal{F},\mathbb{P})$ and given the training set
$\\{(q_{0}^{(i)},c^{(i)})\\}_{i=0,\dots,N}$, compute $u\in$ such that it
minimises

$\displaystyle
S=\int_{Q}\int_{0}^{T}\rho(x,t)\ell_{d}\left(u(t)\right)\,dt\,dx+\sum_{i=0}^{N}\int_{Q}\rho(x,t)L(\pi((q(T)^{(i)}),c^{(i)})\,dx,$
(57)

where $\rho(x,t)$ is the probability density function. The solution should
make $q(t)$ satisfy

$\displaystyle d_{t}q=u\circ q\,dt+\nu\sum\limits_{k}\xi_{k}dW_{t}^{k},$ (58)

and $q(0)=q_{0}^{(i)}\,\text{and}\,\pi(q(T))=c^{(i)}$ for $i=0,\dots,N$.

The solution method for this optimal control problem is mean–field optimal
control, described in [11]. Instead of controlling (58), the equation for the
evolution of probability density function associated with $q(t)$, known as the
Fokker–Planck equation,

$\displaystyle\partial_{t}\rho(x,t)+\nabla\cdot\left(u(x,t)\rho(x,t)\right)-\frac{\nu^{2}}{2}\nabla\cdot(\nabla\rho(x,t))=0,\quad\rho(x,0)=\rho_{0}(x).$
(59)

This way, we can apply Pontryagin’s Maximum Principle in a similar way to the
deterministic optimal control problems. The result is the mean–field games
system

$\displaystyle\rho(x,t)\frac{\delta\ell_{d}}{\delta
u}-\rho(x,t)\nabla\lambda(x,t)=0$ (60)
$\displaystyle\partial_{t}\rho(x,t)+\nabla\cdot\left(u(x,t)\rho(x,t)\right)-\frac{\nu^{2}}{2}\nabla\cdot(\nabla\rho(x,t))=0$
(61)
$\displaystyle\partial_{t}\lambda(x,t)+u(x,t)\nabla\lambda(x,t)+\frac{\nu^{2}}{2}\nabla\cdot(\nabla\lambda(x,t))=\ell_{d}$
(62)

where it is basically a coupled system of Fokker–Planck equation, (61), and
stochastic Hamilton–Jacobi–Bellman equation, (62). In [32], the network
parameter $u$ is substituted by
$u=\omega+\frac{\nu^{2}}{2}\nabla\mathrm{log}(\rho)$ and this is done to
reduce the system of equations (60), (61) and (62) to the following equation

$\displaystyle\partial_{t}\frac{\delta\ell_{d}}{\delta\omega}+\mathrm{ad}^{\ast}_{\omega}\frac{\delta\ell_{d}}{\delta\omega}=\rho\nabla\frac{\delta\ell_{d}}{\delta\rho}$
(63) $\displaystyle\partial_{t}\rho+\nabla\cdot\left(\omega\rho\right)=0$

where
$\mathrm{ad}^{\ast}:\mathfrak{X}\times\mathfrak{X}^{\ast}\to\mathfrak{X}^{\ast}$
is the coadjoint map and in this case it is defined in coordinates of
$\mathbb{R}^{N}$ as

$\displaystyle\left(\mathrm{ad}_{\omega}^{\ast}\frac{\delta\ell_{d}}{\delta\omega}\right)_{i}=\left(\partial_{x^{j}}\left(\left(\frac{\delta\ell_{d}}{\delta\omega}\right)_{i}\omega^{j}\right)+\left(\frac{\delta\ell_{d}}{\delta\omega}\right)_{j}\partial_{x^{i}}\omega^{j}\right).$
(64)

One thing to note that the approach of [32] relies on the Lagrangian
formulation of mechanics and we cannot apply Hamilton–Jacobi theory directly,
but in order to do so, we need to translate it into the Hamiltonian point of
view.

######  Remark 3.3.

It is worth pointing out that when $\nu\to 0$ in (63), it is in effect the
celebrated hydrodynamics formulation of the Monge–Kantorovitch mass problem of
optimal transport [9]. And for the optimal control problem, with $\nu\to 0$ in
(60), (61) and (62), it was shown that it reduces to a Euler’s fluid equation
in [14].

###  3.3 Deep Learning as a Lie-Poisson System

In this section, we formulate the reduced mean–field games (63) to a
Lie–Poisson system, and in doing so, we translate it into the language of
Lie–Poisson Hamilton–Jacobi theory. This allows us to construct a Poisson
mapping from the target of the training dataset to the initial conditions.

This first problem that we need to overcome is that the Hamiltonian for the
mean–field games (60), (61) and (62) is a covariant Hamiltonian and theorem
3.1 is not adapted for such Hamiltonians. For covariant Hamiltonians, there
has been attempt to introduce a general Hamilton–Jacobi theory, such as the
work done in [47, 80, 101]. An alternative approach in [100], where boundary
integrals of the variational principle are used to derive generating functions
for partial differential equations. Having said that, the route taken here is
to expand the work with a larger configuration manifold, that is, the group of
pseudospectral symbols instead of the group diffeomorphism. In order to use
Lie–Poisson Hamilton–Jacobi theory, we discretise the group of pseudospectral
symbols in space only. This is inspired by the work of [87] where the group of
diffeomorphisms is discretised as $GL(N)$ and EPDiff equation is written in
terms of vector fields on dual of the discretised Lie algebra which is a
subalgebra of $\mathfrak{gl}(N)$.

For the deep learning problem, the Hamiltonian for the system (63) is

$\displaystyle H$
$\displaystyle=\frac{1}{2}\int_{Q}\rho\|\omega\|^{2}d\mathrm{Vol}(x)+\frac{\nu^{2}}{2}\int_{Q}\rho\left\langle\omega,\nabla\mathrm{log}(\rho)\right\rangle
d\mathrm{Vol}(x)$ (65)
$\displaystyle+\frac{\nu^{4}}{8}\int_{Q}\rho\|\nabla\mathrm{log}(\rho)\|^{2}d\mathrm{Vol}(x)$

and $H$ is a function defined on the dual of the Lie algebra of a semidirect
product
$\mathfrak{s}^{\ast}=\mathfrak{g}_{\mathrm{INT}}^{\ast}\oplus\mathrm{Dens}(Q)$.
The corresponding group is $S=G_{\mathrm{INT}}\oplus\mathrm{Dens}(Q)$. From
the definitions of group actions on itself, tangent and cotangent lifts
defined in section 3.1, here $G_{\mathrm{INT}}$ acts on $\mathrm{Dens}(Q)$ by
(push-forward), i.e.

$\displaystyle\zeta(g)\rho=(g)_{\ast}\rho,\quad g\in G_{\mathrm{INT}}$ (66)

Another important map that we need to define is the algebra
$\mathfrak{g}_{\mathrm{INT}}$ representation on $\mathrm{Dens}(Q)$, denoted by
$\zeta^{\prime}:\mathfrak{g}_{\mathrm{INT}}\to\mathrm{End}\left(\mathrm{Dens}(Q)\right)$,
is given by

$\displaystyle\zeta^{\prime}(U)\rho=-\pounds_{U}\rho,$ (67)

where $U\in\mathfrak{g}_{\mathrm{INT}}$ and $\pounds_{U}$ denotes the Lie
derivative along the direction of $U$, and in coordinates
$\zeta^{\prime}(U)\rho=\mathrm{div}\left(\rho U\right)$. Once we have the
representation on the Lie algebra, we derive the representation of the dual of
the Lie algebra $\mathfrak{g}_{\mathrm{INT}}^{\ast}$, i.e. the mapping
$\left(\zeta^{\prime}_{U}\right)^{\ast}:\mathrm{Dens}(Q)^{\ast}\to\mathfrak{g}^{\ast}_{\mathrm{INT}}$.
It is obtained by using the pairing and for any $\rho\in\mathrm{Dens}(Q)$,
$\eta\in\mathrm{Dens}(Q)^{\ast}$ and $U\in\mathfrak{g}_{\mathrm{INT}}$ we have

$\displaystyle\left\langle\eta\diamond\rho,U\right\rangle$
$\displaystyle=\left\langle\left(\zeta^{\prime}_{\eta}\right)^{\ast}\rho,U\right\rangle=-\left\langle\eta,\zeta_{U}\rho\right\rangle$
(68) $\displaystyle=-\int\eta\nabla\cdot(U\rho)\,dx=-\int\rho\nabla\eta\cdot
U\,dx,$

and this means
$\eta\diamond\rho=\left(\zeta^{\prime}_{\eta}\right)^{\ast}\rho=\rho\nabla\eta$.
The operator
$\diamond:\mathrm{Dens}(Q)\times\mathrm{Dens}(Q)^{\ast}\to\mathfrak{g}_{\mathrm{INT}}^{\ast}$
is often used in describing motion of various hydrodynamics, introduced in
[46], and it used to map Clebsch variables to the dual of Lie algebra [44].

Going back to training of deep neural networks on Poisson manifolds and for
this particular problem we choose the symbol to be

$\displaystyle U(x,\xi)=u(x)\,\xi^{-1}\in\mathfrak{g}_{\mathrm{INT}},$ (69)

and we require that the dual space $\mathfrak{g}_{\mathrm{INT}}^{\ast}$ should
be the space of pseudodifferential symbols of order $2$, i.e.

$\displaystyle M(x,\xi)=\sum_{i=0}^{2}m_{i}(x)\xi^{i}$ (70)

The importance of these conditions on the Lie algebra is that for the
functionals $F$ and $H$ substituting them into the Lie–Poisson bracket (48),
we obtain

$\displaystyle\left\\{F,H\right\\}=\int m(x)\left[\frac{\delta F}{\delta
m}\nabla\left(\frac{\delta H}{\delta m}\right)-\frac{\delta H}{\delta
m}\nabla\left(\frac{\delta F}{\delta m}\right)\right]\,dx.$ (71)

The Lie–Poisson bracket of functionals on the dual of the Lie algebra of
semidirect product
$\mathfrak{s}^{\ast}=\mathfrak{g}_{\mathrm{INT}}^{\ast}\oplus\mathrm{Dens}(Q)^{\ast}$
is

$\displaystyle\\{F,H\\}_{\pm}(M,\rho)=\pm\left\langle M,\left[\frac{\delta
F}{\delta M},\frac{\delta H}{\delta
M}\right]\right\rangle\mp\left\langle\rho,\pounds_{\frac{\delta F}{\delta
M}}\frac{\delta H}{\delta\rho}-\pounds_{\frac{\delta H}{\delta M}}\frac{\delta
F}{\delta\rho}\right\rangle,$ (72)

where $F,H:\mathfrak{s}^{\ast}\to\mathbb{R}$ and the first term of the
Lie–Poisson bracket is given by (48). In explicit terms,

$\displaystyle\\{F,H\\}_{\pm}(M,\rho)$
$\displaystyle=\pm\int_{Q}m(x)\left[\frac{\delta F}{\delta
m}\nabla\left(\frac{\delta H}{\delta m}\right)-\frac{\delta H}{\delta
m}\nabla\left(\frac{\delta F}{\delta m}\right)\right]\,dx$ (73)
$\displaystyle\mp\int_{Q}\rho\left[\left(\frac{\delta H}{\delta
m}\cdot\nabla\right)\frac{\delta F}{\delta\rho}-\left(\frac{\delta F}{\delta
m}\cdot\nabla\right)\frac{\delta H}{\delta\rho}\right]\,dx.$

######  Remark 3.4.

This Lie–Poisson bracket is the same bracket for ideal compressible fluid [75,
76].

Since the Lie–Poisson bracket is in terms of $(M,\rho)\in\mathfrak{s}^{\ast}$,
and that necessitates the need to write the Hamiltonian (65) in those
coordinates and for that we use $M=\rho\omega$,

$\displaystyle H(M,\rho)$
$\displaystyle=\frac{1}{2}\int_{Q}\frac{\|M\|^{2}}{\rho}d\mathrm{Vol}(x)+\frac{\nu^{2}}{2}\int_{Q}\left\langle
M,\nabla\mathrm{log}(\rho)\right\rangle d\mathrm{Vol}(x)$ (74)
$\displaystyle+\frac{\nu^{4}}{8}\int_{Q}\rho\|\nabla\mathrm{log}(\rho)\|^{2}d\mathrm{Vol}(x)$

Plugging in $F=(M,\rho)$, we obtain the Lie–Poisson on semidirect product
space

$\displaystyle\partial_{t}M$ $\displaystyle=\\{M,H\\}_{+}$ (75)
$\displaystyle\partial_{t}\rho$ $\displaystyle=\\{\rho,H\\}_{+},$

or in coordinates od $\mathbb{R}^{N}$, it is written as

$\displaystyle\partial_{t}m_{i}$
$\displaystyle=m_{i}\partial_{x^{i}}\frac{\delta H}{\delta
m_{i}}+\partial_{x^{i}}\left(m_{i}\frac{\delta H}{\delta
m_{i}}\right)+\rho\partial_{x^{i}}\frac{\delta H}{\delta\rho},$
$\displaystyle\text{for}\,i=1,\dots,N$ (76) $\displaystyle\partial_{t}\rho$
$\displaystyle=-\partial_{x^{j}}\left(\rho\frac{\delta H}{\delta
m_{j}}\right),$

where
$\mathrm{ad}^{\ast}:\mathfrak{g}_{\mathrm{INT}}\times\mathfrak{g}_{\mathrm{INT}}^{\ast}\to\mathfrak{g}_{\mathrm{INT}}^{\ast}$
and it is the same as in (64).

######  Proposition 3.5.

The deep learning dynamics can be written in Lie–Poisson form

$\displaystyle\partial_{t}M$
$\displaystyle=\pm\mathrm{ad}^{\ast}_{\frac{\delta H}{\delta
M}}M\mp\frac{\delta H}{\delta\rho}\diamond\rho$ (77)
$\displaystyle\partial_{t}\rho$ $\displaystyle=\mp\pounds_{\frac{\delta
H}{\delta M}}\rho,$

where
$\mathrm{ad}^{\ast}:\mathfrak{g}_{\mathrm{INT}}\times\mathfrak{g}_{\mathrm{INT}}^{\ast}\to\mathfrak{g}_{\mathrm{INT}}^{\ast}$
and it is the same as in (64).

###### Proof.

Let $(\eta,\varrho),(u,a)\in\mathfrak{s}$ and
$(M,\rho)\in\mathfrak{s}^{\ast}$, from the definition of the coadjoint map and
the diamond operator, we then have

$\displaystyle\left\langle\left(\mathrm{ad}^{\ast}_{u}M-a\diamond\rho,-\pounds_{u}\rho\right),\left(\eta,\varrho\right)\right\rangle$
$\displaystyle=\left\langle\mathrm{ad}^{\ast}_{u}M-a\diamond\rho,\eta\right\rangle+\left\langle-\pounds_{u}\rho,\varrho\right\rangle$
$\displaystyle=\left\langle
M,\mathrm{ad}_{u}\eta\right\rangle+\left\langle\rho,\pounds_{u}\varrho-\pounds_{\eta}a\right\rangle$
$\displaystyle=\left\langle
M,[u,\eta]\right\rangle+\left\langle\rho,\pounds_{u}\varrho-\pounds_{\eta}a\right\rangle$
$\displaystyle=\left\langle M,[u,\eta]\right\rangle-\left\langle\rho\nabla
a,\eta\right\rangle+\left\langle\rho,\pounds_{u}\varrho\right\rangle$
$\displaystyle=\left\langle\\{M,u\\},\eta\right\rangle+\left\langle\pounds_{u}\rho,\varrho\right\rangle$
$\displaystyle=\left\langle\left(\\{M,u\\},\\{\rho,u\\}\right),\left(\eta,\varrho\right)\right\rangle.$

From the equation we get

$\displaystyle\left(\mathrm{ad}^{\ast}_{u}M-a\diamond\rho,-\pounds_{u}\rho\right)=\left(\\{M,u\\},\\{\rho,u\\}\right)$

and thus we can write equation (75) as

$\displaystyle\partial_{t}M$
$\displaystyle=\pm\mathrm{ad}^{\ast}_{\frac{\delta H}{\delta
M}}M\mp\frac{\delta H}{\delta\rho}\diamond\rho$
$\displaystyle\partial_{t}\rho$ $\displaystyle=\mp\pounds_{\frac{\delta
H}{\delta M}}\rho.$

∎

Now, we can apply the Lie–Poisson Hamilton–Jacobi theory for semidirect
products as

$\displaystyle\partial_{t}S_{L}((A,\rho))$
$\displaystyle+h\left(TR^{\ast}_{A}\mathbf{d}_{A}S_{L}+(\zeta^{\prime}_{\rho})^{\ast}\mathbf{d}_{\rho}S_{L},\mathbf{d}_{\rho}S_{L}\right)=0$
(78)

The solution of equation (35) generates a Poisson map that defines the flow of
the reduced Hamiltonian $h$. The Poisson map transforms the initial momentum
map $\mu_{0}\to\mu(t)$, where

$\displaystyle(\mu_{0},\rho_{0})=$ $\displaystyle(-\mathbf{d}_{A}S_{L}\circ
A,\zeta(A)^{\ast}\mathbf{d}_{\rho}S_{L}),$ (79)
$\displaystyle\text{and}\quad(\mu(t),\rho(t))=$
$\displaystyle\left(A\circ\mathbf{d}_{A}S_{L}+(\zeta^{\prime}_{\rho})^{\ast}\mathbf{d}_{\rho}S_{L},\mathbf{d}_{\rho}S_{L}\right).$

Now we decipher the two previous equations in (79). The first component
$\mathbf{d}_{A}S_{L}\circ A$ is the multiplication in $\Psi\mathrm{DS}$ and
the second, as mentioned, is the action of $\Psi G$ acting on
$\mathrm{Dens}(Q)$ by pullback. The term
$(\zeta^{\prime}_{u})^{\ast}\mathbf{d}_{u}S_{L}$ is the diamond operator
defined in (68). The terms $(\mathbf{d}_{A}S_{L},\mathbf{d}_{u}S_{L})$ can be
eliminated from the second equation in (79), by taking the left action of
$(A,u)^{-1}$ of $(\mu_{0},\rho_{0})$. This results in

$\displaystyle(\mu(t),\rho(t))$ $\displaystyle=(-A\circ\mu_{0}\circ
A^{-1}+(\zeta^{\prime}_{\rho})^{\ast}\rho_{0}(A\circ x),\rho_{0}(A\circ x))$
(80) $\displaystyle=\mathrm{Ad}_{(A,u)^{-1}}^{\ast}(\mu_{0},\rho_{0}),$

where
$\mathrm{Ad}_{(A,u)^{-1}}^{\ast}:S\times\mathfrak{s}^{\ast}\to\mathfrak{s}^{\ast}$
is the coadjoint operator.

Despite the Lie–Poisson Hamilton–Jacobi theory introduces a level of
complexity to the matter; the power of the method comes out when it is
employed for constructing numerical algorithms. Up until now, we have worked
exclusively with continuous-time representation of deep learning, but for it
to applicable we need to discretise it. The discrete Lie–Poisson
Hamilton–Jacobi theory would yield an $N$ dimensional system that preserves
the geometric structure of the continuous equations. Furthermore, the flow or
solution trajectories, generated by the Hamiltonian vector-field corresponding
to (73) belong to the coadjoint orbit. _The coadjoint orbit_ ,
$\mathcal{O}_{(\mu,\rho)}$, is the set of elements of $\mathfrak{s}^{\ast}$
that are moved by the elements of the group $S=\Psi G\times\mathrm{Dens}(Q)$,
or mathematically,

$\displaystyle\mathcal{O}_{(\mu,\rho)}:=\left\\{\mathrm{Ad}^{\ast}_{(A,u)^{-1}}(\mu,\rho)\,|\,(A,u)\in\Psi
G\right\\},$ (81)

and $(\mu,\rho)\in\mathfrak{s}^{\ast}$. The coadjoint orbit is in fact a sub-
manifold of $\mathfrak{s}^{\ast}$ and it is endowed with a symplectic
structure [73]. The Lie–Poisson Hamilton–Jacobi constructs a numerical
solution that naively belongs to $\mathcal{O}_{(\mu,\rho)}$ and this makes it
worth studying. The preservation of the coadjoint orbit was further looked at
in [74].

What is interesting regarding the discrete equations is that they are also a
Hamiltonian system and this is a consequence of applying reduction to
generating function on cotangent bundle. In [33] it was shown how to derive
the modified Hamiltonian for the discrete equations obtained from Lie–Poisson
Hamilton–Jacobi theory, albeit the systems considered in [33] are stochastic.
However, the same steps can be followed for the deterministic case.

The fact that mean–field games system (60), (61) and (62) is written as a
coadjoint motion, then according to theorem $2.7$ in [35], the optimal
solution is restricted to $\mathcal{O}_{(\mu,\rho)}$. This means
$(\mu_{0},\rho_{0})$ as well as $(\mu_{T},\rho_{T})$ are on the same coadjoint
orbit $\mathcal{O}_{(\mu,\rho)}$ and that means the curve
$(\mu(t),\rho(t))\in\mathfrak{s}^{\ast}$ is a minimising curve of the
functional and it also connects the two densities $\rho_{0}$ and $\rho_{T}$.

##  4 Hamilton–Jacobi Integrators for Deep Learning

In this section we show the application of generating function in constructing
a structure–preserving integrator for training deep neural network. Symplectic
transformations, that are generated by the solution of the Hamilton–Jacobi
equation, preserve the structure and it is the basis for construction of the
numerical scheme. This is done by applying a numerical method for solving the
Hamilton–Jacobi equation and then use that with (7) to obtain the approximated
solution. From Pontryagin’s Maximum Principle, the necessary conditions for
optimality of the solution of an optimal control problem satisfy a Hamiltonian
system (12). Thus, it is fitting that a geometric integrator is used to solve
optimal control problem to guarantee that the numerical solution is
representative of the control system. First, we describe how optimal is used
to formulate the training of deep neural network in section 2.3 and, here, we
adhere to the approach of [67, 10, 20]. The resultant Hamilton’s equations are
then solved numerically using structure–preserving integrator.

The idea of using generating functions for numerical integrators was first
introduced in [52, 21]. Recall that the solution of a Hamiltonian system is a
symplectomorphism and one time-step integration is also a symplectic
transformation $(q^{k},p^{k})\to(q^{k+1},p^{k+1})$ generated by the solution
of Hamilton–Jacobi equation at time $\Delta t$. Using a generating function of
type-II and (7), we obtain the following numerical scheme

$\displaystyle q^{k+1}=q^{k}+\frac{\partial S_{2}}{\partial
p}(q^{k},p^{k+1},\Delta t),\quad p^{k+1}=p^{k}-\frac{\partial S_{2}}{\partial
q}(q^{k},p^{k+1},\Delta t)$ (82)

and such scheme preserves the symplectic form
$\omega=\mathbf{d}q^{k}\wedge\mathbf{d}p^{k}=\mathbf{d}q^{k+1}\wedge\mathbf{d}p^{k+1}$.
The generating function, $S_{2}$ is the approximation of solution of
Hamilton–Jacobi equation and it can be expressed as a truncated power series

$\displaystyle S_{2}(q,p,t)=\sum\limits_{i=1}^{m}t^{i}S_{2,i}(q,p,t),$ (83)

where the coefficients of the power series are functions in term of the
Hamiltonian $H$ and its partial derivatives

$\displaystyle S_{2,1}=H(q,p),\quad
S_{2,2}=\frac{1}{2}\left\langle\frac{\partial H}{\partial q},\frac{\partial
H}{\partial p}\right\rangle,\quad S_{2},3=\dots$

The order of scheme (82) is determined by the number of terms, $m$, in the
approximation (83).

######  Remark 4.1.

When $m=1$, i.e. $S=\Delta tH$, the scheme (82) becomes the symplectic Euler
scheme

$\displaystyle q^{k+1}=q^{k}+\Delta t\frac{\partial H}{\partial
p}(q^{k},p^{k+1},\Delta t),\quad p^{k+1}=p^{k}-\Delta t\frac{\partial
H}{\partial q}(q^{k},p^{k+1},\Delta t)$

The application of structure–preserving integrators based on generating
function can be used to approximate the solution of optimal control problem,
now we use that to solve problem 2.3. Using Pontryagin’s Maximum Principle,
the control Hamiltonian $H:T^{\ast}Q\times\mathcal{U}\to\mathbb{R}$ is

$\displaystyle
H(q_{0},\dots,q_{N},p_{0},\dots,p_{N},\theta)=\sum\limits_{i=0}^{N}\left\langle
p_{i},\sigma(q_{i},\theta)\right\rangle-\frac{\gamma}{2}\|\theta\|^{2}$ (84)

where $\gamma\in\mathbb{R}_{\geq 0}$ and the last term is the regularisation
term. We also require that

$\displaystyle\frac{\delta
H}{\delta\theta}=\frac{1}{\gamma}\sum\limits_{i=0}^{N}p_{i}\frac{\partial\sigma}{\partial
z}q_{i}^{T}-\theta=0$

and substitute the expression for $\theta$ into $H$, we obtain a Hamiltonian
in terms of canonical coordinates $(q,p)$ on cotangent bundle $T^{\ast}Q$.
Then the Hamiltonian, $H:T^{\ast}Q\to\mathbb{R}$, becomes a function defined
on the cotangent space rather than the Pontryagin’s bundle and it is given by

$\displaystyle H(q_{0},\dots,q_{N},p_{0},\dots,p_{N},t)$
$\displaystyle=\sum\limits_{i=0}^{N}\left\langle
p_{i},\sigma\left(q_{i},\frac{1}{\gamma}\sum\limits_{j=0}^{N}p_{j}\frac{\partial\sigma}{\partial
z}q_{j}^{T}\right)\right\rangle$ (85)
$\displaystyle-\frac{1}{2\gamma}\left\|\sum\limits_{j=0}^{N}p_{j}\frac{\partial\sigma}{\partial
z}q_{j}^{T}\right\|^{2}-\frac{1}{2\gamma}\left\|\sum\limits_{j=0}^{N}p_{j}\frac{\partial\sigma}{\partial
z}\right\|^{2}$

(a)

(b)

(c)

(d)

Figure 2: The accuracy of the neural networks predictions is shown in figures
2(a) and 2(c), for the two spirals and two circles datasets respectively, over
$5000$ training iterations. On the other hand, the residual, computed by
$\sum\limits_{i=0}^{N}\|\pi(q^{(}50)_{i})-c^{i}\|^{2}$, is shown in figures
2(b) and 2(d), for the two spirals and two circles datasets respectively.

(a)

(b)

(c)

Figure 3: The transformation of two spirals training set
$\\{(q_{0}^{(i)},c^{(i)}_{N_{t}})\\}_{i=0,\dots,N}$ with $N=2000$ at the
$20^{th}$ layer (the second column from the left), $30^{th}$ layer (the third
column from the left), and the final layer (rightmost column). Figure 3(a),
shows the network discretised using Euler scheme, while figure 3(b), shows the
results obtained using $4^{th}$ order Runge–Kutta method . Finally, figure
3(c) is when the symplectic method using Euler scheme for the generating
function. The solid black line separates between the two classes/categories.

(a)

(b)

(c)

Figure 4: The trajectories of the two spirals data are plotted, where the
solid dots are the initial data and the hollow circles are the final data. The
light coloured lines are the individual paths each data point followed when
iterated through the discretised continuous residual network. Figure 4(a)
shows the paths when the network is discretised using Euler scheme, while
figure 4(b) and 4(c) show the trajectories of the $4^{th}$ order Runge–Kutta
and symplectic integrator, respectively.

######  Problem 4.2.

Given the training set $\\{(q_{0}^{(i)},c^{(i)}_{N_{t}})\\}_{i=0,\dots,N}$,
compute the sequence $(q^{k},p^{k})$ such that it satisfies

$\displaystyle q_{i}^{k+1}$ $\displaystyle=q_{i}^{k}+\Delta
t\sigma\left(\frac{1}{\gamma}\sum\limits_{j=0}^{N}\frac{\partial\sigma}{\partial
z}p_{j}^{k+1}(q_{j}^{k})^{T}q_{i}^{k}+\frac{1}{\gamma}\sum\limits_{j=0}^{N}\frac{\partial\sigma}{\partial
z}p_{j}^{k+1}\right),$ (86) $\displaystyle p_{i}^{k+1}$
$\displaystyle=p_{i}^{n}-\frac{\Delta t}{\gamma}\frac{\partial\sigma}{\partial
z}p_{i}^{k+1}\left(\sum\limits_{j=0}^{N}\frac{\partial\sigma}{\partial
z}p_{j}^{k+1}(q_{j}^{k})^{T}\right)$

for $k=0,\dots,N_{l}-1$ and the boundary conditions $q^{0}=q_{0}^{(i)}$ and
$\pi(q^{N_{t}})=c^{(i)}$ for each $i=0,\dots,N$.

In order to iterate through (86), we require $p^{0}$, however at this point,
it is not available. Further, the value of $p^{0}$ must ensure that
$q^{0}=q_{0}^{(i)}$ and $\pi(q^{N_{t}})=c^{(i)}$. There are a number of ways
in which one can choose a suitable initial value for the costate equation and
one of them is shooting method [61]. The idea is to pick an arbitrary $p^{0}$
and then evaluate how different $\pi(q^{N_{t}})$ from $c^{(i)}$ and use this
error to modify $p^{0}$. The process is then is iterated to refine the choice
of $p^{0}$. An alternative approach is to use the fact that we can derive the
$p^{N_{l}}$ from the cost functional. Once we have that, we can use that to
define (86) as the Karush–Kuhn–Tucker condition for nonlinear programming. For
such constraint optimisation problems, iterative methods such as sequential
quadratic programming [41], or even using Newton’s iterative method.

To evaluate the numerical scheme presented in this section, we train the
network for two separate classification problems: two spirals, and two
concentric circles classification. That said, this is just a toy example as a
sanity check. Both datasets contain data from only two classes, that means for
all classification space $\mathcal{C}=\\{0,1\\}$. The two spirals dataset
contains $2000$ data points for training and $2000$ other data points for
testing. Whereas, the two concentric circles dataset uses, for each training
and testing, $1000$ data points. The number of layers used for both datasets
is $N_{t}=50$ and $\Delta t=0.075$. The training was done over $5000$
iterations.

In order to put the structure–preserving scheme into perspective, other
numerical schemes for ordinary differential equations are used for neural
networks: Euler method, and fourth order explicit Runge–Kutta method. The
methods are compared using the error function in (57),
$L=\sum\limits_{i=0}^{N}\|\pi((q(50))^{(i)})-c^{(i)}\|^{2}$, and the accuracy
at each training steps. The accuracy and residuals for the two datasets are
shown in figure 2. Figure 3 shows snap shots of the data, only for the two
spirals dataset, propagating through the layers of each discretisation of the
continuous ResNet. At the final layer, we see that, for all implementations,
the majority of data points of each class are carried in a way they become
sorted. Figure 4, shows the trajectories that each data point takes for all
three discretisations.

###  4.1 Deep Neural Networks as Discrete Lie-Poisson System

Before discretising (79), we first need to review a particular discretisation
compatible with the Lie–Poisson Hamilton–Jacobi theory. The discretisation is
done on the algebra of pseudodifferential symbols $\Psi\mathrm{DS}$ instead of
discretising the Hamilton–Jacobi equation. The philosophy here is similar to
[84, 87], where the volume preserving group of diffeomorphisms is discretised
and approximated as a subgroup of $GL(N)$ Lie group instead.

Let $Q$ denote a smooth manifold, its discrete analogue is denoted by

$\displaystyle\mathsf{Q}:=\left\\{(x_{1}^{i},x_{2}^{i},\dots,x_{N}^{i})\,:\,i=0,\dots,N_{t},\quad(x_{1}^{i},x_{2}^{i},\dots,x_{N}^{i})\in
Q\right\\},$

and in this discretisation, the space between each point in each dimension is
a constant $\mathsf{D}_{x}$. Consider the pseudodifferential operator (39),
when applied to a function $u\in C^{\alpha}\left(Q\right)$, and instead of
evaluating the Fourier transform on $Q$, it is evaluated on $\mathsf{Q}$. With
out the loss of generality, in $1$ dimensions, the pseudodifferential operator
on $u$ can be written as

$\displaystyle P_{d}u(x_{j})=\frac{1}{(2\pi)^{n}}\sum\limits_{j,k}\int
e^{i(x_{k}-x_{j})\xi}a_{k}(x_{j})\xi^{k}u_{j}d\xi.$ (87)

From the previous expression, the following polynomial

$\displaystyle p(x_{j},\xi)=\sum\limits_{k\leq m}a_{k,j}\xi^{k},$ (88)

where $a_{k,j}=a_{k}(x_{j})$, is called _discrete symbol of order
$\mathbf{m}$_. We denote the algebra of pseudodifferential symbols of order
$m$ on $\mathsf{Q}$ by $\Psi\mathrm{DS}_{d,m}$. Basic operations and
definitions carry also to pseudodifferential operators on $\mathsf{Q}$.

######  Definition 4.3.

Let $A\in\Psi\mathrm{DS}_{d,m}$ and $B\in\Psi\mathrm{DS}_{d,n}$, then the
multiplication of $A\circ B$ on $\mathsf{Q}$ is the operator $C$ with symbol
of order $m+n$ defined by

$\displaystyle c_{j}(\xi)=$
$\displaystyle\sum\limits_{|\alpha|<N}\frac{1}{\alpha!}\partial^{\alpha}_{\xi}a_{j}(\xi)\mathsf{D}_{x}^{\alpha}{b}_{j}(\xi)$
(89)

where $\mathsf{D}_{x}^{\alpha}$ is the discrete difference operator in the $x$
direction of order $\alpha$.

######  Definition 4.4.

Let $A\in\Psi\mathrm{DS}_{d,m}$ and $B\in\Psi\mathrm{DS}_{d,n}$, the
commutator bracket for the $\Psi\mathrm{DS}_{d}$ is

$\displaystyle\left[A,B\right](x_{j})=\sum_{n\geq
1}\frac{1}{\alpha!}\partial^{\alpha}_{\xi}a_{j}(\xi)\mathsf{D}_{x}^{\alpha}{b}_{j}(\xi)-\mathsf{D}_{x}^{\alpha}{a}_{j}(\xi)\partial^{\alpha}_{\xi}b_{j}(\xi)$
(90)

Another important map that is used extensively is the trace map
$\mathrm{Tr}:\Psi\mathrm{DS}_{d}\to\mathbb{R}$, and here for
$A\in\Psi\mathrm{DS}_{d,m}$, it is defined by

$\displaystyle\mathrm{Tr}\left(A\right)=\sum\limits_{k}a_{-1,k}.$ (91)

Having defined the trace mapping, this enable us to define the pairing between
the algebra $\Psi\mathrm{DS}_{d}$ and its dual $\Psi\mathrm{DS}_{d}^{\ast}$
defined as a map
$\langle\cdot,\cdot\rangle:\Psi\mathrm{DS}_{d}\times\Psi\mathrm{DS}^{\ast}_{d}\to\mathbb{R}$
described by

$\displaystyle\left\langle A,B\right\rangle=\mathrm{Tr}(A\circ B),$ (92)

where $A\in\Psi\mathrm{DS}_{d}$, and $B\in\Psi\mathrm{DS}^{\ast}_{d}$.

######  Remark 4.5.

This is useful when the the space $\Psi\mathrm{DS}_{d}$ is decomposed into two
subspaces as in (43) to help identify the the Lie algebra with its dual.

To directly apply Lie–Poisson Hamilton–Jacobi theory of proposition 3.1, we
first need to give the explicit expression for the group action and its
cotangent lift on pseudodifferential algebra on $\mathsf{Q}$. As we are only
interested in a subalgebra of $\Psi\mathrm{DS}_{d}$, we only state the
algebraic expressions of the subalgebra of interest. From this point onwards,
we use the notation $\mathfrak{g}_{\mathrm{INT}}(\mathsf{Q})$ to denote the
algebra of pseudodifferential integral operators on the discrete set
$\mathsf{Q}$ and $\mathfrak{g}_{\mathrm{INT}}^{\ast}(\mathsf{Q})$ its dual.
Analogous to the continuous case, the Lie algebra
$\mathfrak{g}_{\mathrm{INT}}(\mathsf{Q})$ has an associated Lie group
$G_{\mathrm{INT}}(\mathsf{Q})$ and elements of this Lie group, $g_{j}(s)$,
satisfy the following differential equation

$\displaystyle\frac{d}{ds}g_{j}(s)=A_{j}\circ g_{j}(s),$

where $A_{j}\in\mathfrak{g}_{\mathrm{INT}}^{\ast}(\mathsf{Q})$ and
$g_{j}(0)=\mathrm{Id}$. The subscript $j$ denotes that action is evaluated at
$x_{j}\in\mathsf{Q}$.

Computing the tangent and cotangent lifts the actions of
$G_{\mathrm{INT}}(\mathsf{Q})$ on $\mathfrak{g}_{\mathrm{INT}}(\mathsf{Q})$
and $\mathfrak{g}^{\ast}_{\mathrm{INT}}(\mathsf{Q})$ for all orders is beyond
the scope of this project. However, we only consider the order $-1$ of
pseudodifferential symbols, because it is the symbol that would yield the
desired Lie–Poisson bracket (73).

######  Definition 4.6.

Consider the discrete symbol
$P_{j}=p_{1,j}\xi^{-1}+p_{2,j}\xi^{-2}+p_{3,j}\xi^{-3}\in\mathfrak{g}_{\mathrm{INT}}(\mathsf{Q})$
and the group action
$g_{j}(s)=1+u_{1,j}(s)\xi^{-1}+u_{2,j}(s)\xi^{-2}+u_{3,j}(s)\xi^{-3}$ and its
tangent lift of the left action on $P_{j}$ is

$\displaystyle TL_{g_{j}}P_{j}$
$\displaystyle=p_{1,j}\xi^{-1}+p_{2,j}\xi^{-2}+\left(u_{1,j}p_{1,j}+p_{3,j}\right)\xi^{-3}$
$\displaystyle+\left(u_{1,j}p_{2,j}-2u_{1,j}\mathsf{D}_{x}p_{1,j}+u_{2,j}p_{1,j}\right)\xi^{-4}+\mathcal{O}(\xi^{-5}),$

and the tangent lift of the right action on $P_{j}$ is

$\displaystyle TR_{g_{j}}P_{j}$
$\displaystyle=p_{1,j}\xi^{-1}+p_{2,j}\xi^{-2}+\left(u_{1,j}p_{1,j}+p_{3,j}\right)\xi^{-3}$
$\displaystyle+\left(u_{1,j}p_{2,j}+u_{2,j}p_{1,j}-p_{1,j}\mathsf{D}_{x}u_{1,j}\right)\xi^{-4}+\mathcal{O}(\xi^{-5}).$

Meanwhile, the cotangent lift of the action of $G_{\mathrm{INT}}(\mathsf{Q})$
on $\mathfrak{g}^{\ast}_{\mathrm{INT}}(\mathsf{Q})$ are computed according to
(53) and (53).

######  Definition 4.7.

Consider the discrete symbol $\alpha_{j}=\alpha_{1,j}\xi+\alpha_{2,j}\xi^{2}$
on the dual space of $\mathfrak{g}^{\ast}_{\mathrm{INT}}(\mathsf{Q})$, then
the cotangent lift of the left group action
$g_{j}(s)=1+u_{1,j}(s)\xi^{-1}+u_{2,j}(s)\xi^{-2}$ on $\alpha_{j}$ is

$\displaystyle TL^{\ast}_{g_{j}}\alpha_{j}$
$\displaystyle=\alpha_{2,j}\xi^{2}+\left(u_{1,j}\alpha_{2,j}+\alpha_{1,j}\right)\xi+u_{1,j}\alpha_{1,j}+u_{2,j}\alpha_{2,j}+2\alpha_{2,j}\mathsf{D}_{x}u_{1,j}$
$\displaystyle+\left(u_{2,j}\alpha_{1,j}+\alpha_{1,j}\mathsf{D}_{x}u_{1,j}+1.0\alpha_{2,j}\mathsf{D}_{x}^{2}u_{1,j}+2\alpha_{2,j}\mathsf{D}_{x}u_{2,j}\right)\xi^{-1}$
$\displaystyle+\left(\alpha_{1,j}\mathsf{D}_{x}u_{2,j}+1.0\alpha_{2,j}\mathsf{D}_{x}^{2}u_{2,j}\right)\xi^{-2}+\mathcal{O}(\xi^{-3}),$

in parallel, the cotangent lift of the right action is

$\displaystyle TR^{\ast}_{g_{j}}\alpha_{j}$
$\displaystyle=\alpha_{2,j}\xi^{2}+\left(u_{1,j}\alpha_{2,j}+\alpha_{1,j}\right)\xi+u_{1,j}\alpha_{1,j}-u_{1,j}\mathsf{D}_{x}\alpha_{2,j}+u_{2,j}\alpha_{2,j}$
$\displaystyle+\left(-u_{1,j}\mathsf{D}_{x}\alpha_{1,j}+1.0u_{1,j}\mathsf{D}_{x}^{2}\alpha_{2,j}+u_{2,j}\alpha_{1,j}-2u_{2,j}\mathsf{D}_{x}\alpha_{2,j}\right)\xi^{-1}$
$\displaystyle+\left(u_{1,j}\mathsf{D}_{x}^{2}\alpha_{1,j}-2u_{2,j}\mathsf{D}_{x}\alpha_{1,j}+3.0u_{2,j}\mathsf{D}_{x}^{2}\alpha_{2,j}\right)\xi^{-2}+\mathcal{O}(\xi^{-3}).$

In addition to the tangent and cotangent lifts, we need to define the discrete
version of the representation of the action of $G_{\mathrm{INT}}(\mathsf{Q})$
on $\mathrm{Dens}(Q)$. One issue that arises is then in (80), the action is
appears in Lagrangian coordinates, but the discrete setting is based on the
Eulerian point of view. There are two ways we can resolve this issue: Using
interpolation of the probability density function $\rho$ and numerically
integrating the continuity equation in (75). For the former, it is inspired by
smoothed particle hydrodynamics [82], where the elements of $\mathrm{Dens}(Q)$
are approximated by

$\displaystyle\rho(x,t)\approx\sum\limits_{j}\rho_{j}(t)W_{j}(x).$ (93)

Here, $W_{j}(x)$ are weights functions and one possible candidate of
$W_{j}(x)$ is the Gaussian function

$\displaystyle
W_{j}(x)=e^{-\alpha\|x-x_{j}\|^{2}},\quad\text{with}\quad\alpha\in\mathbb{R}_{>0}.$

The second approach is to use a numerical scheme for approximating
$\rho(x,t)$, where the (75) is solved in the interval $[t,t+\Delta t]$.

######  Definition 4.8.

Consider $g_{j}(s)\in G_{\mathrm{INT}}(\mathsf{Q})$ on
$\mathrm{Dens}(\mathsf{Q})$, the function (66), the map
$\zeta:G_{\mathrm{INT}}(\mathsf{Q})\to\mathrm{Dens}(\mathsf{Q})$ defined by

$\displaystyle\zeta(g_{j})\rho=\rho_{0}(g_{j}\circ
x)=\sum\limits_{j}\rho_{j}(t)W_{j}\left(x+\int_{0}^{t}u_{j}\,dt\right).$ (94)

On the other hand, it can be defined as one time iteration of the
Fokker–Planck equation

$\displaystyle\zeta(g_{j})\rho=\rho^{k}+\left(\Delta
t\right)\mathsf{D}_{x}\left(\rho_{j}^{k}u_{j}\right).$ (95)

######  Remark 4.9.

From a numerical stand point, using the pullback operator given in definition
(96) is more stable and suitable to the continuity equation, however, the one
given by definition (95) guarantees that we recover the continuity equation
when taking $\rho^{k+1}-\rho^{k}$.

The final term needed to explicitly apply the discrete Lie–Poisson
Hamilton–Jacobi theory is the representation of the Lie group action on the
dual of the Lie algebra.

######  Definition 4.10.

Consider $g_{j}(s)\in\Psi G_{d}$ on $\mathrm{Dens}^{\ast}(\mathsf{Q})$, the
discrete version of the map
$(\zeta^{\prime})^{\ast}:\mathfrak{g}_{\mathrm{INT}}(\mathsf{Q})\to\mathrm{Dens}^{\ast}(\mathsf{Q})$
defined by

$\displaystyle(\zeta^{\prime}_{f})^{\ast}\rho=\rho^{k}_{j}\mathsf{D}_{x}f,$
(96)

for all $f\in\mathrm{Dens}(\mathsf{Q})$.

######  Proposition 4.11.

The Hamiltonian system (75), written in the Lie–Poisson Hamilton–Jacobi form
(79), is solved numerically by

$\displaystyle S_{L}$ $\displaystyle=S_{0}+\Delta
tH\left(TL^{\ast}_{A}\mathbf{d}_{A}S_{0},\mathbf{d}_{\rho}S_{0}\right)$ (97)
$\displaystyle(\mu^{k},\rho^{k})$
$\displaystyle=\left(\xi+u_{j}\xi^{-1}+\mathsf{D}_{x}u_{j}\xi^{-2}+(\Delta
t)\frac{\delta H}{\delta A}\xi+(\Delta t)u_{j}\frac{\delta H}{\delta
A}\xi^{-1}\right.$ $\displaystyle\left.+(\Delta t)\frac{\delta H}{\delta
A}\mathsf{D}_{x}u_{j}\xi^{-2},\rho_{k}-\left(\Delta
t\right)\mathsf{D}_{x}\left(\rho_{k,j}u_{j}\right)+(\Delta t)\frac{\delta
H}{\delta\rho}\right)$ (98) $\displaystyle(\mu^{k+1},\rho^{k+1})$
$\displaystyle=\left(\xi+u_{j}\xi^{-1}+(\Delta t)\frac{\delta H}{\delta
A}\xi+(\Delta t)u_{j}\frac{\delta H}{\delta A}\xi^{-1}-2(\Delta
t)u_{j}\mathsf{D}_{x}\frac{\delta H}{\delta A}\xi^{-2}\right.$
$\displaystyle\left.+3(\Delta t)u_{j}\Delta^{2}_{x}\frac{\delta H}{\delta
A}\xi^{-3}+\rho_{k+1,j}\mathsf{D}_{x}f,\rho+(\Delta t)\frac{\delta
H}{\delta\rho}\right),$ (99)

with the initial condition $S_{0}=\mathrm{res}(A)+\|\rho\|^{2}$.

One thing that needs to be pointed out, that this formulation introduces a
number of open problems regarding the inputs and target data. The use of
probability density functions eliminates the role of inputs and targets and
requires us to find suitable initial and final joint distributions in order to
incorporate the training set in learning. Nevertheless, one application which
does not require picking suitable joint distributions is normalizing flow
[90]. In statistical inference, one requires to sample random data from
various distributions and sometimes such distributions are complicated, thus
renders standard techniques inefficient. In [90], it was proposed invertible
maps, implemented using neural networks, to be applied successively as a way
of mapping the target distribution to a simpler one, such as the normal
Gaussian distribution. Once that is finished, sampling from the target
distribution is done by simply sampling from the simple distribution and then
map the sample points using the learned transformation to points that belong
to the target distribution. The formulation of presented here lends itself to
such application. Indeed, in [32], a stochastic ResNet whose parameters
satisfy (77) was used to learn four different non-Gaussian distributions and
the results are shown in figure 5. Although, the numerical scheme used in [32]
uses a class of structurre-preserving integrators known as _multisymplectic
variational integrators_

Figure 5: This numerical experiment is reproduction of the result in [90]. A
stochastic ResNet is used for a normalizing flow example. The idea is to find
network parameters which would map a normal Gaussian distribution to a target
distribution, that is not necessarily Gaussian. Each row is a separate
distribution and the left most column is the target distribution. The other
columns are snapshots of the progression a normal distribution through the
layers of the stochastic ResNet. A mean–field control is used to solve the
stochastic optimal control problem used to train the network and the
parameters evolve according to (77). A multisymplectic numerical scheme for
deep learning that solves (77), presented in [32], was used to obtain this
reproduction.

##  5 Connection to Integrable Partial Differential Equations

In many cases, substituting the optimality condition $\partial_{u}H=0$ into
the Pontryagin’s Hamiltonian yields, an integrable system [50, 51]. This line
of investigation is interesting for deep learning, as it can be beneficial for
studying various behaviours. Apart from knowing the solution of the integrable
systems, they do posses interesting properties, such as the bi-Hamiltonian
structure. The existence of such structure means that the problem has another
Hamiltonian function and according to Magri’s lemma [71], we can construct an
infinite number of conservation laws by simply knowing the two Poisson
structures.

Here, we use Lie–Poisson formulation to show that deep learning training, with
certain conditions on the neural network’s parameters, is equivalent to
solving the nonlinear Schrödinger equation. One very interesting fact
regarding the Hamilton–Jacobi equation is that it works as a bridge from
classical to quantum mechanics [36]. This is done by writing the generating
function in terms of the wave function, $\psi$, i.e.

$\displaystyle\lambda(x,t)=\frac{1}{i}\mathrm{log}\psi(x,t),$ (100)

and substituting it into the Hamilton–Jacobi equation and the result is the
Schrödinger equation. For the system of equations (60), (61) and (62), instead
of (100), the following wave function is used

$\displaystyle\psi=\sqrt{\rho(x,t)}e^{i\lambda(x,t)}$ (101)

where $\rho$ satisfy (61) and $\lambda$ satisfy (62), and this transformation
is called _Madelung transformation_ and it casts the nonlinear Schrödinger
equation as the _quantum Euler equation_ [70]. This transformation has gained
attention in recent years [31, 54, 55] as it defines a symplectic
transformation from the Hilbert space of complex-valued wave functions
$H^{\infty}(\mathbb{R}^{n},\mathbb{C})$ to the dual of semidirect product Lie
algebra $\mathfrak{X}^{\ast}(\mathbb{R}^{n})\times
H^{\infty}(\mathbb{R}^{n})$.

######  Theorem 5.1.

Consider the family of Schrödiner equations

$\displaystyle
i\hbar\partial_{t}\psi+\frac{\hbar^{2}}{2m}\Delta\psi+f(|\psi|^{2})\psi,$
(102)

where $\hbar$ is the Planck’s constant, $m$ is the mass and
$f:\mathbb{R}_{>0}\to\mathbb{R}$ and whose Hamiltonian on the Hilbert space
for complex-valued wave functions
$H_{NLS}:H^{\infty}(Q,\mathbb{C})\to\mathbb{R}$ is given by

$\displaystyle
H(\psi)=\frac{\hbar}{2m}\|\nabla\psi\|^{2}+\int_{Q}F(|\psi|^{2})\,dx.$ (103)

With $f=0$, it is equivalent to the training of deep residual neural networks

$\displaystyle
J(\theta)=\int_{0}^{T}\ell_{d}\left(q(t),\theta(t)\right)\,dt+\sum_{i=0}^{N}L(\pi((q(T)^{(i)}),c^{(i)}),$
(104)

constrained to

$\displaystyle d_{t}q^{(i)}(t)=\theta(t)\circ q^{(i)}(t)\,dt+\nu dW_{t},\quad
q^{(i)}(0)=q^{(i)}_{0},\,\pi(q^{(i)}(T))=c^{(i)},$ (105)

with $i=1,\dots,N$, where $(q^{(i)}_{0},c^{(i)})$ form the dataset used in
supervised training. The main assumption imposed here is that

$\displaystyle\mathrm{div}(\theta)-\frac{\nu^{2}}{2}\mathrm{div}\left(\nabla\mathrm{log}(\rho)\right)=0$
(106)

###### Proof.

The direct equivalence between the two problems is not direct, so the first
step of the proof is to cast the training problem as a mean–field games
problem. Instead of controlling (113), the probability density function is
controlled. The dynamical constraint for the probability density function is
the Fokker–Planck equation

$\displaystyle\partial_{t}\rho+\mathrm{div}\left(\rho\sigma(x,\theta)\right)-\frac{\nu^{2}}{2}\Delta\rho=0,$
(107)

with $\rho(x,0)=\rho_{0}$ and $\rho(x,T)=\rho_{T}$. Before solving the optimal
control problem for the probability density, as mentioned in section 3.2, a
change of variable is performed. The network parameters are decomposed into

$\displaystyle\theta(x,t)=\omega(x,t)+\frac{\nu^{2}}{2}\nabla\mathrm{log}{\rho},$
(108)

and substituting it into (107), we obtain the continuity equation

$\displaystyle\partial_{t}\rho+\mathrm{div}\left(\rho\omega(x,t)\right)=0.$
(109)

The Pontryagin’s Hamiltonian then becomes

$\displaystyle
H(\rho,\lambda)=\frac{1}{2}\int_{Q}\rho\|\omega\|^{2}d\mathrm{Vol}(x)+\frac{\nu^{2}}{2}\int_{Q}\rho\left\langle\omega,\nabla\mathrm{log}(\rho)\right\rangle
d\mathrm{Vol}(x)+\frac{\nu^{4}}{8}\int_{Q}\rho\|\nabla\mathrm{log}(\rho)\|^{2}d\mathrm{Vol}(x)$

From the assumption (106) we have the condition that $\mathrm{div}\omega=0$ as
a result, the Hamiltonian further simplifies to

$\displaystyle
H(\rho,\lambda)=\frac{1}{2}\int_{Q}\rho\|\omega\|^{2}d\mathrm{Vol}(x)+\frac{\nu^{4}}{8}\int_{Q}\rho\|\nabla\mathrm{log}(\rho)\|^{2}d\mathrm{Vol}(x).$
(110)

The Hamiltonian function (110) is obtained by substituting the Madelung
transformation (101), which is a mapping from $H^{\infty}(Q,\mathbb{C})$ to
$\mathfrak{X}^{\ast}(Q)\times\mathrm{Dens}(Q)$, into the Hamiltonian,
$H_{NLS}$, for the Schrödinger equation (102), we have

$\displaystyle H(\lambda,\rho)$
$\displaystyle=\frac{\hbar}{2m}\left\langle\frac{\nabla\rho}{2\sqrt{\rho}}e^{i\frac{\lambda}{\sqrt{\hbar}}}+\frac{i}{\sqrt{\hbar}}\nabla\lambda\sqrt{\rho}e^{i\frac{\lambda}{\sqrt{\hbar}}},\frac{\nabla\rho}{2\sqrt{\rho}}e^{-i\frac{\lambda}{\sqrt{\hbar}}}-\frac{i}{\sqrt{\hbar}}\nabla\lambda\sqrt{\rho}e^{-i\frac{\lambda}{\sqrt{\hbar}}}\right\rangle$
(111)
$\displaystyle=\frac{\hbar}{8}\frac{\left\|\nabla\rho\right\|^{2}}{{\rho}}+\frac{1}{2}\rho\left\|\nabla\lambda\right\|^{2}$

if $\hbar=\nu^{4}$ and using (60), i.e. $\nabla\lambda=\omega$, we obtain the
Hamiltonian (110). ∎

######  Corollary 5.2.

Consider the deep neural network training problem

$\displaystyle
J(\theta)=\int_{0}^{T}\ell_{d}\left(q(t),\theta(t)\right)\,dt+\sum_{i=0}^{N}L(\pi((q(T)^{(i)}),c^{(i)}),$
(112)

constrained to the ensemble of stochastic differential equations

$\displaystyle d_{t}q^{(i)}(t)=\theta(t)\circ q^{(i)}(t)\,dt+\nu dW_{t},\quad
q^{(i)}(0)=q^{(i)}_{0},\,\pi(q^{(i)}(T))=c^{(i)},$ (113)

and the pair of initial and final conditions are the training dataset
$(q^{(i)}_{0},c^{(i)})$, where $i=1,\dots,N$. The stochastic control problem
is equivalent to the Lie–Poisson system with Hamiltonian
$H:\Psi\mathrm{DS}(\mathbb{C})^{\ast}\to\mathbb{R}$ given by

$\displaystyle H(L)=\frac{i}{3}\mathrm{Tr}(L^{3}),$ (114)

where $L\in\Psi\mathrm{DS}(\mathbb{C})^{\ast}$ is defined as

$\displaystyle L=\partial_{x}+\psi\partial_{x}^{-1}\psi^{\ast}.$ (115)

With the Hamiltonian defined on $\Psi\mathrm{DS}(\mathbb{C})^{\ast}$, the dual
algebra of complex-valued pseudodifferential symbols, we can simply apply
theorem 3.1.

The he nonlinear Schrödinger equation is not only integrable equation that
arises from deep learning, in fact, when equation (77), when $\nu\to 0$, and
the network’s parameters regularising norm $\ell$ is chosen as $H^{1}$ Sobolev
norm, the equation becomes

$\displaystyle\partial_{t}m$ $\displaystyle=\pm\mathrm{ad}^{\ast}_{u}m,$
$\displaystyle m=u-\alpha^{2}\partial_{x}^{2}u$

or when expressed in coordinates on $\mathbb{R}^{N}$, it becomes

$\displaystyle\partial_{t}m_{i}$
$\displaystyle+\partial_{x^{j}}\left(m_{i}u^{j}\right)+\left(m_{j}\right)\partial_{x^{i}}u^{j}=0$
(116) $\displaystyle m_{i}$
$\displaystyle=u_{i}-\alpha^{2}\partial^{2}_{x^{i}}u_{i},$

and it is the _dispersionless Camassa–Holm equation_ [18] and it is an example
of a Lie–Poisson system. An interesting property of this partial differential
equation is that its weak solution is also a Hamiltonian system

$\displaystyle H(q,p,t)$
$\displaystyle=\sum\limits_{i=0,k=0}^{N_{b},N_{b}}p^{(i)}p^{(k)}e^{-\frac{1}{2}\|q^{(i)}-q^{(k)}\|^{2}}$
(117)

The physics that are described by this Hamiltonian is an approximation to the
N-particle solution, or _peakon train_ , of Euler–Poincaré equation with
$H^{1}$ Sobolev norm [60, 81]. What is meant by the peakon or peaked soliton
is the a type of soliton, but with a discontinuity or a cusp at the crest of
the wave. For the actual peakon train, the Gaussian in the Hamiltonian (117),
is replaced by Green’s function $G(x,y)=e^{-|x-y|}x$ for shifted Poisson
differential operator $I_{d}-\partial^{2}_{x}$.

Further, the Hamilton’s equation for the peakon train do have a Lax pair, and
it is an example of an integrable system.

$\displaystyle\frac{d}{dt}L=[L,P],$ (118)

where $L$ an $N\times N$ matrix, known as Lax matrix, with entries
$L_{i,j}=\sqrt{p^{(i)}p^{(k)}}e^{-|q^{(i)}-q^{(k)}|}$ and $P$ is also an
$N\times N$ matrix and its entries are
$P_{i,j}=-2\sqrt{p^{(i)}p^{(k)}}\mathrm{sign}(q^{(i)}-q^{(k)})e^{-|q^{(i)}-q^{(k)}|}$.
One advantage of knowing the Lax matrix is that taking trace of its power,
i.e. $\mathrm{Tr}(L^{k})$ for $k=2,\dots,N$, results in the constant of motion
[18].

##  6 Conclusion & Future Direction

We have introduced, in this paper, an equivalent formulation of deep learning
as a hydrodynamics system. This made possible by substituting the deep neural
network by an ordinary differential equation and formulating the training of
deep learning as an optimal control problem as it is equivalent to solving a
Hamiltonian system. First, it was the use of mean–field type control for the
stochastic optimal control problem and that approach allows for the
application of Pontryagin’s Maximum Principle. The Hamilton’s equations that
are the optimality condition for deep learning, in this case, are partial
differential equations, and instead of using them directly, we have used an
equivalent system of equations that is similar to quantum fluid equation. Here
we have demonstrated that this optimality condition is a Lie–Poisson system.
Interestingly, the Poisson structure for deep learning coincides with the
Poisson structure for the ideal compressible fluid.

What is interesting about the Lie–Poisson formulation is that it links deep
learning with equations from fluid dynamics, in some specific cases, to the
nonlinear Schrödinger equation. Such connection have not been well studied,
and it is a future research direction as the nonlinear Schrödinger equation,
in one-dimension at least, has a hierarchy of integrable equations.
Potentially, this could help deriving conserved quantities for deep learning.
The importance of conserved quantities is they allow us to better analyse the
performance of the system and construct Lyapunov functions for studying the
stability of the system; this approach forms the basis of Arnold’s method for
stability.

Despite the fact Lie–Poisson systems are Hamiltonian and with a Poisson
structure, it is challenging to solve numerically and for this
structure–preserving integrators are required. One type of
structure–preserving integrators we have explored in this paper is integrators
based on the Hamilton–Jacobi theory and its symmetry reduced version the
Lie–Poisson Hamilton–Jacobi theory. The main reason for this choice is it can
be used for both Hamiltonian systems on symplectic cotangent bundles and
Poisson manifolds. Moreover, one of the advantages of this approach is that
can be easily be made high-order without the need to impose additional
constraints to enforce the conservation of the symplectic structure. This
forms the basis of a training algorithm, and we have shown the algorithm
behaves such similar methods using Euler method and $4^{th}$ order Runge–Kutta
method.

The application is not straight forward and as a result there has been a
number of modifications are needed. On top of that, the underlying dynamics
are defined on infinite dimensional semidirect product space. For that reason,
we have extended the Lie–Poisson Hamilton–Jacobi theory to accommodate
semidirect product space and we have also formulated the hydrodynamics problem
on the group of pseudodifferential symbols. The last step is required as it
allows us to apply Lie–Poisson Hamilton–Jacobi theory directly to the quantum
fluid equation for deep learning and derive structure–preserving numerical
schemes.

Stability of Deep Learning: As mentioned, the Arnold’s method of stability can
be applied to the nonlinear Schrödinger equation as a surrogate way to study
the stability of deep learning, however, other methods that reply on the
geometry of the problem can be used. One of such methods is the energy–Casimir
method of [45, 5], where it is used to study the nonlinear stability of
Lie–Poisson systems’ relative equilibria and it is applied to various ideal
compressible fluid and even plasma physics. The notion of nonlinear stability
for a trajectory $(M(t),\rho(t))\in\mathfrak{s}^{\ast}$ is that there exists
$\epsilon>0$ such that $\|(M_{0},\rho_{0})\|<\epsilon$, then
$\|((M(t)-M_{0},\rho(t)-\rho_{0})\|<\delta$ for all of $t$, where $\delta>0$.
In other words, the trajectory, whose initial condition belongs to a specific
set, would remain close to the set of initial values. The main ingredient is
the Casimir function, which is a map $C:\mathfrak{s}^{\ast}\to\mathbb{R}$ such
that it is the centre of the Lie algebra meaning for every function
$F:\mathfrak{s}^{\ast}\to\mathbb{R}$, the Lie–Poisson bracket is equal to
zero, i.e. $\\{C,F\\}=0$. The Hamiltonian is then augmented with the Casimir
function to form a new Hamiltonian $H_{C}=H+\Phi(C)$ and first variation of
the Hamiltonian is evaluated at the equilibrium and then equated to zero. From
there we find the function $\Phi$. The next step is to determine the convexity
of the second variation of the Hamiltonian $H_{C}$. Once all criteria are met,
one can determine the nonlinear stability of the system.

The major challenge with the energy–Casimir method is that it requires a
knowledge of the Casimir method, and this motivated the development of
energy–momentum method of [95]. Instead, the augmented Hamiltonian is
$H_{\xi}=H+\left\langle\mathbf{J},\xi\right\rangle$, where
$\mathbf{J}:T^{\ast}Q\to\mathfrak{s}^{\ast}$ and $\xi\in\mathfrak{s}$. The Lie
algebra element $\xi$ plays the role of the vector field associated with the
relative equilibria and it renders the first variation $\delta H_{\xi}=0$.
Same as energy–Casimir method, the definiteness of $\delta^{2}H_{\xi}$
indicates the stability of the system.

Another route for studying stability of Hamiltonian systems is to find the
action-angle variables, denoted by
$(I_{1},\dots,I_{n},\theta_{1},\dots,\theta_{n})$ In this coordinate system,
the Hamiltonian $H$ depends only on the actions $I_{1},\dots,I_{n}$ and the
result of that is the configuration manifold is an n-torus $\mathbb{T}^{n}$.
These variables are needed in order to study the perturbation of the
Hamiltonian and we can assess the stability by KAM theory [59, 83, 6]. That
said, KAM theory works on integrable systems and, at this point it is yet to
be investigates, it is not known whether or not the deep learning system is
integrable or not. What we know is that for a specific case, we can the
Hamiltonian obtained from Pontryagin’s Maximim Principle is integrable, which
leads us to the next possible research direction.

Generating Function For Optimal Control: The issue with using
structure–preserving integrators and gradient methods is that they require an
initial guess for the costate and control parameters. Such guesses can affect
the behaviour of the Hamiltonian system and convergence in general and in
addition, iterative schemes for such problems can be time consuming and
computationally taxing. A different approach that is related to dynamic
programming was proposed in [37] and further investigated in [86, 38, 92, 39].
The key ingredient is to compute the generating function for symplectomorphism
that maps between the initial and final conditions. This further makes it
possible to find the solution of the costate in terms of training dataset,
thus sidestepping the need to guess the initial/final condition for the
costate. Such research direction is worth investigating for training of deep
neural networks.

The idea is to treat initial and terminal conditions, $q(0)=q_{0}\in Q$ and
$q(T)=q_{T}\in Q$ respectively, as old and new coordinates, $(q_{0},p_{0})$
and $(q_{T},p_{T})$. To see how generating functions and optimal transport are
related, we start by looking at dynamic programming approach to optimal
control. Given the cost functional

$\displaystyle\mathbb{S}(q)=\int_{0}^{T}\ell(\theta)\,dt+\Phi(q(T),q_{T}),$
(119)

we define the _value function_ as

$\displaystyle
V(q(t),t)=\min\limits_{\theta}\left\\{V(q(t+dt),t+dt)+\int_{t}^{t+dt}\ell(\theta)\,dt\right\\},$
(120)

and it satisfies the _Hamilton–Jacobi–Bellman equation_ , which can be
regarded as a generalisation of the Hamilton–Jacobi equation. As shown in
(Park), the value function is

$\displaystyle V(q(t),t)=-S(q(t),q_{T},t)+\Phi(q(T),q_{T})$ (121)

which mean that it can be obtained by solving the Hamilton–Jacobi equation for
$S$. That said, one small technical issue: it is not practical to use $S$
directly and that is because at $t=0$ it is hard to find initial conditions
that satisfy the identity transformation. For this reason, we use type-II
generating function $S_{2}$, or type-III generating function $S_{3}$ instead.
Once $V$ is known, the parameters that would train the network are then given
by

$\displaystyle\theta=\arg\max\limits_{\theta^{\ast}}\left(H\left(q,\frac{\partial
V}{\partial q},\theta^{\ast},t\right)\right)$ (122)

That said, finding a solution for the generating function is challenging on
its own right, however, in [37] it was shown that it can be solved using power
series expansion of generating function whose coefficients are obtained by
solving a system of ordinary differential equations. The derivation of the
system of ordinary differential equations is done by applying the Taylor
expansion to the Hamiltonian and the generating function. Once that is done,
both expansions are substituted into the Hamilton–Jacobi equation and as a
result we obtain a power series whose coefficients are differential equations.
Setting the coefficients of the combined series to be equal to zero, we obtain
the system of equations whose solution is the coefficients of the generating
function $S$.

The $N$ terms Taylor expansion of the Hamiltonian, around the null point of
$T^{\ast}Q$, is given by

$\displaystyle\widetilde{H}(q,p,t)=\sum_{i_{q_{1}}=0}^{N}\cdots\sum_{i_{q_{1}}=0}^{N}\frac{(q_{1})^{i_{q_{1}}}\cdots(q_{n})^{i_{q_{n}}}(p_{1})^{i_{p_{1}}}\cdots(p_{n})^{i_{p_{n}}}}{{i_{q_{1}}}!\cdots
i_{p_{1}}!}\,\left(\frac{\partial^{i_{q_{1}}+\cdots+i_{p_{n}}}}{\partial
q_{1}^{i_{q_{1}}}\cdots\partial p_{n}^{i_{p_{n}}}}H\right)(0,\ldots,0,t),$

likewise, the generating function, is written as a polynomial with $N$ terms,

$\displaystyle
F(q,P,t)=\sum_{i_{q_{1}}=0}^{N}\cdots\sum_{i_{q_{1}}=0}^{N}\frac{(q_{1})^{i_{q_{1}}}\cdots(q_{n})^{i_{q_{n}}}(p_{1})^{i_{p_{1}}}\cdots(p_{n})^{i_{p_{n}}}}{{i_{q_{1}}}!\cdots
i_{p_{1}}!}\,F_{i_{q_{1}},\ldots,i_{p_{n}}}(t),$

where $F_{i_{q_{1}},\ldots,i_{p_{n}}}(t)\in\mathcal{C}(Q)$ are the
coefficients of polynomial $F$ and they are differentiable functions on $Q$.
Substituting the expansions $F(q,P,t)$ and $\widetilde{H}(q,p,t)$ into the
Hamilton–Jacobi equation, we obtain a system of differential equations for the
coefficients of $F(q,P,t)$. Once the solution is obtained, that yields the
generating function $F(q,P,t)$, which is then used to obtain (121), which in
turn allows us to obtain the weights of the network (122).

This formulation is in the continuous-time space and for large systems, it is
practical to use discrete time formulation. Mirroring the approach in the
previous section, a discrete version of the generating function approach for
solving two-boundary point problems could be beneficial. The Hamilton–Jacobi
equation becomes a discrete one and it was introduced in [85]. Given the set
of discrete points of the generating function $\\{S^{k}\\}_{k=1}^{N}$ which
satisfy the discrete Hamilton–Jacobi equation

$\displaystyle S^{k+1}(q_{k+1})=S^{k}(q_{k+1})+D_{2}S^{k+1}(q_{k+1})\cdot
q_{k+1}+H^{+}_{d}(q_{k},D_{2}S^{k+1}(q_{k+1})),$ (123)

with the symplectic mapping $(q_{k},p_{k})\to(q_{k+1},p_{k+1})$ defined by

$\displaystyle q_{k+1}=q_{k}+D_{1}S^{k+1}(q_{k+1}),$ (124) $\displaystyle
p_{k+1}=p_{k}-D_{2}S^{k+1}(q_{k+1}).$ (125)

Once $\\{S^{k}\\}_{k=1}^{N}$ are obtained, we can then use (121) to obtain $V$
and from it $\theta$ of (122).

## References

  * [1] L. Abia and J. M. Sanz-Serna.  Partitioned Runge-Kutta Methods for Separable Hamiltonian Problems.  Mathematics of Computation, 60(202):617–634, 1993.  Publisher: American Mathematical Society.
  * [2] R. Abraham, J. E. Marsden, and T. Ratiu.  Manifolds, Tensor Analysis, and Applications.  Applied Mathematical Sciences. Springer-Verlag, New York, 2 edition, 1988.
  * [3] A. Agrachev, M. Zelikin, Y. Sachkov, Y. Sachkov, and Û. Sačkov.  Control Theory from the Geometric Viewpoint.  Control theory and optimization. Springer, 2004.
  * [4] A. Arnaudon, N. Ganaba, and D. D. Holm.  The stochastic energy-Casimir method.  Comptes Rendus Mécanique, 346(4):279–290, Apr. 2018.
  * [5] A. Arnaudon, N. Ganaba, and D. D. Holm.  The stochastic energy-casimir method.  Comptes Rendus Mécanique, 346(4):279–290, 2018.
  * [6] V. I. Arnol’d.  Proof of a theorem of A. N. Kolmogorov on the invariance of quasi-periodic motions under small perturbations of the Hamiltonian.  Russian Mathematical Surveys, 18:9–36, 1963.
  * [7] V. I. Arnold.  Sur la géométrie différentielle des groupes de lie de dimension infinie et ses applicationsa l’hydrodynamique des fluides parfaits.  Ann. Inst. Fourier, 16(1):319–361, 1966.
  * [8] Y. Barash, J. A. Calarco, W. Gao, Q. Pan, X. Wang, O. Shai, B. J. Blencowe, and B. J. Frey.  Deciphering the splicing code.  Nature, 465(7294):53–59, May 2010.  Number: 7294 Publisher: Nature Publishing Group.
  * [9] J.-D. Benamou and Y. Brenier.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  Numerische Mathematik, 84(3):375–393, 2000.
  * [10] M. Benning, E. Celledoni, M. J. Ehrhardt, B. Owren, and C.-B. Schönlieb.  Deep learning as optimal control problems: Models and numerical methods.  Journal of Computational Dynamics, 6(2):171–198, 2019.
  * [11] A. Bensoussan, J. Frehse, and P. Yam.  The Mean Field Type Control Problems.  In A. Bensoussan, J. Frehse, and P. Yam, editors, Mean Field Games and Mean Field Type Control Theory, SpringerBriefs in Mathematics, pages 15–29. Springer, New York, NY, 2013.
  * [12] M. Blaszak.  Multi-Hamiltonian Theory of Dynamical Systems.  Theoretical and Mathematical Physics. Springer Berlin Heidelberg, 2012\.
  * [13] A. M. Bloch and P. E. Crouch.  Optimal Control, Optimization, and Analytical Mechanics.  In J. Baillieul and J. C. Willems, editors, Mathematical Control Theory, pages 268–321. Springer, New York, NY, 1999.
  * [14] A. M. Bloch, D. D. Holm, P. E. Crouch, and J. E. Marsden.  An optimal control formulation for inviscid incompressible ideal fluid flow.  In Proceedings of the 39th IEEE Conference on Decision and Control (Cat. No. 00CH37187), volume 2, pages 1273–1278. IEEE, 2000.
  * [15] A. M. Bloch, I. I. Hussein, M. Leok, and A. K. Sanyal.  Geometric structure-preserving optimal control of a rigid body.  Journal of Dynamical and Control Systems, 15(3):307–330, July 2009\.
  * [16] M. Born and E. Wolf.  Principles of optics: Electromagnetic theory of propagation, interference and diffraction of light.  Elsevier Science, 2013.
  * [17] F. Bullo and A. D. Lewis.  Supplementary chapters for _Geometric Control of Mechanical Systems_ [FB/ADL:04], jan 2005.
  * [18] R. Camassa and D. D. Holm.  An integrable shallow water equation with peaked solitons.  Physical Review Letters, 71(11):1661–1664, Sept. 1993.  Publisher: American Physical Society.
  * [19] C. Carathéodory.  Die Grundlagen der geometrischen Optik.  In C. Carathéodory, editor, Geometrische Optik, Ergebnisse der Mathematik und Ihrer Grenzgebiete, pages 15–35. Springer, Berlin, Heidelberg, 1937.
  * [20] E. Celledoni, M. J. Ehrhardt, C. Etmann, R. I. McLachlan, B. Owren, C.-B. Schönlieb, and F. Sherry.  Structure preserving deep learning.  arXiv:2006.03364 [cs, math, stat], June 2020.  arXiv: 2006.03364.
  * [21] P. J. Channell and C. Scovel.  Symplectic integration of Hamiltonian systems.  Nonlinearity, 3(2):231–259, May 1990.
  * [22] H. Chen, O. Engkvist, Y. Wang, M. Olivecrona, and T. Blaschke.  The rise of deep learning in drug discovery.  Drug Discovery Today, 23(6):1241–1250, June 2018.
  * [23] T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud.  Neural ordinary differential equations.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 6571–6583. Curran Associates, Inc., 2018.
  * [24] C. Cotter, D. Holm, and P. Hydon.  Multisymplectic formulation of fluid dynamics using the inverse map.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 463(2086):2671–2687, 2007.
  * [25] E. Dupont, A. Doucet, and Y. W. Teh.  Augmented neural odes.  arXiv preprint arXiv:1904.01681, 2019.
  * [26] W. E.  A Proposal on Machine Learning via Dynamical Systems.  Communications in Mathematics and Statistics, 5(1):1–11, mar 2017\.
  * [27] W. E, J. Han, and Q. Li.  A mean-field optimal control formulation of deep learning.  Research in the Mathematical Sciences, 6(1):10, dec 2018.
  * [28] D. G. Ebin and J. Marsden.  Groups of diffeomorphisms and the motion of an incompressible fluid.  Annals of Mathematics, 92(1):102–163, 1970.
  * [29] K. Engø and S. Faltinsen.  Numerical Integration of Lie–Poisson Systems While Preserving Coadjoint Orbits and Energy.  SIAM Journal on Numerical Analysis, 39(1):128–145, Jan. 2001.  Publisher: Society for Industrial and Applied Mathematics.
  * [30] C. Farabet, C. Couprie, L. Najman, and Y. LeCun.  Learning Hierarchical Features for Scene Labeling.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(8):1915–1929, Aug. 2013.  Conference Name: IEEE Transactions on Pattern Analysis and Machine Intelligence.
  * [31] D. Fusca.  The Madelung transform as a momentum map.  Journal of Geometric Mechanics, 9(2):157, 2017.
  * [32] N. Ganaba.  Multisymplectic formulation of deep learning.  In preparation, 2021.
  * [33] N. Ganaba.  Stochastic Lie-Poisson-Hamilton-Jacobi Theory.  In preparation, 2021.
  * [34] E. Gawehn, J. A. Hiss, and G. Schneider.  Deep Learning in Drug Discovery.  Molecular Informatics, 35(1):3–14, 2016.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/minf.201501008.
  * [35] F. Gay-Balmaz and T. S. Ratiu.  Clebsch optimal control formulation in mechanics.  Journal of Geometric Mechanics, 3(1):41, 2011.
  * [36] H. Goldstein and C. Poole.  Classical Mechanics.  Addison-Wesley series in physics. Addison-Wesley Publishing Company, 1980\.
  * [37] V. Guibout and D. Scheeres.  Formation Flight with Generating Functions: Solving the Relative Boundary Value Problem.  In AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Monterey, California, Aug. 2002. American Institute of Aeronautics and Astronautics.
  * [38] V. M. Guibout.  The hamilton-jacobi theory for solving two-point-boundary value-problems: Theory and numerics with application to spacecraft formation flight.  In Optimal Control, and the Study of Phase Space Structure. University of Michigan at Ann Arbor, Ph. D. Dissertation, 2004.
  * [39] V. M. Guibout and D. J. Scheeres.  Solving Relative Two-Point Boundary Value Problems: Spacecraft Formulation Flight Transfers Application.  Journal of Guidance, Control, and Dynamics, 27(4):693–704, 2004\.
  * [40] E. Haber and L. Ruthotto.  Stable architectures for deep neural networks.  Inverse Problems, 34(1):014004, dec 2017.
  * [41] S.-P. Han.  Superlinearly convergent variable metric algorithms for general nonlinear programming problems.  Mathematical Programming, 11(1):263–282, Dec. 1976.
  * [42] G. Hinton, L. Deng, D. Yu, G. E. Dahl, A. R. Mohamed, N. Jaitly, A. Senior, V. Vanhoucke, P. Nguyen, T. N. Sainath, and B. Kingsbury.  Deep Neural Networks for Acoustic Modeling in Speech Recognition: The Shared Views of Four Research Groups.  IEEE Signal Processing Magazine, 29(6):82–97, Nov. 2012.  Conference Name: IEEE Signal Processing Magazine.
  * [43] D. D. Holm.  Euler’s fluid equations: Optimal control vs optimization.  Physics Letters A, 373(47):4354–4359, Nov. 2009.
  * [44] D. D. Holm and B. A. Kupershmidt.  Poisson brackets and clebsch representations for magnetohydrodynamics, multifluid plasmas, and elasticity.  Physica D: Nonlinear Phenomena, 6(3):347–363, Apr. 1983.
  * [45] D. D. Holm, J. E. Marsden, T. Ratiu, and A. Weinstein.  Nonlinear stability of fluid and plasma equilibria.  Physics Reports, 123(1):1 – 116, 1985.
  * [46] D. D. Holm, J. E. Marsden, and T. S. Ratiu.  The Euler-Poincaré Equations and Semidirect Products with Applications to Continuum Theories.  Advances in Mathematics, 137(1):1–81, July 1998.
  * [47] P. Horava.  On a covariant Hamilton-Jacobi framework for the Einstein-Maxwell theory.  Classical and Quantum Gravity, 8(11):2069–2084, Nov. 1991.
  * [48] O. Junge, J. E. Marsden, and S. Ober-Blöbaum.  Discrete mechanics and optimal control.  IFAC Proceedings Volumes, 38(1):538–543, 2005.
  * [49] V. Jurdjevic.  Geometric Control Theory.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 1996.
  * [50] V. Jurdjevic.  Integrable Hamiltonian systems on symmetric spaces: Jacobi, Kepler and Moser.  arXiv e-prints, 1103:arXiv:1103.2818, Mar. 2011.
  * [51] V. Jurdjevic.  Optimal Control and Geometry: Integrable Systems.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2016.
  * [52] F. Kang, W. Hua-mo, Q. Meng-zhao, and W. Dao-liu.  Construction of canonical difference schemes for hamiltonian formalism via generating functions.  Journal of Computational Mathematics, 7(1):71–96, 1989.
  * [53] J. B. Keller.  Corrected bohr-sommerfeld quantum conditions for nonseparable systems.  Annals of Physics, 4(2):180–188, June 1958.
  * [54] B. Khesin, G. Misiolek, and K. Modin.  Geometric hydrodynamics via Madelung transform.  Proceedings of the National Academy of Sciences, 115(24):6165–6170, jun 2018.
  * [55] B. Khesin, G. Misiolek, and K. Modin.  Geometric Hydrodynamics of Compressible Fluids.  arXiv:2001.01143 [math-ph], jan 2020.  arXiv: 2001.01143.
  * [56] B. Khesin and R. Wendt.  The geometry of infinite-dimensional groups, volume 51.  Springer Science & Business Media, 2008.
  * [57] B. Khesin and I. Zakharevich.  Poisson-Lie group of pseudodifferential symbols.  Communications in Mathematical Physics, 171(3):475–530, Aug. 1995\.
  * [58] M. B. Kobilarov and J. E. Marsden.  Discrete Geometric Optimal Control on Lie Groups.  IEEE Transactions on Robotics, 27(4):641–655, aug 2011.
  * [59] A. N. Kolmogorov.  On conservation of conditionally periodic motions for a small change in hamilton’s function.  In Dokl. Akad. Nauk SSSR, volume 98, pages 527–530, 1954.
  * [60] S. Kouranbaeva.  The Camassa-Holm equation as a geodesic flow on the diffeomorphism group.  Journal of Mathematical Physics, 40(2):857–868, Jan. 1999.
  * [61] D. Kraft.  On Converting Optimal Control Problems into Nonlinear Programming Problems.  In K. Schittkowski, editor, Computational Mathematical Programming, NATO ASI Series, pages 261–280, Berlin, Heidelberg, 1985\. Springer.
  * [62] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  ImageNet classification with deep convolutional neural networks.  Communications of the ACM, 60(6):84–90, May 2017.
  * [63] T. Lee, M. Leok, and N. H. McClamroch.  Computational Geometric Optimal Control of Rigid Bodies.  Communications in Information & Systems, 8(4):445–472, Jan. 2008\.
  * [64] T. Lee, M. Leok, and N. H. McClamroch.  Optimal Attitude Control of a Rigid Body Using Geometrically Exact Computations on SO(3).  Journal of Dynamical and Control Systems, 14(4):465–487, Oct. 2008\.
  * [65] M. Leok and J. Zhang.  Discrete Hamiltonian variational integrators.  IMA Journal of Numerical Analysis, 31(4):1497–1532, Oct. 2011.
  * [66] M. K. Leung, H. Y. Xiong, L. J. Lee, and B. J. Frey.  Deep learning of the tissue-regulated splicing code.  Bioinformatics, 30(12):i121–i129, 2014.
  * [67] Q. Li, L. Chen, C. Tai, and E. Weinan.  Maximum principle based algorithms for deep learning.  The Journal of Machine Learning Research, 18(1):5998–6026, 2017\.
  * [68] Q. Li and S. Hao.  An optimal control approach to deep learning and applications to discrete-weight neural networks.  In J. Dy and A. Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 2985–2994. PMLR, 10–15 Jul 2018.
  * [69] J. Ma, R. P. Sheridan, A. Liaw, G. E. Dahl, and V. Svetnik.  Deep Neural Nets as a Method for Quantitative Structure–Activity Relationships.  Journal of Chemical Information and Modeling, 55(2):263–274, Feb. 2015.  Publisher: American Chemical Society.
  * [70] E. Madelung.  Quantentheorie in hydrodynamischer Form.  Zeitschrift für Physik, 40(3):322–326, Mar. 1927.
  * [71] F. Magri.  A simple model of the integrable hamiltonian equation.  Journal of Mathematical Physics, 19(5):1156–1162, 1978.
  * [72] J. Marsden, G. Misiolek, J. Ortega, M. Perlmutter, and T. Ratiu.  Hamiltonian Reduction by Stages.  Lecture Notes in Mathematics. Springer Berlin Heidelberg, 2007.
  * [73] J. Marsden and A. Weinstein.  Coadjoint orbits, vortices, and Clebsch variables for incompressible fluids.  Physica D: Nonlinear Phenomena, 7(1):305–323, may 1983.
  * [74] J. E. Marsden, S. Pekarsky, and S. Shkoller.  Discrete euler-poincaré and lie-poisson equations.  Nonlinearity, 12(6):1647–1662, oct 1999.
  * [75] J. E. Marsden, T. Ratiu, R. Schmid, R. Spencer, and A. J. Weinstein.  Hamiltonian systems with symmetry, coadjoint orbits and plasma physics.  Atti della Accademia delle scienze di Torino, 117(1):289–340, 1983\.
  * [76] J. E. Marsden, T. Ratiu, and A. Weinstein.  Semidirect Products and Reduction in Mechanics.  Transactions of the American Mathematical Society, 281(1):147–177, 1984.
  * [77] J. E. Marsden and T. S. Ratiu.  Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems, volume 17.  Springer Science & Business Media, 2013.
  * [78] J. E. Marsden and M. West.  Discrete mechanics and variational integrators.  Acta Numerica, 10:357–514, 2001.
  * [79] R. I. McLachlan.  Explicit Lie-Poisson integration and the Euler equations.  Physical Review Letters, 71(19):3043–3046, Nov. 1993.  Publisher: American Physical Society.
  * [80] M. McLean and L. K. Norris.  Covariant field theory on frame bundles of fibered manifolds.  Journal of Mathematical Physics, 41(10):6808–6823, Sept. 2000.
  * [81] G. Misiołek.  A shallow water equation as a geodesic flow on the Bott-Virasoro group.  Journal of Geometry and Physics, 24(3):203–208, Feb. 1998.
  * [82] J. J. Monaghan.  Smoothed particle hydrodynamics.  Annual review of astronomy and astrophysics, 30(1):543–574, 1992\.
  * [83] J. Moser.  On invariant curves of area-preserving mappings of an annulus.  Nachrichten der Akademie der Wissenschaften in Göttingen. II. Mathematisch-Physikalische Klasse, 1962:1–20, 1962.  Publisher: Vandenhoeck & Ruprecht, Göttingen.
  * [84] P. Mullen, K. Crane, D. Pavlov, Y. Tong, and M. Desbrun.  Energy-preserving integrators for fluid animation.  In ACM Transactions on Graphics (TOG), volume 28, page 38. ACM, 2009\.
  * [85] T. Ohsawa, A. Bloch, and M. Leok.  Discrete hamilton-jacobi theory.  SIAM J. Control. Optim., 49:1829–1856, 2011.
  * [86] C. Park and D. J. Scheeres.  Solutions of the optimal feedback control problem using hamiltonian dynamics and generating functions.  In 42nd IEEE International Conference on Decision and Control (IEEE Cat. No. 03CH37475), volume 2, pages 1222–1227. IEEE, 2003.
  * [87] D. Pavlov, P. Mullen, Y. Tong, E. Kanso, J. E. Marsden, and M. Desbrun.  Structure-preserving discretization of incompressible fluids.  Physica D: Nonlinear Phenomena, 240(6):443–458, 2011.
  * [88] L. D. Persio and M. Garbelli.  Deep Learning and Mean-Field Games: A Stochastic Optimal Control Perspective.  Symmetry, 13(1):14, Jan. 2021.
  * [89] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko.  The mathematical theory of optimal processes.  1962\.
  * [90] D. Rezende and S. Mohamed.  Variational inference with normalizing flows.  In International Conference on Machine Learning, pages 1530–1538. PMLR, 2015.
  * [91] L. Ruthotto and E. Haber.  Deep Neural Networks Motivated by Partial Differential Equations.  Journal of Mathematical Imaging and Vision, sep 2019.
  * [92] D. J. Scheeres, C. Park, and V. M. Guibout.  Solving optimal control problems with generating functions.  Advances in the Astronautical Sciences, 116:1185–1205, 2003.
  * [93] F. Seide, G. Li, and D. Yu.  Conversational Speech Transcription Using Context-Dependent Deep Neural Networks.  page 4.
  * [94] J. Simo, N. Tarnow, and K. Wong.  Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics.  Computer Methods in Applied Mechanics and Engineering, 100(1):63–116, 1992.
  * [95] J. C. Simo, D. Lewis, and J. E. Marsden.  Stability of relative equilibria. part i: The reduced energy-momentum method.  Archive for Rational Mechanics and Analysis, 115:15–59, 1991.
  * [96] Y. Sun, L. Zhang, and H. Schaeffer.  NeuPDE: Neural Network Based Ordinary and Partial Differential Equations for Modeling Time-Dependent Data.  arXiv:1908.03190 [cs, stat], Aug. 2019.  arXiv: 1908.03190.
  * [97] H. J. Sussmann.  Geometry and Optimal Control.  In J. Baillieul and J. C. Willems, editors, Mathematical Control Theory, pages 140–198. Springer, New York, NY, 1999.
  * [98] I. Sutskever, O. Vinyals, and Q. V. Le.  Sequence to sequence learning with neural networks.  In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2, NIPS’14, page 3104–3112, Cambridge, MA, USA, 2014. MIT Press.
  * [99] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich.  Going Deeper with Convolutions.  arXiv:1409.4842 [cs], Sept. 2014.  arXiv: 1409.4842.
  * [100] J. Vankerschaver, C. Liao, and M. Leok.  Generating functionals and Lagrangian partial differential equations.  Journal of Mathematical Physics, 54(8):082901, Aug. 2013.
  * [101] J. von Rieth.  The Hamilton-Jacobi theory of De Donder and Weyl applied to some relativistic field theories.  JMP, 25(4):1102–1115, Apr. 1984.
  * [102] E. Weinan.  A proposal on machine learning via dynamical systems.  Communications in Mathematics and Statistics, 5(1):1–11, 2017.
  * [103] E. Weinan, J. Han, and Q. Li.  A mean-field optimal control formulation of deep learning.  Research in the Mathematical Sciences, 6(1):10, 2019.
  * [104] H. Y. Xiong, B. Alipanahi, L. J. Lee, H. Bretschneider, D. Merico, R. K. Yuen, Y. Hua, S. Gueroussov, H. S. Najafabadi, T. R. Hughes, et al.  The human splicing code reveals new insights into the genetic determinants of disease.  Science, 347(6218), 2015.
  * [105] H. Yoshida.  Construction of higher order symplectic integrators.  Physics Letters A, 150(5):262–268, Nov. 1990.
  * [106] V. Zeitlin.  Finite-mode analogs of 2D ideal hydrodynamics: Coadjoint orbits and local canonical structure.  Physica D: Nonlinear Phenomena, 49(3):353–362, Apr. 1991.
  * [107] G. Zhong and J. E. Marsden.  Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators.  Physics Letters A, 133(3):134–139, 1988.