Lie-Poisson Neural Networks (LPNets): Data-Based Computing of
  Hamiltonian Systems with Symmetries
An accurate data-based prediction of the long-term evolution of Hamiltonian
systems requires a network that preserves the appropriate structure under each
time step. Every Hamiltonian system contains two essential ingredients: the
Poisson bracket and the Hamiltonian. Hamiltonian systems with symmetries, whose
paradigm examples are the Lie-Poisson systems, have been shown to describe a
broad category of physical phenomena, from satellite motion to underwater
vehicles, fluids, geophysical applications, complex fluids, and plasma physics.
The Poisson bracket in these systems comes from the symmetries, while the
Hamiltonian comes from the underlying physics. We view the symmetry of the
system as primary, hence the Lie-Poisson bracket is known exactly, whereas the
Hamiltonian is regarded as coming from physics and is considered not known, or
known approximately. Using this approach, we develop a network based on
transformations that exactly preserve the Poisson bracket and the special
functions of the Lie-Poisson systems (Casimirs) to machine precision. We
present two flavors of such systems: one, where the parameters of
transformations are computed from data using a dense neural network (LPNets),
and another, where the composition of transformations is used as building
blocks (G-LPNets). We also show how to adapt these methods to a larger class of
Poisson brackets. We apply the resulting methods to several examples, such as
rigid body (satellite) motion, underwater vehicles, a particle in a magnetic
field, and others. The methods developed in this paper are important for the
construction of accurate data-based methods for simulating the long-term
dynamics of physical systems.
Christopher Eldred, François Gay-Balmaz, Sofiia Huraka, Vakhtang Putkaradze
# Lie-Poisson Neural Networks (LPNets): Data-Based Computing of Hamiltonian Systems with Symmetries

###### Abstract

An accurate data-based prediction of the long-term evolution of Hamiltonian systems requires a network that preserves the appropriate structure under each time step. Every Hamiltonian system contains two essential ingredients: the Poisson bracket and the Hamiltonian. Hamiltonian systems with symmetries, whose paradigm examples are the Lie-Poisson systems, have been shown to describe a broad category of physical phenomena, from satellite motion to underwater vehicles, fluids, geophysical applications, complex fluids, and plasma physics. The Poisson bracket in these systems comes from the symmetries, while the Hamiltonian comes from the underlying physics. We view the symmetry of the system as primary, hence the Lie-Poisson bracket is known exactly, whereas the Hamiltonian is regarded as coming from physics and is considered not known, or known approximately. Using this approach, we develop a network based on transformations that exactly preserve the Poisson bracket and the special functions of the Lie-Poisson systems (Casimirs) to machine precision. We present two flavors of such systems: one, where the parameters of transformations are computed from data using a dense neural network (LPNets), and another, where the composition of transformations is used as building blocks (G-LPNets). We also show how to adapt these methods to a larger class of Poisson brackets. We apply the resulting
methods to several examples, such as rigid body (satellite) motion, underwater vehicles, a particle in a magnetic field, and others. The methods developed in this paper are important for the construction of accurate data-based methods for simulating the long-term dynamics of physical systems.

keywords: Neural equations, Data-based modeling, Long-term evolution, Hamiltonian Systems, Poisson brackets +
Footnote †: journal: Journal of Computational Physics

## 1 Introduction

### Relevant work

Machine Learning (ML) approaches to data assimilation and modeling have been successful in interpreting large amounts of unstructured data. However, for many scientific and industrial purposes, direct applications of ML methods for big data have been developed without understanding that the underlying engineering and physics are challenging. To address these problems, _Physics Informed Neural Networks (PINNs)_ have been developed (Raissi et al., 2019). As it is impossible to provide a complete overview of the PINNs literature here, we refer the reader to recent reviews (Cuomo et al., 2022; Karniadakis et al., 2021) for references and more thorough review of the literature. In that approach, one assumes that the motion of the system \(\mathbf{u}(t)\) is approximated by a law of motion

\[\dot{\mathbf{u}}=\mathbf{f}(\mathbf{u},t), \tag{1}\]

where \(\mathbf{u}\) and \(\mathbf{f}\) can be either finite-dimensional, forming a system of ODEs, or infinite dimentional, forming a system of PDEs. In the 'vanilla' PINNs approach, the initial conditions, boundary conditions, and the ODEs itself form a part of the loss function optimized during the learning procedure. The last part - computing the difference \(\dot{\mathbf{u}}-\mathbf{f}(\mathbf{u},t)\) - can be achieved for a neural network approximating \(\mathbf{u}\), that is, approximating the mapping \((\mathbf{u}_{0},t)\rightarrow\mathbf{u}\). Similarly, for PDEs in a \((1+1)\) dimension, for example, the solution \(u(x,t)\) is a mapping from the position \(x\) and time \(t\) to the space of solutions. In PINNs, one approximates this mapping by a neural network. For a given set of weights and activating functions in neural network, at the given set of points in space and time \((t_{i},\mathbf{x}_{i})\), one can compute the _exact_ values of temporal and spatial derivatives using the method of automatic differentiation (Baydin et al., 2018). One then constructs the Mean-Square-Error (MSE) of the PDE approximation of solution and boundary conditions, taken as the cost function, and optimizes network weights to minimize that cost. The advantage of PINNs is their computational efficiency: speedups more than 100,00x
for evaluations of solutions of complex systems like weather have been reported (Bi et al., 2023; Pathak et al., 2022), although the quantification of speedup still needs to be understood (Karniadakis et al., 2021). PINNs are thus extremely useful in practical applications, and have been implemented as the open source using _Python_ and _Julia_, as well as Nvidia's _Modulus_.

In spite of the successes, there are still many uncertainties in the applications of PINNs. For example, it is known that PINNs may struggle with systems having dynamics with widely different time scales (Karniadakis et al., 2021), or having trouble in finding the true minimum for optimization (Krishnapriyan et al., 2021). Finally, systems with very little friction, like Hamiltonian systems, are difficult to describe by the standard PINNs, since evaluation errors created on every time step tend to accumulate with no chance of dissipating quickly enough. This paper develops a method of simulating a type of non-canonical Hamiltonian system, namely, Lie-Poisson systems. These systems are extremely important in building models of physical systems in Eulerian or body coordinates. The method is also shown to generalize to some other types of non-canonical Poisson brackets. The remainder of the literature review is dedicated to work on the use of data-based methods in simulating Hamiltonian systems. In order to make the discussion more concrete, we introduce a brief background in Hamiltonian and Poisson systems. More details on this can be found in Section 2.

There are several approaches applying physics-informed data methods to computations of Hamiltonian systems. Most of the work has been focused on _canonical_ Hamiltonian systems, _i.e._, the systems where the law of motion (1) has the following particular structure: \(\mathbf{u}\) is \(2n\)-dimensional, \(\mathbf{u}=(\mathbf{q},\mathbf{p})\), and there is a function \(H(\mathbf{q},\mathbf{p})\), the Hamiltonian, such that \(\mathbf{f}\) in (1) becomes

\[\mathbf{f}=\mathbb{J}\nabla_{\mathbf{u}}H\,,\quad\mathbb{J}=\left(\begin{array} []{cc}0&\mathbb{I}_{n}\\ -\mathbb{I}_{n}&0\end{array}\right)\,, \tag{2}\]

with \(\mathbb{I}_{n}\) the \(n\times n\) identity matrix, leading to the _canonical_ Hamilton equations for \((\mathbf{q},\mathbf{p})\):

\[\dot{\mathbf{q}}=\frac{\partial H}{\partial\mathbf{p}}\,,\quad\dot{\mathbf{p} }=-\frac{\partial H}{\partial\mathbf{q}}\,. \tag{3}\]

The evolution of an arbitrary phase space function \(F(\mathbf{q},\mathbf{p})\) along a solution of (3) is then described by the canonical Poisson bracket :

\[\frac{dF}{dt}=\{F,H\}=\frac{\partial F}{\partial\mathbf{q}}\frac{\partial H} {\partial\mathbf{p}}-\frac{\partial H}{\partial\mathbf{q}}\frac{\partial F}{ \partial\mathbf{p}}. \tag{4}\]
The bracket (4) is a mapping sending two smooth functions of \((\mathbf{q},\mathbf{p})\) into a smooth function of the same variables. This mapping is bilinear, antisymmetric, acts as a derivation on both functions, and satisfies the Jacobi identity: for all functions \(F,G,H\)

\[\{\{F,G\},H\}+\{\{H,F\},G\}+\{\{G,H\},F\}=0\,. \tag{5}\]

Brackets that satisfy all the required properties, _i.e._, are bilinear, antisymmetric, act as a derivation and satisfy Jacobi identity (5), but are not described by the canonical equations (4), are called general Poisson brackets (also known as non-canonical Poisson brackets). The corresponding equations of motion are called non-canonical Hamiltonian (or Poisson) systems. Often, these brackets have a non-trivial null space leading to the conservation of certain quantities known as the Casimir constants, or simply _Casimirs_. The Casimirs are properties of the Poisson bracket and are independent of a particular realization of a given Hamiltonian. This paper will focus on the data-based approaches for computations of an important class of non-canonical Poisson systems.

There is an avenue of thought that focuses on learning the actual Hamiltonian for the system from data, or, in the more general case, the Poisson bracket. This approach was explicitly implemented for canonical Hamiltonian systems in (Greydanus et al., 2019) under the name of _Hamiltonian Neural Networks (HNN)_, which approximated the Hamiltonian function \(H(\mathbf{q},\mathbf{p})\) fitting evolution of the particular data sequence through equations (3). It was shown that embedding the dynamics with the knowledge of the data allows a much more accurate and robust way to approximate the solution compared to a general Neural Network (NN). This work was further extended to include the adaptive learning of parameters and transitions to chaos (Han et al., 2021). The mathematical background guaranteeing the existence of Hamiltonian function sought in HNN was derived in (David and Mehats, 2021). An alternative method of learning the equations is given by the _Lagrangian Neural Networks (LNNs)_(Cranmer et al., 2020) which approximates the solutions of Euler-Lagrange equations, _i.e._, the equations in the coordinate-velocity space \((\mathbf{q},\dot{\mathbf{q}})\)_before_ Legendre-transforming to the momentum-coordinate representation \((\mathbf{q},\mathbf{p})\) given by equations (3). More generally, learning a vector field for the non-canonical Poisson brackets was suggested in (Sipka et al., 2023). The main challenge in that work was enforcing Jacobi identity (5) for the learned equation structure.

In these and other works on the topic, one learns the vector field governing the system, with the assumption that the vector field can be solved using appropriate numerical methods. However, one needs to be aware that care must be
taken in computing the numerical solutions for Hamiltonian systems, especially for long-term computations, as regular numerical methods lead to distortion of quantities that should be conserved, such as total energy and, when appropriate, the momenta. In order to compute long-term evolution of systems obeying the Hamiltonian vector fields, whether exact or approximated by the Hamiltonians derived from neural networks, one can use variational integrator methods (Hall and Leok, 2015; Leok and Shingel, 2012; Marsden and West, 2001) that conserve momenta-like quantities with machine precision. However, these integrators may be substantially more computationally intensive compared to non-structure preserving methods.

In this paper, we focus on an alternative approach, namely, exploring the learning transformations in phase space that satisfy appropriate properties. Recall that if \(\mathbf{\phi}(\mathbf{u})\) is a map in the phase space of equation (3) with \(\mathbf{u}=(\mathbf{q},\mathbf{p})\), then this map is called symplectic if

\[\left(\frac{\partial\mathbf{\phi}}{\partial\mathbf{u}}\right)^{T}\mathbb{J}\left( \frac{\partial\mathbf{\phi}}{\partial\mathbf{u}}\right)=\mathbb{J}\,. \tag{6}\]

A well-known result of Poincare states that the flow \(\mathbf{\phi}_{t}(\mathbf{u})\) of the canonical Hamiltonian system (3), sending initial conditions \(\mathbf{u}\) to the solution at time \(t\), is a symplectic map (Arnol'd, 2013; Marsden and Ratiu, 2013). Several authors pursued the idea of searching directly for the symplectic mappings obtained from the data, instead of finding actual equations of the canonical systems and then solving them.

Perhaps the first work in the area of computing the symplectic transformations directly was done in (Chen et al., 2020), where _Symplectic Recurring Neural Networks (SRNNs)_ were designed. The SRNNs computed an approximation to the symplectic transformation from data using the appropriate formulas for symplectic numerical methods. An alternative method of computation of canonical Hamilton equations for non-separable Hamiltonians was done in (Xiong et al., 2020), building approximations to symplectic steps using a symplectic integrator suggested in (Tao, 2016). This technique was named _Non-Separable Symplectic Neural Networks (NSSNNs)_.

A more direct computation of symplectic mappings was done using three different methods in (Chen and Tao, 2021; Jin et al., 2020). The first approach derived in (Jin et al., 2020) computes the dynamics through the composition of symplectic maps of certain type, which was implemented under the name of _SympNets_. Another approach (Chen and Tao, 2021) derives the mapping directly using a generating function approach for canonical transformations, implemented as _Generating Function Neural Networks (GFNNs)_. The approach using GFNNs al
lows for an explicit estimate of the error in a long-term simulation. In contrast, error analysis in SympNets focuses on the local approximation error. Finally, Burby et al. (2020) developed _HenonNets_, Neural Networks based on the Henon mappings, capable of accurately learning Poincare maps of a Hamiltonian systems while preserving the symplectic structure. SympNets, GFNNs and HenonNets showed the ability to accurately simulate long-term behavior of simple integrable systems like a pendulum or a single planet orbit; and satisfactory long-term simulation for chaotic systems like the three-body plane problem. Thus, learning symplectic transformations directly from data shows great promise for long-term simulations of Hamiltonian systems.

The method of SympNets was further extended for non-canonical Poisson systems by transforming the non-canonical form to local canonical coordinates using the Lie-Darboux theorem and subsequently using SympNets (Jin et al., 2022) by assuming that the dynamics occurs within a neighborhood in which the Poisson structure has constant rank. This method was named _Poisson Neural Networks (PNNs)_. While this method can in principle treat any Poisson system by learning the transformation to the canonical variables and its inverse, we shall note that there are several difficulties associated with this approach:

1. It is known that the Lie-Darboux transformation, in general, is only local: a global function transforming the system to canonical coordinates may not exist, although such transformation did exist for all examples presented in (Jin et al., 2022). When such non-locality happens, on would need to define several transformations in overlapping domains and ensure the smoothness between them. It is not clear how the network-based Lie-Darboux function will perform in that case.
2. If an _absolutely accurate_ Lie-Darboux transformation coupled with a symplectic integrator was used to transform the coordinates to canonical form in (Jin et al., 2022), it would of course preserve all the Casimirs. However, any numerical errors in determining that transformation will yield corresponding errors in the Casimir evolution.
3. Since the errors in Casimir evolution are determined by the errors in the Lie-Darboux mapping, it is also not clear how these errors will accumulate over long-term evolution of the system.

The preservation of Casimirs is especially important for predicting the probabilistic properties for the long-term evolution of many trajectories in the Poisson system (Dubinkina and Frank, 2007). In particular, even when the errors in each
individual component of the solution may accumulate over time, the fact that the solution stays exactly on the Casimir surface will play an essential role in the probability distribution in phase space.

The SympNets and PNN approach was further extended in (Bajars, 2023) where volume-preserving neural networks _LocSympNets_ and their symmetric extensions _SymLocSympNets_ were derived, based on the composition of mappings of certain type. A consistent good accuracy of long-term solutions obtained by LocSympNets and SymLocSympNets was demonstrated for several problems, including a discretized linear advection equation, rigid body dynamics and a particle in magnetic field. Although the methods of (Bajars, 2023) did not explicitly appeal to the Poisson structure of equations, the efficiency of the methods was demonstrated as applied to several problems that are essentially Poisson in nature, such as rigid body motion and the motion of particle in a magnetic field. However, the extension of the theory to more general problems was hampered by the fact that the completeness of activation matrices suggested in (Bajars, 2023) was not yet known.

The limitations of the methods of (Jin et al., 2022) and (Bajars, 2023) come from the fact that they relate to a general system of equations, where it is assumed that very little is known about the general system apart the fact that it is Hamiltonian, or Poisson. On the other hand, there is a large class of physical systems where the Poisson bracket is known exactly; in particular, _Lie-Poisson_ systems. In these approaches, the bracket does not come as a consequence of the equations of motion, but from general considerations of the Lie group symmetries of the system. In that case, the Poisson bracket has an explicit expression stemming from the expression for the Lie algebra bracket, and is called the Lie-Poisson bracket. The choice of actual Hamiltonian is a secondary step coming from physics. For example, the motion of the rigid body comes from the invariance of the Lagrangian/Hamiltonian in the body frame with respect to rigid rotations, _i.e._ SO(3) symmetry, see the introduction to the theory in Section 2 below. Table 1 provides an incomplete list of physical examples admitting a description through Lie-Poisson or closely related brackets.

In this paper, we develop the methods of constructing the activation maps directly from the brackets, predicting dynamics built out of maps computed by particular explicit solutions of Lie-Poisson systems. The method is applicable to _all_ finite-dimensional Lie-Poisson systems, as well as any Poisson system where explicit integration of appropriate equations for appropriate transformations are available. The advantage of utilizing explicit integration in Lie-Poisson equations to drastically speed up calculations on every time step was already noticed
in (McLachlan, 1993), although the application was limited to Hamiltonians of certain form depending on the Poisson bracket. Our method is also applicable to arbitrary Hamiltonians; in fact, the Hamiltonian does not need to be known for the construction of the Neural Network. However, it is assumed that the underlying symmetry and the appropriate Lie-Poisson bracket are known. This is indeed often the case, as the detailed expression of the Hamiltonian in terms of the variables

\begin{table}
\begin{tabular}{|c|c|} \hline Problem & Reference \\ \hline Rigid body & (Holm et al., 2009) \\  & (Marsden and Ratiu, 2013) \\ \hline Heavy top & (Holm et al., 2009) \\  & (Marsden and Ratiu, 2013) \\ \hline Underwater vehicles & (Leonard, 1997) \\  & (Leonard and Marsden, 1997) \\  & (Holmes et al., 1998) \\ \hline Plasmas & (Morrison, 1980), \\  & (Marsden and Weinstein, 1982), \\  & (Holm et al., 1985), \\  & (Holm and Tronci, 2010) \\ \hline Fluids & (Marsden and Weinstein, 1983), \\  & (Marsden et al., 1984), \\  & (Holm et al., 1985), \\  & (Morrison, 1998), \\  & (Morrison et al., 2006), \\ \hline Geophysical fluid dynamics & (Weinstein, 1983), \\  & (Holm, 1986), \\  & (Salmon, 2004) \\ \hline Complex and nematic fluids & (Holm, 2002), \\  & (Gay-Balmaz and Ratiu, 2009), \\  & (Gay-Balmaz and Tronci, 2010) \\ \hline Molecular strand dynamics & (Ellis et al., 2010), \\  & (Gay-Balmaz et al., 2012) \\ \hline Fluid-structure interactions & (Gay-Balmaz and Putkaradze, 2019) \\ \hline Hybrid quantum-classical dynamics & (Gay-Balmaz and Tronci, 2022), \\  & (Gay-Balmaz and Tronci, 2023) \\ \hline \end{tabular}
\end{table}
Table 1: A short description of some physical problems that can be written in Lie-Poisson form and related Poisson brackets.
is often only approximate and is driven by the modeling choices for the particular physical system (Tonti, 2013), in contrast to the Lie-Poisson bracket itself.

The novel contributions of the paper are as follows:

1. We show that a large class of Poisson systems, namely Lie-Poisson systems obtained by symmetry reduction, allow explicit construction of maps for every particular Lie-Poisson bracket.
2. These maps are global, in contrast to Lie-Darboux coordinates that are only local and may need to be re-computed depending on the position of solution in the configuration manifold.
3. By construction, these maps preserve Casimirs of each particular bracket exactly, which is not necessarily the case for PNNs, LocSympNets, SymLocSympNets and any other methods known to us.

## 2 An introduction to the Lie-Poisson equations

### General introduction: from Lagrangian to Hamiltonian description

We start with a brief introduction of the origin of Lie-Poisson systems to introduce some notations and show why this particular type of approach is essential for many physical problems. We will consider only finite-dimensional systems here not to make the consideration too abstract. Suppose a mechanical system is described by the coordinates \(\mathbf{q}\) and velocities \(\dot{\mathbf{q}}\), with \(\mathbf{q}\) lying on some configuration manifold \(Q\) of dimension \(n\). Hamilton's action principle is based on the Lagrangian function \(L(\mathbf{q},\dot{\mathbf{q}})\) (possibly depending on time) and on the action \(S=\int_{t_{0}}^{t_{f}}L(\mathbf{q},\dot{\mathbf{q}})\mathrm{d}t\), and imposes the condition that the variations of the action must vanish on the variations of \(\mathbf{q}\) that are fixed on the boundaries \(t=t_{0},t_{f}\)

\[\delta S=\delta\int_{t_{0}}^{t_{f}}L(\mathbf{q},\dot{\mathbf{q}})\mathrm{d}t= 0\,,\quad\delta\mathbf{q}(t_{0})=\delta\mathbf{q}(t_{f})=0\,. \tag{7}\]

In the Lagrangian approach, one takes the variations of (7) and gets the Euler-Lagrange equations, which are second order equations in \(\mathbf{q}\). In the Hamiltonian approach, one introduces the momenta \(\mathbf{p}=\frac{\partial L}{\partial\dot{\mathbf{q}}}\) and assumes that this relation can be inverted for each \(\mathbf{q}\), giving the velocities as \(\dot{\mathbf{q}}=\dot{\mathbf{q}}(\mathbf{q},\mathbf{p})\). One then defines the Hamiltonian function \(H(\mathbf{p},\mathbf{q})=\mathbf{p}\cdot\dot{\mathbf{q}}(\mathbf{q},\mathbf{p} )-L(\mathbf{q},\dot{\mathbf{q}}(\mathbf{q},\mathbf{p}))\), and the Euler-Lagrange equations of motion are equivalent to the _canonical_ Hamilton equations:

\[\dot{\mathbf{q}}=\frac{\partial H}{\partial\mathbf{p}}\,,\quad\dot{\mathbf{p}} =-\frac{\partial H}{\partial\mathbf{q}}\,. \tag{8}\]
Any function \(F(\mathbf{q},\mathbf{p})\) evolves according to the _canonical bracket_

\[\frac{dF}{dt}=\{F,H\}=\frac{\partial F}{\partial\mathbf{q}}\cdot\frac{\partial H }{\partial\mathbf{p}}-\frac{\partial F}{\partial\mathbf{p}}\cdot\frac{\partial H }{\partial\mathbf{q}}\,. \tag{9}\]

One can see that the bracket (9) satisfies the following properties that are valid for any functions \(F(\mathbf{q},\mathbf{p})\), \(G(\mathbf{q},\mathbf{p})\) and \(H(\mathbf{q},\mathbf{p})\):

1. Antisymmetry: \(\{F,H\}=-\{H,F\}\),
2. Linearity in each component: \(\{aF+bG,H\}=a\{F,H\}+b\{G,H\}\), \(a,b\in\mathbb{R}\),
3. Leibniz rule (acts as a derivation): \(\{FG,H\}=F\{G,H\}+G\{F,H\}\),
4. Jacobi identity: \(\{F,\{G,H\}\}+\{H,\{F,G\}\}+\{G,\{H,F\}\}=0\).

After defining the total phase space \(\mathbf{u}=(\mathbf{q},\mathbf{p})\), the canonical bracket (4) and Hamilton's equations of motion can be written in coordinates as

\[\{F,H\}=\frac{\partial F}{\partial\mathbf{u}}\cdot\mathbb{J}^{-1}\frac{ \partial H}{\partial\mathbf{u}}\,,\quad\dot{\mathbf{u}}=\mathbb{J}^{-1}\frac{ \partial H}{\partial\mathbf{u}}\,,\quad\mathbb{J}=\left(\begin{array}{cc}0& \mathbb{I}_{n}\\ -\mathbb{I}_{n}&0\end{array}\right)\,, \tag{10}\]

where \(\mathbb{I}_{n}\) is the \(n\times n\) unit matrix. From the definition of \(\mathbb{J}\) above, \(\mathbb{J}^{-1}=\mathbb{J}^{T}=-\mathbb{J}\).

### General Poisson brackets

A general Poisson bracket \(\{F,H\}\) satisfies all four properties of the canonical Poisson bracket above, but cannot be necessarily expressed in coordinates as (10). Instead, the bracket and corresponding equations of motion are described more generally in local coordinates as

\[\{F,H\}=\frac{\partial F}{\partial\mathbf{u}}\cdot\mathbb{B}(\mathbf{u})\frac {\partial H}{\partial\mathbf{u}}\,,\qquad\dot{\mathbf{u}}=\mathbb{B}(\mathbf{ u})\frac{\partial H}{\partial\mathbf{u}}\,. \tag{11}\]

In order for the bracket (11) to be Poisson, the matrix \(\mathbb{B}(\mathbf{u})\), also known as the Poisson tensor, must be antisymmetric and satisfy a condition involving both the matrix \(\mathbb{B}(\mathbf{u})\) and its derivatives. If such bracket can be found, the system is called (non-canonical) Poisson.

A special attention should be paid to the case when the matrix \(\mathbb{B}\) is degenerate. In that case, there very often1 are special functions for that particular
bracket, called _Casimirs_, which are conserved for _any_ Hamiltonian. By definition, a Casimir function \(C\) satisfies

\[\{F,C\}=0\quad\text{for all $F$}\,. \tag{12}\]

Any evolution must occur in such a way that all Casimirs of the system are conserved. Geometrically, the motion is only possible on the intersection of Casimir level sets, no matter what the Hamiltonian for the system is.

As we have seen, the motion of mechanical systems without friction is governed by canonical Poisson brackets, therefore the appearance of a non-canonical bracket (11) is somewhat mysterious. It turns out that in presence of symmetries, one can drastically reduce the degrees of freedom of the canonical Hamiltonian system by expressing the dynamics in terms of reduced coordinates, such as body or spatial coordinates. In such reduced variables, the dynamics is still governed by a Poisson bracket, but which is no more canonical. One of the most important example is the case of a system whose configuration is a Lie group: \(Q=G\). Let us consider the case of a rigid body. The configuration manifold is the Lie group of rotation matrices in \(\mathbb{R}^{3}\), also known as \(G=SO(3)\). This group consists of all \(3\times 3\) orthogonal matrices \(\Lambda\) (_i.e._, \(\Lambda^{T}\Lambda=\mathbb{I}_{3}\) with the determinant equal to 1). The Lagrangian depends on the variables \((\Lambda,\dot{\Lambda})\), and is just the kinetic energy. If one were to write either the Euler-Lagrange equations or the canonical Hamilton equations for \(\Lambda\) obtained by parameterizing \(SO(3)\) using local representation of rotation matrices, such as rotation angles, one would end up with complicated and unwieldy equations. Instead, the equations of motion for the rigid body can be efficiently written in terms of the angular velocity in the body frame \(\mathbf{\Omega}\) using the tensor of inertia \(\mathbb{I}\) as:

\[\mathbb{I}\dot{\mathbf{\Omega}}=-\mathbf{\Omega}\times\mathbb{I}\mathbf{ \Omega}\,. \tag{13}\]

As it turns out, the equations (13) can be understood from the point of view of symmetry. The rigid body kinetic energy is invariant with respect to _left_ rotations \(\Lambda\to R\Lambda\), where \(R\in SO(3)\) is a fixed rotation matrix. We can thus express the kinetic energy in terms of the antisymmetric \(3\times 3\) matrices \(\widehat{\Omega}=\Lambda^{T}\dot{\Lambda}\) which take on the role of angular velocities. One can compute the vector \(\mathbf{\Omega}\) as \(\widehat{\Omega}_{ij}=\epsilon_{ijk}\mathbf{\Omega}_{k}\), with \(\epsilon_{ijk}\) being the Levi-Civita symbol. The equations of motion (13) can be written in terms of momenta \(\mathbf{\Pi}\), with the utilization of the rigid body bracket:

\[\dot{\mathbf{\Pi}}=-\frac{\partial H}{\partial\mathbf{\Pi}}\times\mathbf{\Pi }\,,\quad\{F,H\}:=-\mathbf{\Pi}\cdot\left(\frac{\partial F}{\partial\mathbf{ \Pi}}\times\frac{\partial H}{\partial\mathbf{\Pi}}\right)\,. \tag{14}\]

The bracket in (14), as it turns out, satisfies the four properties of a Poisson bracket and is expressed in terms of the Lie bracket on the Lie algebra of \(SO(3)\), given
by the vector product, stemming from the invariance of the system with respect to \(SO(3)\) rotations.

These ideas can be generalized as follows, see (Marsden and Ratiu, 2013, Chap.10). For a general system defined on a Lie group \(G\), one naturally has the associated Lie algebra \(\mathfrak{g}\) with Lie bracket denoted as \([\alpha,\beta]\) for \(\alpha,\beta\in\mathfrak{g}\). If \(\{e_{a}\}\), \(a=1,\ldots,n\) is a basis of \(\mathfrak{g}\), then the Lie bracket is locally expressed in terms of the structure constants \(C^{d}_{ab}\) such that

\[[e_{a},e_{b}]=C^{d}_{ab}e_{d}. \tag{15}\]

Let us denote by \(\langle\mu,\alpha\rangle\) the duality pairing between vectors \(\alpha\) in the Lie algebra \(\mathfrak{g}\) and co-vectors (or momenta) \(\mu\) in the dual space \(\mathfrak{g}^{*}\) to \(\mathfrak{g}\). The partial derivatives of functions \(F,H:\mathfrak{g}^{*}\rightarrow\mathbb{R}\) with respect to \(\mu\) thus belong to \(\mathfrak{g}\), and one can define the _Lie-Poisson bracket_ derived from the Lie bracket as follows:

\[\{F,H\}=\pm\left\langle\mu,\left[\frac{\partial F}{\partial\mu}\,,\frac{ \partial H}{\partial\mu}\right]\right\rangle\,. \tag{16}\]

We refer to A for the explanation of the \(\pm\) sign and for the relation between this bracket and the canonical bracket. In terms of coordinates, the bracket (16) is a particular case of (11) with the matrix \(\mathbb{B}(\mu)\) defined as

\[\mathbb{B}_{ab}(\mu)=\pm C^{d}_{ab}\mu_{d} \tag{17}\]

and the Lie-Poisson equations are expressed in coordinates as

\[\dot{\mu}_{a}=\pm C^{d}_{ab}\mu_{d}\frac{\partial H}{\partial\mu_{b}}\,. \tag{18}\]

One can verify that the Lie-Poisson bracket (16) satisfies all the conditions of a general Poisson bracket. While the Lie-Poisson brackets appear because of the fundamental considerations of the symmetries of the physical system, the Hamiltonian is a modelling choice related to physics. Thus, we assume that if a system possesses a Lie-Poisson bracket, it is known explicitly _a priori_ and does not need to be computed or determined from the data. However, the dependence of the Hamiltonian on its arguments and parameters is not known and must be determined.

## 3 LPNets as Poisson maps for a general system

We are now ready to consider the theory of structure-preserving neural networks. We follow the ideas originally introduced in (Jin et al., 2020) as SympNets. The idea behind SympNets is to learn in the space of available symplectic
transformations for a canonical system. Further, (Jin et al., 2022) introduced a generalization of this method to an arbitrary Poisson system by using the Lie-Darboux theorem, stating that every Poisson system can be locally transformed into a canonical form, implementing these ideas a PoissonNets. This theory was further extended in (Patel et al., 2022; Zhang et al., 2022) for thermodynamics systems using the GENERIC metriplectic bracket approach.

Mathematically speaking, previous works in this avenue of thinking (Jin et al., 2022, 2020) and (Bajars, 2023) seek to determine a mapping \(\phi_{h}(\mathbf{y}_{i})\) in the next computational point \(\mathbf{y}_{i+1}\), satisfying as many restrictions preserving the structure of the actual flow as possible. The authors developed universal mappings that are symplectic in the canonical coordinates \((\mathbf{q},\mathbf{p})\), and sought the solution as a combination of these mappings. Since the mappings are universal for all finite-dimensional Hamiltonian systems, one then needed to prove the convergence and completeness of the combinations of these mapping, especially when coupled to the Lie-Darboux theorem which maps the Poisson system to a canonical form.

We build on the ideas of SympNets and Poisson nets by directly constructing the transformations related to the Poisson systems (_Poisson maps_) for the _known_ Lie-Poisson bracket, and using them as elements for construction of LPNets.

We first recall the useful concept of _Poisson maps_.

**Definition 3.1**.: _(Marsden and Ratiu, 2013, SS10.3) Let \((P_{i},\{\cdot,\cdot\}_{i})\), \(i=1,2\) be two Poisson manifolds. A mapping \(f:P_{1}\to P_{2}\) is Poisson if it preserves the Poisson bracket, i.e., for all functions \(F,G:P_{2}\to\mathbb{R}\)_

\[\{F,G\}_{2}\circ f=\{F\circ f,G\circ f\}_{1}\,. \tag{19}\]

A critical piece of information for our further progress is contained in the following result.

**Theorem 3.2** (Hamiltonian flows are Poisson).: _(Marsden and Ratiu, 2013, Thm. 10.3.1.) Consider a Poisson manifold \((P,\{\cdot,\cdot\})\) and a Hamiltonian \(H:P\to\mathbb{R}\). Let \(\phi_{t}(\mathbf{u}_{0})\) be the flow of the Poisson system associated with \(H\), see (11), which maps the initial conditions \(\mathbf{u}_{0}\) to the solution \(\mathbf{u}(t)\) at time \(t\). Then, the mapping \(\phi_{t}\) is Poisson, i.e., it satisfies:_

\[\{F\circ\phi_{t},G\circ\phi_{t}\}=\{F,G\}\circ\phi_{t}\,, \tag{20}\]

_for all functions \(F,G:P\to\mathbb{R}\)._
Following this theorem, we can design Poisson maps for a particular Poisson bracket using a sequence of flows created by simplified Hamiltonians for the Lie-Poisson dynamics. The advantage of this approach is that the resulting transformations will be Poisson and will be able to approximate any flow locally for the particular system considered. The disadvantage of our approach is inherently intertwined with the advantages: the mappings have to be constructed explicitly for every Lie-Poisson bracket. Fortunately, as we show below, this is possible as these mappings are derived as solutions of a linear system of ODEs. The method presented here is reminiscent of the Hamiltonian splitting methods used in numerical analysis (McLachlan, 1993; McLachlan and Quispel, 2002), reformulated for the purpose of data-based computations and the use of neural networks to find appropriate parameters of the Hamiltonians.

## 4 A general application to finite-dimensional Lie groups

The key to this paper lies in considering the evolution of the momentum \(\mu\) coming from the equations (18) for particular expressions for the Hamiltonian, namely, Hamiltonians linear in momenta

\[H(\mu)=\langle\alpha,\mu\rangle=\alpha^{a}\mu_{a}\,, \tag{21}\]

where \(\alpha^{a}\) are some constants that are to be found based on the learning procedure. The flow generated by the Hamiltonian (21) is given by a linear equation in \(\mu\) that can be written in two equivalent ways (choosing the \(+\) sign in (18)):

\[\dot{\mu}_{a}=C^{d}_{ab}\alpha^{b}\mu_{d}:=\mathbb{M}(\alpha)^{d}_{a}\mu_{d}= \mathbb{N}(\mu)_{ab}\alpha^{b}\,. \tag{22}\]

The number of possible dimensions of \(\alpha\) is exactly equal to the dimension of the momentum space. However, the "effective" dimension of this space may be less, and is related to the dimension of the Lie algebra \(n\) minus the dimension of the kernel of the operator \(\mathbb{N}(\mu):\mathfrak{g}\to\mathfrak{g}^{*}\), see Remark 4.1. The operators \(\mathbb{M}(\alpha)\) acting on the space of momenta \(\mu\) and \(\mathbb{N}(\mu)\) acting on the space of \(\alpha\) can be described in the coordinate-free form as \(\mathrm{ad}^{*}_{\alpha}\,\square\) and \(\mathrm{ad}^{*}_{\square}\,\mu\), respectively, see A.

Let us assume that there are \(d\) independent Casimir functions \(C_{j}(\mu)\), \(j=1,...,d\). From (12) and (17) such functions satisfy \(C^{d}_{ab}\frac{\partial C_{j}}{\partial\mu_{b}}\mu_{d}=0\) for all \(\mu\), i.e., \(\partial C_{j}/\partial\mu\) must belong to the kernel of the operator \(\mathbb{N}(\mu)\). We shall consider the effective dimension of the space of all possible \(\alpha\) to be exactly \(n-d\), where \(n\) is the dimension of the momentum space and \(d\) is the number of independent Casimirs. Thus, locally the vectors \(\alpha\) in this effective space form exactly the right number
of local tangent vectors to the intersection of Casimir surfaces to reach any point locally using a flow generated by the Poisson map. The next step is to compute that map exactly.

**Remark 4.1** (On the effective number of dimensions).: _In general, the number of null directions of \(\mathbb{N}(\mu_{0})\), denoted as \(k(\mu_{0})\) may depend on the momenta \(\mu_{0}\in\mathfrak{g}^{*}\) (or, more precisely, on the actual coadjoint orbit the solution is on), whereas the number of independent Casimirs \(d\) is fixed. The dimension of the image of the map \(\xi\in\mathfrak{g}\mapsto\mathbb{N}(\mu_{0})\xi=\mathrm{ad}_{\xi}^{*}\,\mu_{0} \in\mathfrak{g}^{*}\) is the dimension of the coadjoint orbit of \(\mu_{0}\), denoted as \(\mathcal{O}_{\mu_{0}}\). Thus, in general, we have \(k(\mu_{0})+\dim(\mathcal{O}_{\mu_{0}})=n\) and we have \(d\leq k(\mu_{0})\) for almost all \(\mu_{0}\), but not necessarily \(d=k(\mu_{0})\)2. However, the data are exceptionally unlikely to lie on orbits with high codimension \(k(\mu)>d\), and we are going to assume that the effective number of dimensions of \(\alpha\) is \(n-d\). If the data was obtained from one of the orbits with high codimension, one would amend formula (26) below using the information about that exceptional orbit. Note that our method preserves the general form of coadjoint orbits in all cases, whether the orbit is exceptional or not. This fact could be useful in data-based computations of exceptional orbits. It is an interesting question which we will address in the follow-up work._

Footnote 2: Note that there are examples of Lie algebras whose generic coadjoint orbits have codimension strictly bigger than the number of independent Casimirs.

Equations (22) are linear differential equations in \(\mu\), and the solutions of these equations can be (in principle) found in explicit form. This solution can be written in compact form as

\[\mathbb{T}(t,\alpha)\mu_{0}=e^{\mathbb{M}(\alpha)t}\mu_{0}\,. \tag{23}\]

The mappings \(\mathbb{T}(t,\alpha)\) defined in (23) satisfy all the requirements as the building blocks for the neural networks.

Let us denote the mapping \(\mathbb{T}_{a}(t,\alpha^{a})\), \(a=1,...,n\), to be the map originating from the \(a\)-th component of \(\alpha\) to have the value of \(\alpha^{a}\) and all other values being zero. We divide each time step into \(n\) substeps \(h_{a}\) with \(\sum_{a}h_{a}=h\). Notice that the number of steps is \(n\) and not \(n-d\), where \(d\) is the number of Casimirs, as we explain later.

These mappings satisfy the following conditions:

1. The mappings \(\mathbb{T}_{a}(h_{a},\alpha^{a})\) are Poisson, see (20), for the Lie-Poisson bracket,
2. The mappings conserve all Casimirs of the system, 3 Footnote 3: And, in fact, these maps also preserve all types of coadjoint orbits, _i.e._, they accurately represent the coadjoint action \(\mathrm{Ad}_{g}^{*}\mu\). This information would be useful if we were to compute exceptional orbits, which we will not do here.
3. Any two points close enough to each other on the same Casimir surface can be connected using a combination of mappings \(\mathbb{T}_{a}(h_{a},\alpha^{a})\).

The procedure of constructing LPNets is as follows:

1. Find explicit solutions for (23) by solving equations (22).
2. Define \(\bar{\alpha}=(\alpha^{1},\ldots,\alpha^{n})\) and \[\mathbb{T}(\bar{\alpha})=\mathbb{T}_{n}(h_{n},\alpha^{n})\circ\mathbb{T}_{n-1 }(h_{n-1},\alpha^{n-1})\circ\ldots\circ\mathbb{T}_{1}(h_{1},\alpha^{1})\,.\] (24) For the set of \(N\) data pairs \((\mu_{i}^{0},\mu_{i}^{f})\), \(i=1,\ldots N\) set up a minimization procedure to find the set of numbers \(\bar{\alpha}\) for each \(\mu_{i}^{0}\), minimizing the "individual" square loss \[\bar{\alpha}_{i}=\arg\min\left|\mathbb{T}(\bar{\alpha}_{i})\mu_{i}^{0}-\mu_{i }^{f}\right|^{2}\] (25) In most of the examples we consider here, we can find \(\bar{\alpha}\) analytically; however, in more general situations, one could also utilize a root-finding or a minimization procedure finding \(\bar{\alpha}_{i}\) for every pair of data points.
3. Since \(\bar{\alpha}\) are only defined up to a vector normal to the Casimir surface, in order to make the data consistent, we need to project out the appropriate components of the Casimirs. We will need to find coefficients \(p_{j}\), \(j=1,\ldots,d\) such that the projection of \(\bar{\alpha}\) on the gradients of the Casimirs vanishes: \[\bar{\alpha}_{i}\rightarrow\bar{\alpha}_{i}-\sum_{j=1}^{d}p_{j}\left.\frac{ \partial C_{j}}{\partial\mu}\right|_{\mu=\mu_{i}^{0}}\,,\,\,\mathrm{with}\, \left\langle\bar{\alpha}_{i}\,,\,\,\frac{\partial C_{j}}{\partial\mu}\right|_ {\mu=\mu_{i}^{0}}\!\right\rangle=0\,,j=1,...,d,\] (26) where \(C_{j}\), \(j=1,...,d\) are a set of \(d\) independent Casimir functions.
4. Create a neural network approximating the mapping \(\mu_{i}\rightarrow\bar{\alpha}_{i}\). This function will be denoted as \(\bar{\alpha}=NN(\mu)\).
Since the operators \(\mathbb{T}_{a}\), in general, do not commute, the composition order should be fixed ahead of time and also be the same for all data pairs. A different choice of the composition order will lead to a different (but of course equivalent) set of parameters \(\alpha\). The choice of the composition order of \(\mathbb{T}_{a}\) has to be preserved in the prediction step as well. Also, we do not consider different time sub-steps \(h_{a}\) here, taking them all to be equal. It is possible that the choice of \(h_{a}\) can be made to improve accuracy or convergence properties of the scheme.

The solution starting at a given initial condition \(\mu_{0}\) can then be evaluated using the neural network. If \(\mu=\mu_{j}\) at the time step \(j\), then the value of the solution at the next step is computed as

\[\mu_{j+1}=\mathbb{T}(\bar{\alpha}_{j}^{\rm est})\mu_{j},\quad\bar{\alpha}_{j}^ {\rm est}=NN(\mu_{j})\,. \tag{27}\]

Note that we never compute the actual Hamiltonian, its gradients, or the equations of motion. Instead, we just compute the composition of Poisson transformations reproducing the dynamics in phase space of some unknown Poisson system - with the known Lie-Poisson bracket.

We will now apply this procedure to several particular examples of physical systems which have a Lie-Poisson bracket. We shall call this method Local _LPNets_, or simply LPNets, as the exponential representation of the map is only defined locally. The advantage of this method is that when operating on Lie groups, we are guaranteed to be able to reach any point locally with an exponential map, so the completeness is automatic. However, this mathematical simplicity has to be compensated by the necessity to apply a neural network to learn the mappings from data points in the neighborhood of a trajectory. After our description of LPNets, we develop a more complex procedure which derives Lie-Poisson activation modules that directly extend the work of (Bajars, 2023; Jin et al., 2022, 2020). We call these methods global Lie-Poisson networks, or G-LPNets. Using the example of a rigid body motion, we show that G-LPNets provide a promising simple, accurate and highly computationally effective neural network.

A more general derivation performed in the coordinate-free form and in the general language of modern geometric approach is presented in A.

## 5 Test cases for LPNets

We choose to use the same test cases as in (Bajars, 2023; Jin et al., 2022), except for an additional test for \(SE(3)\) (the underwater vehicle). No detailed comparisons are performed between our results and those in (Bajars, 2023; Jin et al.,
2022). This is because the accuracy of each method depends strongly on the structure of the neural network, the accuracy and distribution of data used in training, and specific learning procedures employed, which prevents a fair comparison between methods. Instead, we will outline the performance of LPNets on these problems, and contrast with general results from (Bajars, 2023; Jin et al., 2022). In all cases, LPNets conserve the Casimir functions to machine-precision, unlike the other approaches.

In what follows, we will develop the prediction of parameters for LPNets using a dense Neural Network structure. Naturally, such a dense network will require appropriate number of data points to avoid overfitting. In Section 6, we show how to reduce the size of that network and achieve the accuracy in the whole phase space with a particular structure of a network which we will call G-LPNets.

### Rigid body dynamics

_Motion of a rigid body as dynamics on the Lie group \(S\,O(3)\)_. The (Lie-)Poisson bracket, the Hamiltonian, and the corresponding equations of motion for momenta \(\mathbf{\Pi}\) (measured from the body frame of reference), are (Holm et al., 2009):

\[\begin{split}\{F,G\}&=-\mathbf{\Pi}\cdot\left(\frac{ \partial F}{\partial\mathbf{\Pi}}\times\frac{\partial G}{\partial\mathbf{\Pi}}\right) \\ H(\mathbf{\Pi})&=\frac{1}{2}\mathbf{\Pi}\cdot\mathbb{I}^{- 1}\mathbf{\Pi}\\ \dot{\mathbf{\Pi}}&=-\mathbb{I}^{-1}\mathbf{\Pi}\times\mathbf{ \Pi}\,.\end{split} \tag{28}\]

Suppose there is a sequence of pairs of initial and final points of transformation \((\mathbf{\Pi}_{i}^{0},\mathbf{\Pi}_{i}^{f})\), \(i=1,\ldots,N\), coming from some information that we call _ground truth_. If that sequence comes from a single trajectory of length \(N\), _i.e._ has the form \(\mathbf{\Pi}_{0},\mathbf{\Pi}_{1},\ldots,\mathbf{\Pi}_{N}\) we take \(\mathbf{\Pi}_{i}^{0}=\mathbf{\Pi}_{i}\), \(\mathbf{\Pi}_{0}^{f}=\mathbf{\Pi}_{i+1}\). However, our method does not explicitly assume the existence of a single trajectory for learning.

To find the Poisson map approximating the motion of the system at every time step \(i\), let us consider Hamiltonians \(H_{i}\) linear in momenta, _i.e._ having the form \(H_{i}(\mathbf{\Pi})=\mathbf{A}_{i}\cdot\mathbf{\Pi}\), where \(\mathbf{A}_{i}\) is some unknown constant vector that is different for every pair of points. The Hamiltonian flow generated by that Hamiltonian is described by

\[\dot{\mathbf{\Pi}}=-\frac{\partial H_{i}}{\partial\mathbf{\Pi}}\times\mathbf{\Pi}=- \mathbf{A}_{i}\times\mathbf{\Pi}\,. \tag{29}\]

Note that we can also consider a Hamiltonian of the form \(H(\mathbf{\Pi})=f(\mathbf{A}\cdot\mathbf{\Pi})\). The dynamics induced by this Hamiltonian preserve the quantity \(\mathbf{A}\cdot\mathbf{\Pi}\), as one
can see from (29). Thus, that extension corresponds to a simple redefinition of \(\mathbf{A}\) by scaling and does not bring extra insight into the problem.4 Thus, the motion defined by (29) is a rotation of a vector \(\mathbf{\Pi}\) about the axis \(\mathbf{A}_{i}\) with a constant angular velocity. The flow preserves the Lie-Poisson bracket since it is a Hamiltonian flow with the same bracket. The dynamics (28) also preserve the Casimir \(C(\mathbf{\Pi})=|\mathbf{\Pi}|^{2}\) exactly. After the time \(t=h\), the dynamics (28) rotates the vector of angular momentum \(\mathbf{\Pi}\) by an angle \(\phi(\mathbf{A}_{i})=|\mathbf{A}_{i}|h\) around the axis \(\mathbf{n}_{\mathbf{A}_{i}}=\mathbf{A}_{i}/|\mathbf{A}_{i}|\). In other words, if \(\mathbb{R}(\mathbf{n},\phi)\) is the matrix of rotation around the axis \(\mathbf{n}\) by the angle \(\phi\), then (28) transforms the momentum as \(\mathbf{\Pi}\rightarrow\mathbb{R}(\mathbf{n}_{\mathbf{A}_{i}},\phi(\mathbf{ A}_{i}))\mathbf{\Pi}\).

Footnote 4: Of course, \(\mathbf{A}\cdot\mathbf{\Pi}\) is not a constant for the general system (28) - that is only true for the particular choice of the Hamiltonian \(H(\mathbf{\Pi})=f(\mathbf{A}\cdot\mathbf{\Pi})\).

Everywhere in this paper, the ground truth is obtained by BDF integrator in Python's _Scipy_ package, with relative and absolute tolerances being set at \(10^{-13}\) and \(10^{-14}\), respectively. Note that this numerical method is not expected to conserve Casimirs so our method will actually be more precise than the ground truth. We note that while Lie-Poisson integrators preserving Lie-Poisson structure do exist (Marsden et al., 1999), the precise implementation needs to be rederived in an explicit form for each particular Lie-Poisson system. We found it to be more appropriate to use a high accuracy algorithm that is common to all problems for a fair comparison, rather than build an algorithm that is tailored to every particular problem. The accuracy of the BDF algorithm is more than sufficient for our ground truth calculation and providing comparisons with the previous works. Thus, to be consistent, we used a high precision non-symplectic integrator for all cases, carefully checking its accuracy in all applications.

_Data preparation._ We should consider three angles of rotations. Given a sequence of begin and end pairs of momenta (\(\mathbf{\Pi}_{i}^{0},\mathbf{\Pi}_{i}^{f}\)), we compute \(\mathbf{A}_{i}\) as a function of \(\mathbf{\Pi}_{i}^{0}\) as follows:

\[\mathbf{A}_{i}=\frac{1}{h}\mathbf{\Pi}_{i}^{0}\times\mathbf{\Pi}_{i}^{f} \tag{30}\]

and the angle \(\theta_{i}\) of rotations from vector \(\mathbf{\Pi}_{0}^{i}\) and \(\mathbf{\Pi}_{f}^{i}\) as the shortest motion along the sphere \(\mathbf{\Pi}=\)const. The cross product contains information for both the direction normal to both \(\mathbf{\Pi}_{i}^{0}\) and \(\mathbf{\Pi}_{i}^{f}\), and the angle of rotation, as described above. The factor \(1/h\) is introduced to normalize the output data to be of order 1.

We could of course get \(\mathbf{A}_{i}\) by finding the match of three subsequent rotations around, say, Euler angles, and subtracting the corresponding rotation about the axis \(\mathbf{\Pi}_{i}^{0}\) or \(\mathbf{\Pi}_{i}^{f}\), applying equation (25) directly. However, the notation of vector
cross product, only available for \(SO(3)\), provides a simple and efficient alternative to this more complex procedure.

_Neural network._ Make a standard neural network learning from the data for the mapping \(\mathbf{\Pi}\to\mathbf{A}\). The neural network will have the three components of \(\mathbf{\Pi}\) as inputs and the three components of \(\mathbf{A}\) as outputs. Here and everywhere else, we utilize the package _Tensorflow5_. In the results shown in Figure 1, the learning is done on \(N=1000\) pairs produced by the single trajectory originating at \(\mathbf{\Pi}_{0}=(1/\sqrt{2},-1/\sqrt{2},1)\), with the interval between trajectory points given by \(h=0.1\).

Footnote 5: [https://www.tensorflow.org](https://www.tensorflow.org)

The neural network has three hidden layers of 16 neurons with the sigmoid activation function, with 659 trainable parameters. Out of 1000 pairs as the input data, 80% are used for training and 20% for evaluation, with the loss measured as the mean square discrepancy between \(\mathbb{R}(\mathbf{n}_{\mathbf{A}_{i}},\phi(\mathbf{A}_{i}))\mathbf{\Pi}_{i}^{0}\) and \(\mathbf{\Pi}_{i}^{f}\). Adam optimization algorithm is used with the learning rate starting at \(10^{-3}\). The loss and validation loss reach the values of approximately \(3.8\cdot 10^{-9}\) and \(4.1\cdot 10^{-9}\) after \(10^{5}\) epochs.

_Prediction._ The predicted trajectory start \(\mathbf{\Pi}_{0}\) is taken to coincide with the endpoint of the learning trajectory, with the values of \(\mathbf{\Pi}_{0}\simeq(0.43,1.33,-0.21)\). On each step, once \(\mathbf{\Pi}_{j-1}\) is known, Neural Network creates the prediction for the rotation axis \(\mathbf{A}_{j}/h\) and the angle \(\phi_{j}=|\mathbf{A}_{j}|/h\). That prediction of the neural network is used to produce the mappings \(\mathbf{\Pi}_{j}=\mathbb{R}(\mathbf{A}_{j},\phi_{j})\mathbf{\Pi}_{j-1}\), where \(j=1\dots m\) after the desired number of steps. That prediction by LPNets is compared with the ground truth prediction obtained by high accuracy ODE solver as described above. We perform 10000 time steps to reach \(t=1000\) and present the results in Figure 1 and Figure 2. Note that for clarity, only the first 2000 time steps up to \(t=200\) are shown in the individual momenta plots, the left panel of Figure 1. The right panel of that Figure shows that all available data coincide perfectly.

In Figure 2, we present the results for the conservation of the Hamiltonian \(H(\mathbf{\Pi})\) and the Casimir \(C(\mathbf{\Pi})=|\mathbf{\Pi}|^{2}\). The Hamiltonian is preserved to about 0.01% relative accuracy. The Casimir in the ground truth solution is preserved to about \(10^{-9}\) accuracy. In the LPNets solution, the Casimir is preserved to machine precision, far exceeding possible accuracy for ground truth for this variable. On the right side of this panel, we plot the error of the solution as the \(L_{2}\) norm of deviation between two solutions (ground truth and LPNets) for each \(t\). The deviation is growing roughly linearly in time to the values of about 0.5 after the time \(t=1000\).
This may come as a surprise given the excellent agreement on the right hand side of the Figure 1 for all \(t\). This phenomenon has a simple explanation: high accuracy of conservation laws presented on Figure 2 forces the solution to exist on the intersection of \(H=\)const (an ellipsoid) and \(C=\)const (a sphere) with the high accuracy.

Comparison with the previous literatureThe case of rigid body motion was considered in (Bajars, 2023). The appropriate comparison is the learning of the single trajectory of the rigid body dynamics contained in SS4.2.3 of that paper. The error of a single trajectory (ground truth vs. a solution obtained by the neural network solution) is of the same order as our results. The error in Hamiltonian is somewhat better in our case, with the relative error being about \(10^{-4}\) (\(0.01\%\)) vs. \(0.008\) in (Bajars, 2023). However, that number depends on the particular realization of both neural networks and learning data, as we outlined above. The value of the Casimir \(C(\mathbf{\Pi})=|\mathbf{\Pi}|^{2}\) is conserved to machine precision in our case, whereas it would follow the general accuracy of computations in the previous work on the subject.

Learning general dynamics of the rigid bodyIn the above calculation, we have followed the method of (Bajars, 2023; Jin et al., 2022) and learned the dynamics continuing a _single_ trajectory. However, it is also possible to extend the LPNets to learn several trajectories simultaneously, and predict the dynamics of the trajectory the method has not seen. We take the initial condition \(\mathbf{\bar{\Pi}}_{0}=(1/\sqrt{2},-1/\sqrt{2},1)\) and compute 20 ground truth trajectories with the initial con

Figure 1: Left: Results of LPNets applied to the motion of a rigid body (red) versus ground truth (blue) for the individual momenta. Right: Parametric plot of the momenta in the phase space. The results are visually indistinguishable.
ditions \(\mathbf{\Pi}_{0}^{j}=\bar{\mathbf{\Pi}}_{0}+\epsilon_{j}\), where \(\epsilon_{j}\) is a uniformly distributed random variable in the cube \([-0.1,0.1]\times[-0.1,0.1]\times[-0.1,0.1]\). Each trajectory generates 200 pairs of ground truth mapping between the momenta at the neighboring points, a total of 4000 data points. A neural network is constructed with a similar structure as the one used for learning a single trajectory, having three inner layers with 32 neurons, each having a sigmoid activation function, and the total of 2339 trainable parameters. The neural network is trained using Adam algorithm with initial time step of \(10^{-3}\) decaying exponentially to \(10^{-4}\), over \(10^{5}\) epochs. The final training and validation losses are slightly below \(10^{-6}\) and \(10^{-5}\) respectively, after 100,000 epochs.

Trajectory prediction using LPNetsA trajectory is then constructed with the initial conditions \(\bar{\mathbf{\Phi}}_{0}=\bar{\mathbf{\Pi}}_{0}\) iterating over 10000 steps up to the time \(t=1000\). Even though all the learning data were taken a finite distance from this solution, the LPNets faithfully reproduces the ground truth, as shown in Figure 3. While the solutions were computed until \(t=1000\), the left panel of Figure 3 only illustrates the evolution of momenta until \(t=200\) for clarity. This example of \(SO(3)\) shows that our method is capable of learning multiple trajectories and understanding the general Poisson dynamics. Of course, one has to keep in mind that in order to achieve good global accuracy, one needs to have quite a dense covering of the configuration manifold with the data points for learning, a task that may be difficult in many dimensions. A compromise presented here that is feasible to implement

Figure 2: Left: Conservation of the Hamiltonian \(H\) (top) and the Casimir \(C\) (bottom), comparing the results of LPNets (red) and ground truth (blue). Notice that LPNets conserve the Casimir exactly (to machine precision) and thus substantially exceeds the ground truth in the conservation of Casimirs. Right: The discrepancy between the results of LPNets and the ground truth. The discrepancy comes mostly from time mismatch, whereas the amplitude of oscillations is conserved with high precision.
is to consider a few trajectories in the neighborhood of the desired trajectory for learning. To illustrate that point, we show the 3D plot of the trajectories and the corresponding data for learning on the right panel of Figure 3. Note that all the data were used on the right panel, and there is no visible deviation between the ground truth and the solutions in 3D. We present the accuracy of the results compared to the ground truth and the preservation of the conserved quantities (Energy and Casimir) on Figure 4. The conservation of the Hamiltonian is satisfied with the relative accuracy of about 0.25% (0.001 in absolute accuracy), and the Casimir is conserved in our case to the absolute accuracy of about \(10^{-11}\), several orders of magnitude better than the ground truth solution.

### Extended pendulum case

In order to compare our results directly with the Poisson Neural Networks developed in (Jin et al., 2022), we consider the extended pendulum test case from that paper. This example shows that our method extends beyond the Lie-Poisson case, as long as an explicit integration for the Poisson equations for linear Hamiltonians is available. This is case here since the Poisson tensor \(\mathbb{B}(\mathbf{y})\) is affine in \(\mathbf{y}\).

Consider a standard pendulum of length 1 and mass 1, with the Hamiltonian

Figure 3: Left: Same results as in Figure 1, with the exception that the Neural Network is trained on several trajectories with random initial conditions different from the desired trajectory. Starting point is taken to be \(\boldsymbol{\Pi}_{0}=(1/\sqrt{2},-1/\sqrt{2},1)\). Right: Trajectories in 3D space, similar to the presentation on the right panel of Figure 1. In addition, the trajectories used for data learning are shown in black, and the blue dot indicates the starting point. Trajectories from LPNets are shown in red, learning trajectories in black, and trajectories considered ground truth are in blue.
\(H=\frac{1}{2}p^{2}-\cos q\). The equations of motion are written as

\[\dot{q}=\frac{\partial H}{\partial p}=p\,,\quad\dot{p}=-\frac{\partial H}{ \partial q}=-\sin q\,. \tag{31}\]

The paper (Jin et al., 2022) then introduces an extra variable \(c\) with equation \(c=\)const, extending the system to the three-dimensional space

\[\frac{d}{dt}\left(\begin{array}{c}p\\ q\\ c\end{array}\right)=\left(\begin{array}{c}-\sin q\\ p+c\\ 0\end{array}\right)=\left(\begin{array}{ccc}0&-1&0\\ 1&0&0\\ 0&0&0\end{array}\right)\nabla_{(p,q,c)}\widetilde{H}\,, \tag{32}\]

with the new Hamiltonian \(\widetilde{H}=\frac{1}{2}p^{2}-\cos q+pc\). The paper (Jin et al., 2022) then makes a transformation of variables

\[\begin{split}&(p,q,c)=\theta(u,v,r)=(u,v,r-u^{2}-v^{2})\,,\\ &(u,v,r)=\theta^{-1}(p,q,c)=(p,q,p^{2}+q^{2}+c)\,.\end{split} \tag{33}\]

In the new variables \((u,v,r)\), the equations of motion become

\[\frac{d}{dt}\left(\begin{array}{c}u\\ v\\ r\end{array}\right)=\left(\begin{array}{ccc}0&-1&-2v\\ 1&0&2u\\ 2v&-2u&0\end{array}\right)\nabla_{(u,v,r)}K \tag{34}\]

Figure 4: Left: Conservation of the Hamiltonian \(H\) (top) and the Casimir \(C\) (bottom), comparing the results of LPNets (red) and ground truth (blue), for the case of global dynamics. Again, LPNets conserves the Casimir exactly (to machine precision) and thus substantially exceeds the ground truth in the conservation of Casimirs. Right: The discrepancy between the results of LPNets and the ground truth.
with the new Hamiltonian \(K(u,v,r)=\frac{1}{2}u^{2}-\cos v+ur-u^{3}-uv^{2}\). The system (34) is Poisson with the corresponding bracket

\[\{F,H\}=(\nabla_{\mathbf{y}}F)^{T}\cdot\mathbb{B}(\mathbf{y})\cdot\nabla_{ \mathbf{y}}H\,,\quad\mathbb{B}(\mathbf{y}):=\left(\begin{array}{ccc}0&-1&-2v \\ 1&0&2u\\ 2v&-2u&0\end{array}\right)\,, \tag{35}\]

where we have denoted \(\mathbf{y}=(u,v,r)^{T}\). The matrix \(\mathbb{B}\) defined in (35) is degenerate, and

\[C(\mathbf{y})=r-u^{2}-v^{2}=y_{3}-y_{1}^{2}-y_{2}^{2}\,, \tag{36}\]

which is just \(c\) in the old variables \((p,q,c)\), is a Casimir of the bracket (35). Indeed, one can readily check that \(\mathbb{B}(\mathbf{y})\cdot\nabla_{\mathbf{y}}C(\mathbf{y})=\mathbf{0}\). Moreover, one also checks that \(\mathbb{B}(\mathbf{y})\) only has a single eigenvalue of \(0\), so \(C(\mathbf{y})\) defined in (36) is the only Casimir. In order to apply the method of LPNets, we take the test Hamiltonian linear in \(\mathbf{y}\) as \(H(\mathbf{y})=\boldsymbol{\alpha}\cdot\mathbf{y}\). For test Hamiltonians of that type, the equations of motion become

\[\dot{\mathbf{y}}=\left(\begin{array}{c}-\alpha_{2}-2y_{2}\alpha_{3}\\ \alpha_{1}+2y_{1}\alpha_{3}\\ 2\alpha_{1}y_{2}-2\alpha_{2}y_{1}\end{array}\right)\,. \tag{37}\]

There are three test Hamiltonians to consider, \(H_{a}(\mathbf{y})=\alpha_{a}y_{a}\), \(a=1,2,3\) (no sum). Then, the equations of motion are

\[\begin{cases}H_{1}=\alpha_{1}y_{1}\Rightarrow\,\dot{y}_{1}=0\,,\dot{y}_{2}= \alpha_{1}\,,\dot{y}_{3}=2\alpha_{1}y_{2}\\ H_{2}=\alpha_{2}y_{2}\Rightarrow\,\dot{y}_{1}=-\alpha_{2}\,,\dot{y}_{2}=0\,, \dot{y}_{3}=-2\alpha_{2}y_{1}\\ H_{3}=\alpha_{3}y_{3}\Rightarrow\,\dot{y}_{1}=-2y_{2}\alpha_{3}\,,\dot{y}_{2}= 2y_{1}\alpha_{3}\,,\dot{y}_{3}=0\,.\end{cases} \tag{38}\]

Equations (38) are easily solved explicitly. Each Hamiltonian \(H_{a}=\alpha_{a}y_{a}\) leads to explicit expressions for an affine transformation \(\mathbf{T}(t,\alpha_{a},\mathbf{y}_{0})\) of the initial condition \(\mathbf{y}_{0}\) to the final solution after time \(t\):

\[H_{a}=\alpha_{a}y_{a}\,\Rightarrow\,\mathbf{y}=\mathbf{T}_{a}(t,\alpha_{a}, \mathbf{y}_{0}) \tag{39}\]
with the transformations \({\bf T}_{a}\) given by

\[{\bf T}_{1}(t,\alpha_{1},{\bf y}_{0})= \left(\begin{array}{c}y_{1}(0)\\ y_{2}(0)+\alpha_{1}t\\ y_{3}(0)+2t\alpha_{1}y_{2}(0)+\alpha_{1}^{2}t^{2}\end{array}\right) \tag{40}\] \[{\bf T}_{2}(t,\alpha_{2},{\bf y}_{0})= \left(\begin{array}{c}y_{1}(0)-\alpha_{2}t\\ y_{2}(0)\\ y_{3}(0)-2t\alpha_{2}y_{1}(0)+\alpha_{2}^{2}t^{2}\end{array}\right)\] \[{\bf T}_{3}(t,\alpha_{3},{\bf y}_{0})= \left(\begin{array}{c}y_{1}(0)\cos(2\alpha_{3}t)-y_{2}(0)\sin(2 \alpha_{3}t)\\ y_{1}(0)\sin(2\alpha_{3}t)+y_{2}(0)\cos(2\alpha_{3}t)\\ y_{3}(0)\end{array}\right)\,.\]

_Data preparation, Part 1: finding the transformations._ Suppose we have pairs of solution \(({\bf y}_{i},{\bf y}_{i+1})\) that are obtained from snapshots of a single or several trajectories of equation (34); the time difference between the snapshots is \(\Delta t=h\). We separate the interval \(\Delta t=h\) between the snapshots into three equal sub-intervals (although other divisions of the time interval are also possible). Then, for each initial point of the pair \({\bf y}_{i}\) we are looking for the sequence of parameters \((\alpha_{1},\alpha_{2},\alpha_{3})\) defining the transformations \({\bf T}_{i}\) as in (40), such that

\[{\bf y}_{i}^{1} ={\bf T}_{1}\left(\frac{h}{3},\alpha_{1},{\bf y}_{i}\right)\,, \tag{41}\] \[{\bf y}_{i}^{2} ={\bf T}_{2}\left(\frac{h}{3},\alpha_{2},{\bf y}_{i}^{1}\right)\,,\] \[{\bf y}_{i}^{3} ={\bf T}_{3}\left(\frac{h}{3},\alpha_{3},{\bf y}_{i}^{2}\right)= {\bf y}_{i+1}\,.\]

In other words, \(\overline{\boldsymbol{\alpha}}_{i}=(\alpha_{1},\alpha_{2},\alpha_{3})({\bf y }_{i})\) should be such that the three transformations \({\bf T}_{1,2,3}\), performed in the corresponding order, map the initial snapshot of the pair \({\bf y}_{i}\) to the final snapshot of the pair \({\bf y}_{i+1}\).

For data that is not precise, the matching of \(\overline{\boldsymbol{\alpha}}_{i}\) to the data can be accomplished numerically using a gradient descent method, with further removal of components of the Casimir gradient \(\nabla C\). Indeed, each set of parameters \(\overline{\boldsymbol{\alpha}}_{i}=(\alpha_{1},\alpha_{2},\alpha_{3})({\bf y }_{i})\) determined in the previous step, is only defined up to the value of \(\nabla C\). Since \(\nabla C\) is a zero eigenvector of the matrix \(\mathbb{B}({\bf y})\), changing

\[\boldsymbol{\alpha}_{i}\rightarrow\boldsymbol{\alpha}_{i}+k\nabla C \tag{42}\]
for arbitrary \(k\) does not change the result of composition of transformations \(\mathbf{T}_{3}\circ\mathbf{T}_{2}\circ\mathbf{T}_{1}\). We thus choose \(k\) in (42) such that the projection of \(\boldsymbol{\alpha}\) on \(\nabla C\) vanishes, _i.e._, take

\[\boldsymbol{\alpha}_{i}^{*}=\boldsymbol{\alpha}_{i}-\nabla C\frac{\boldsymbol{ \alpha}_{i}\cdot\nabla C}{|\nabla C|^{2}}\,. \tag{43}\]

However, in the case the training data is _exact_, for example, is obtained from a numerical simulation with a very high accuracy, we can solve equations for \(\overline{\boldsymbol{\alpha}}_{i}\) analytically. For shortness, we denote the components of \(\mathbf{y}_{i}\) as \((y_{1},y_{2},y_{3})\), dropping the index \(i\), and the components of \(\mathbf{y}_{i+1}\) by \((y_{1,f},y_{2,f},y_{3,f})\), since \(\mathbf{y}_{i+1}\) represents the final value of the interval \(t\in(t_{i},t_{i+1})\). Let us notice that the sequential application of \(\mathbf{T}_{1}\) and \(\mathbf{T}_{2}\) after time \(\Delta t=h/3\) on each sub-step gives the intermediate values

\[\begin{split} y_{1}&\to y_{1}^{*}=y_{1}-\alpha_{2} \Delta t\,,\\ y_{2}&\to y_{2}^{*}=y_{2}+\alpha_{1}\Delta t\,,\\ y_{3}&\to y_{3}^{*}=y_{3}+(\alpha_{1}-\alpha_{2}) \Delta t+(\alpha_{1}^{2}+\alpha_{2}^{2})\Delta t^{2}\,.\end{split} \tag{44}\]

On the third sub-step, \(y_{3}\) doesn't change so \(y_{3}^{*}=y_{3,f}\). Also, at that sub-step, the application of \(\mathbf{T}_{3}\) just induces a rotation of \((y_{1}^{*},y_{2}^{*})\) by the angle \(2\alpha_{3}\Delta t\), so the compatibility conditions to match the data points precisely is

\[y_{3}^{*}=y_{3,f}\,,\quad(y_{1}^{*})^{2}+(y_{2}^{*})^{2}=y_{1,f}^{2}+y_{2,f}^{ 2}\,. \tag{45}\]

Using (44), we can rewrite (45) as

\[y_{3,f}-y_{1,f}^{2}-y_{2,f}^{2}=y_{3}-y_{1}^{2}-y_{2}^{2}\,, \tag{46}\]

so the solution for \(\overline{\boldsymbol{\alpha}}_{i}\) giving exact match between the data points can be found if and only if the Casimir is exactly the same for all values of the data. If the data contains noise, an optimization procedure must be introduced searching for an optimal value of \(\overline{\boldsymbol{\alpha}}_{i}\) at every step. In this paper, we assume that the learning data is exact. In that case, we can set \(\alpha_{3}=0\) on every time step and match the values of the components \(y_{1}\) and \(y_{2}\); the matching of the component \(y_{3}\) is done automatically due to the conservation of the Casimir.

\[\alpha_{1}=\frac{3}{h}\left(y_{2,i+1}-y_{2,i}\right)\,,\quad\alpha_{2}=-\frac {3}{h}\left(y_{1,i+1}-y_{1,i}\right)\,,\quad\alpha_{3}=0\,. \tag{47}\]

_Learning procedure: Neural Network approximation._ The inputs \(\mathbf{y}_{i}\) and outputs \(\overline{\boldsymbol{\alpha}}_{i}\) computed from (47) are used to learn the mapping \(\overline{\boldsymbol{\alpha}}(\mathbf{y})\) using a Neural Network. In order to match (Jin et al., 2022), we use three trajectories with initial conditions
\(\mathbf{y}_{0}=(0,1,1^{2})\), \((0,1.5,1.5^{2}+0.1)\) and \((0,2,2^{2}+0.2)\), with the time step of \(h=0.1\). We use the number of points \(N=1000\) with half of these points used for training and half for validation. The ground truth solution is computed using Python's _Scipy odeint_ routine with _BDF_ (Backward Differentiation Formula), with the tolerances, relative and absolute, set at \(10^{-13}\) and \(10^{-14}\), respectively.

Evaluation and dynamics using Neural NetworkJust as in (Jin et al., 2022), we compute the trajectories starting at the last point of the trajectory \(\mathbf{y}_{N}\) for 1000 points, with half of these points used for training and half for validation. We use a Sequential Keras program, with three inner layers of neurons with sigmoid activation functions. Each layer has the breadth of 16 neurons with the total number of tunable parameters being 642. The network uses the Adam optimizer with the learning rate of \(10^{-3}\), exponentially decaying to \(10^{-5}\) and Mean Square Error loss, computed for 50000 epochs. The resulting MSE achieved is between \(10^{-7}\) and \(10^{-8}\). The results of the simulations are shown in Figure 5. On the left side of this Figure, we show the match between the "exact" solution and the solution obtained by LPNets (notice that the "ground truth" solution is still obtained using a numerical method with a given accuracy). On the right hand side of that figure, we present a three-dimensional plot of the phase space \((u,v,r)=(y_{1},y_{2},y_{3})\) comparing the numerical results in blue with the LPNets results in red. In Figure 6, we show the conservation of the Casimir for all three initial conditions. Note that the ground truth numerical solution (blue line) only conserves the Casimir \(C=r-u^{2}-v^{2}\) up to the accuracy of calculation, accumulating the error of the or

Figure 5: Results of simulation of equations (34) and the corresponding solution of using LPNets procedure. Left panel: three components of the solution \(y_{1}=u\), \(y_{2}=v\) and \(y_{3}=r\) versus time. Right panel: phase space plot of the solution. Blue line: high-precision numerics taken as the exact solution; red line: LPNets.
der \(10^{-9}\) after \(t=100\). In contrast, the Casimir in LPNets (red line) is preserved to machine precision by the very nature of transformations performed on each time step. On the right panel of that Figure, we present the accuracy of the evolution of Energy in LPNets. As one can see, the energy conservation by LPNets is quite satisfactory, yielding the relative error of about \(1\%\) or less for all cases.

Finally, in Figure 7, we present the results for the discrepancy between the ground truth and LPNets solutions. Even though the discrepancy grows, it is mostly due to the fact that the time between the ground truth and the solution is slightly mismatched, which explains why energy is conserved to much higher precision than the solution itself. The accuracy is still quite good and the solution obtained by the LPNets is virtually indistinguishable from the ground truth on the left side of the Figure 5.

Finally, we would like to emphasize that a system of three inner layers having a width of \(16\) in each layer should be viewed as very compact. Higher accuracy can be achieved with more data points and correspondingly wider or deeper network, or having more insights into the structure of the framework which we have assumed to be completely dense. Achieving these efficiencies is an interesting challenge that we will consider in the future.

Comparison with previous results for extended pendulum systemThe visual agreement of the solutions and the ground truth is similar to the one presented in (Jin et al., 2022) for this system. The conservation of the Casimir is not presented in (Jin et al., 2022). We interrogated the data produced by the code made available in

Figure 6: Left: Conservation of the Casimir \(C=r-u^{2}-v^{2}\) in the solutions of equations (34) obtained by the high precision numerics (blue line) and the corresponding solution of using LPNets procedure (red line). Even though the precision of numerics is \(10^{-11}\), LPNets is substantially more precise, as it achieves machine precision of Casimir conservation on each time step. Right: Conservation of the energy of the system for all three cases. The relative accuracy in the conservation of energy is about \(0.5-1\%\).
Figure 7: The discrepancy between the ground truth and the solutions provided by LPNets for the extended pendulum case.
(Jin et al., 2022) and found that the relative errors for Casimir in the transformed Lie-Darboux coordinates are quite small, between \(10^{-7}\) and \(10^{-6}\). The errors in Casimir in the original coordinates \((u,v,r)\) for the parameters presented in (Jin et al., 2022), although still quite small, are several orders of magnitude larger than the transformed coordinates. This could be attributed to the fact that the inverse of Lie-Darboux transformation is not being computed sufficiently accurately. This accuracy of PNNs can certainly be improved by appropriate modification of the original PNN network from (Jin et al., 2022), which is beyond the scope of the paper. In any case, our method conserves the Casimir exactly for the original system with no necessity to search for the Lie-Darboux transformation to the canonical coordinates.

### A particle in a magnetic field

To compare with the second test case computed in (Jin et al., 2022), we study a particle of mass \(m\) and charge \(q\) moving in a magnetic field \(\mathbf{B}(\mathbf{x})\). We assume that the motion of the particle is in \(\mathbf{x}\in\mathbb{R}^{3}\), and the relevant variables are the particle position \(\mathbf{x}=(x_{1},x_{2},x_{3})\) and its' momentum \(\mathbf{p}=(p_{1},p_{2},p_{3})\). The equations of motion for a particle in a magnetic field \(\mathbf{B}(\mathbf{x})\) are:

\[\frac{d}{dt}\left(\begin{array}{c}\mathbf{p}\\ \mathbf{x}\end{array}\right)=\left(\begin{array}{cc}-\frac{q}{m}\widehat{B}( \mathbf{x})&-\mathbb{I}_{3}\\ \mathbb{I}_{3}&0\end{array}\right)\left(\begin{array}{c}\frac{\partial H}{ \partial\mathbf{p}}\\ \frac{\partial H}{\partial\mathbf{x}}\end{array}\right):=\mathbb{B}(\mathbf{x })\nabla_{(\mathbf{p},\mathbf{x})}H \tag{48}\]

Notice that equations (48) are not of Lie-Poisson form. However, explicit solutions for the transformations on each time step can be found here as well. We thus believe that this problem is a useful case for demonstrating the power and applicability of these methods beyond the Lie-Poisson equations.

The Hamiltonian for simulations is taken as

\[H(\mathbf{x},\mathbf{p})=\frac{1}{2m}|\mathbf{p}|^{2}+q\varphi(\mathbf{x})\,. \tag{49}\]

Here, \(\mathbb{I}_{3}\) is, as usual, a \(3\times 3\) unity matrix, and we used the hat map notation

\[\widehat{B}(\mathbf{x})=\left(\begin{array}{ccc}0&-B_{3}(\mathbf{x})&B_{2}( \mathbf{x})\\ B_{3}(\mathbf{x})&0&-B_{1}(\mathbf{x})\\ -B_{2}(\mathbf{x})&B_{1}(\mathbf{x})&0\end{array}\right)\,. \tag{50}\]

Note that \(\widehat{B}(\mathbf{x})\mathbf{v}=\mathbf{B}(\mathbf{x})\times\mathbf{v}\) for all \(\mathbf{x}\in\mathbb{R}^{3}\). Similar to (Bajars, 2023; Jin et al., 2022), we take the following values for parameters, electric potential \(\varphi(\mathbf{x})\) and the
magnetic field \(\mathbf{B}(\mathbf{x})\):

\[\begin{split}& q=1,\quad m=1,\\ &\mathbf{B}=\left(0,0,B_{3}\right),\quad\text{with}\quad B_{3}( \mathbf{x})=\sqrt{x_{1}^{2}+x_{2}^{2}}\,,\\ &\varphi(\mathbf{x})=\frac{1}{100\sqrt{x_{1}^{2}+x_{2}^{2}}}\,. \end{split} \tag{51}\]

One can readily check that (48) possesses no Casimirs since \(\mathbb{B}\) is non-degenerate.

_Reduction of motion and conserved quantities._ In (Jin et al., 2022), the initial conditions were chosen to be \((x_{3}=0,v_{3}=0)\) so the particle would always move on the plane. The motion in that case is four-dimensional, and we only need four test Hamiltonians. We thus take the Hamiltonians linear in velocities \((v_{1},v_{2})\) and coordinates \((x_{1},x_{2})\) and compute the corresponding motion.

Note that when the system (48) is restricted so both \(\mathbf{x}\) and \(\mathbf{p}\) are in the plane, _i.e._ both \(x_{3}=0\) and \(p_{3}=0\), and for the choice of any \(\mathbf{B}=\mathbf{e}_{3}B_{3}(r)\) and \(\varphi=\varphi(r)\), where \(r=\sqrt{x_{1}^{2}+x_{2}^{2}}\), there are two integrals of motion. One is clearly the Hamiltonian (49). Another one can be found by considering the evolution for the angular momentum in the \(\mathbf{e}_{3}\) direction \(M_{3}=\mathbf{e}_{3}\cdot(\mathbf{x}\times\mathbf{p})\). We can observe that

\[\dot{M}_{3}=-qB_{3}(r)\mathbf{x}\cdot\dot{\mathbf{x}}=-qB_{3}r\dot{r} \tag{52}\]

leading to the conservation law that for (51) and \(q=1\) reduces to

\[I=M_{3}+q\int^{r}B_{3}(s)s\mathrm{d}s=x_{1}p_{2}-x_{2}p_{1}+\frac{1}{3}\left(x _{1}^{2}+x_{2}^{2}\right)^{3/2}=\text{const.} \tag{53}\]

Therefore, the system (48), for the choice of functions (51) and reduced to four-dimensional motion, is essentially a two-dimensional motion because of the conservation of the Hamiltonian (49) and (53). Thus, one should expect a limited richness of solution behavior. Nevertheless, it is a good test problem and since it has been used in both recent papers on the subject (Bajars, 2023; Jin et al., 2022), we study this particular case as well.

_Lie-Poisson transformations._ Suppose we have a set of \(N\) data pairs, each pair is obtained by the phase flow from \(\mathbf{y}_{i}=(\mathbf{x}_{i},\mathbf{p}_{i})\) to some value \(\mathbf{y}_{i}^{f}=(\mathbf{x}_{i}^{f},\mathbf{p}_{i}^{f})\). If these pairs are obtained from a single trajectory, then \(\mathbf{y}_{i}^{f}=\mathbf{y}_{i+1}\), although this does not have to be the case - our method is capable of learning from several trajectories.
In order to apply our method, we need to compute the results of phase flows for Hamiltonians linear in coordinates \(\mathbf{x}\) and momenta \(\mathbf{p}\).

We just present the answers for these transformations here for brevity; an interested reader may readily check these formulas. In what follows, we shall use the function \(\mathbf{\Phi}\) of time \(t\), initial conditions \(\mathbf{X}=\mathbf{x}(0)\) and parameters \(\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2})\) defined as follows:

\[\mathbf{\Phi}(t;\mathbf{X},\boldsymbol{\alpha})=\int_{0}^{t}B_{3}(\mathbf{X}+ \boldsymbol{\alpha}s)\mathrm{d}s\,. \tag{54}\]

For the expression of \(B_{3}(\mathbf{x})\) taken in (51), both components of the function \(\mathbf{\Phi}\) can be expressed in terms of a single elementary function coming from taking the quadrature of (54):

\[\begin{split}\Phi(t;\mathbf{X},\boldsymbol{\alpha})& =\frac{1}{2A^{3}}\left(A(X_{1}\alpha_{1}+X_{2}\alpha_{2}+A^{2}t) \sqrt{x_{1}(t)^{2}+x_{2}(t)^{2}}\right.\\ &\qquad+\left.(X_{1}\alpha_{2}-X_{2}\alpha_{1})^{2}\right.\\ &\qquad\times\log\left[A\sqrt{x_{1}(t)^{2}+x_{2}(t)^{2}}+(X_{1} \alpha_{1}+X_{2}\alpha_{2}+A^{2}t)\right]\right)\,,\\ A&:=|\boldsymbol{\alpha}|=\sqrt{\alpha_{1}^{2}+ \alpha_{2}^{2}}\,,\quad x_{i}(t):=X_{i}+\alpha_{i}t,\,i=1,2\,.\end{split} \tag{55}\]

The presence of explicit expression for the function (54) is helpful, but it is not essential; for general \(B_{3}(\mathbf{x})\), one can use the quadrature expressions (54). The transformations for each test Hamiltonian are:

1. \(H_{1}=\alpha_{1}p_{1}+\alpha_{2}p_{2}\) leads to the motion \(\mathbf{y}=\mathbf{T}_{1}(t,\mathbf{p}_{0},\mathbf{x}_{0})\) \[\left(\begin{array}{c}p_{1}(t)\\ p_{2}(t)\\ x_{1}(t)\\ x_{2}(t)\end{array}\right)=\mathbf{T}_{1}(t,\mathbf{p}_{0},\mathbf{x}_{0})= \left(\begin{array}{c}-\Phi(t;\mathbf{x}_{0},\boldsymbol{\alpha})+\Phi(0; \mathbf{x}_{0},\boldsymbol{\alpha})+p_{1}(0)\\ \Phi(t;\mathbf{x}_{0},\boldsymbol{\alpha})-\Phi(0;\mathbf{x}_{0},\boldsymbol{ \alpha})+p_{2}(0)\\ \alpha_{1}t+x_{1}(0)\\ \alpha_{2}t+x_{1}(0)\end{array}\right),\] (56)
2. \(H_{2}=\beta_{1}x_{1}+\beta_{2}x_{2}\) leads to the motion \(\mathbf{y}=\mathbf{T}_{2}(t,\mathbf{p}_{0},\mathbf{x}_{0})\) \[\left(\begin{array}{c}p_{1}(t)\\ p_{2}(t)\\ x_{1}(t)\\ x_{2}(t)\end{array}\right)=\mathbf{T}_{2}(t,\mathbf{p}_{0},\mathbf{x}_{0})= \left(\begin{array}{c}-\beta_{1}t+p_{1}(0)\\ -\beta_{2}t+p_{2}(0)\\ x_{1}(0)\\ x_{2}(0)\end{array}\right).\] (57)

Clearly, \(\mathbf{T}_{2}\) does not alter the coordinates \(x_{1,2}\). Thus, the algorithm of computation is very explicit, as \(\alpha_{1,2}\) and \(\beta_{1,2}\) can be solved without any root-finding procedure in the learning stage.
Data preparation.Divide the interval \(\Delta t=h\) into two equal steps \(h/2\). On each of those steps, for each pair of data points \((\mathbf{y}_{i},\mathbf{y}_{i}^{f})\) where \(\mathbf{y}=(p_{1},p_{2},x_{1},x_{2})\), we compute the parameters \((\alpha_{1,i},\alpha_{2,i},\beta_{1,i},\beta_{2,i})\) as follows:

1. Compute \(\alpha_{1,2}\) at the \(i\)-th data point to match \(x_{1,2}\), and corresponding new intermediate momenta \(p_{1,2}^{*}\) by using the corresponding transformations \(\mathbf{T}_{1,2}\) according to \[\alpha_{1,i} =\frac{2}{h}\left(x_{1,i}^{f}-x_{1,i}\right)\] (58) \[\alpha_{2,i} =\frac{2}{h}\left(x_{2,i}^{f}-x_{2,i}\right)\] \[p_{1,i}^{*} =p_{1}-\Phi(h/2;\mathbf{x}_{0},\boldsymbol{\alpha})+\Phi(0; \mathbf{x}_{0},\boldsymbol{\alpha})\] \[p_{2,i}^{*} =p_{2}+\Phi(h/2;\mathbf{x}_{0},\boldsymbol{\alpha})-\Phi(0; \mathbf{x}_{0},\boldsymbol{\alpha})\,.\]
2. Compute \(\beta_{1,2}\) at the \(i\)-th data point according to \[\beta_{1,i} =\frac{2}{h}\left(p_{1,i}^{*}-p_{1,i}^{f}\right)\] (59) \[\beta_{2,i} =\frac{2}{h}\left(p_{2,i}^{*}-p_{2,i}^{f}\right)\,.\]

Learning using Neural Network.The network will learn the mapping between \(\mathbf{y}\) and the parameters \((\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})\) using \(\mathbf{y}_{i}\) as inputs and corresponding parameters \((\boldsymbol{\alpha}_{i},\boldsymbol{\beta}_{i})\) as outputs.

The network structure consists of an input layer taking four inputs: \((\mathbf{x},\mathbf{p})\), three densely connected hidden layers of 32 neurons, with the sigmoid activation function, and the output layer producing an estimate for four outputs \((\boldsymbol{\alpha},\boldsymbol{\beta})\). The number of trainable parameters in the network is 2404. Adam optimization algorithm uses the step \(10^{-3}\) decaying exponentially to \(10^{-5}\), with the number of epochs equal to \(2\cdot 10^{5}\). The loss function is taken to be the mean square average (MSE). Because of the large number of parameters, to avoid overfitting, we take 20000 data points of a single trajectory separated by the time step \(h=0.1\). For that trajectory, use the initial conditions \(\mathbf{v}_{0}=(1,0.5)\) and \(\mathbf{x}_{0}=(0.5,1)\). From all the data points, 80% are used for learning and 20% for evaluation. At the end of the learning procedure, both the loss and the validation loss drop to values below \(10^{-5}\).

Solution evaluations.We choose the initial condition to coincide with the end of the learning trajectory, which is located at \(\mathbf{x}_{0}\simeq(-0.627,0.985)^{T}\) and
\((0.119,1.112)^{T}\). We present the results of simulations in Figure 8. The ground truth was obtained by solving (48) numerically with the given initial conditions using the BDF solver of \(SciPy\), with the relative tolerance of \(10^{-13}\) and absolute tolerance \(10^{-14}\) for \(t=200\), providing outputs every \(h=0.1\). The ground truth is presented with a solid blue line.

The iterative solution using the evaluation provided by successful iterations of Poisson maps (56) and (57) with the parameters \((\alpha,\beta)\) provided by the neural network, approximates the solution at the time points \(t=ih\), \(i=1,\ldots 2000\). The iteration starts with the same initial conditions as the ground truth solutions. The results of the iterations are presented in Figure 8 with solid red lines. As one can see, the LPNets solution approximates ground truth quite well. To further investigate the results, on the left panel Figure 9 we show the mean-square deviation between the ground truth and LPNets solutions of (48) presented in Figure 8. On the right panel of this Figure, we present the conserved quantities: Hamiltonian (Energy) \(H\) given by (49) and the conserved quantity (53). As one can see, the Hamiltonian is preserved with the relative accuracy of \(1-2\%\) and the integral (53)

Figure 8: Solutions of (48): ground truth (blue lines) vs an iterative solution obtained by LPNets (red lines).
to about \(3-4\%\).

Comparison with previous literatureThe problem of a particle in a magnetic field was considered in both (Jin et al., 2022) and (Bajars, 2023). The conservation law (53) has not been identified in these papers, thus we do not provide a direct comparison. The visual agreement of the solutions is just as good as the one demonstrated in (Bajars, 2023; Jin et al., 2022). The absolute error in solution and the Hamiltonian is also similar to that presented in (Bajars, 2023). There are no Casimirs in the system, so there is no additional advantage over the methods presented in (Jin et al., 2022) or (Bajars, 2023).

Additionally, although the system (48) is Poisson, it is not in the Lie-Poisson form. The fact that LPNets can solve this problem on par with other methods is a consequence of the particular form of the magnetic field \(\mathbf{B(x)}\) and the fact that the corresponding integrals (54) can be computed explicitly. Even though the system (48) is not Lie-Poisson, we found it useful to present the results in order to show possible extensions of the methods of LPNets for more general systems.

### Kirchhoff's equations for an underwater vehicle

The motion of a neutrally buoyant underwater body is described by the Kirchhoff equations, see (Leonard, 1997; Leonard and Marsden, 1997) for the discussion on the Lie-Poisson structure of these equations and the corresponding stability results. When the centre of gravity coincides with the centre of buoyancy, the system simplifies somewhat but still possesses a rich dynamics with several integrable and chaotic cases with a rich history of study (Holmes et al., 1998).

Figure 9: Left: the value of Hamiltonian (energy) obtained from (49) and the conserved quantity (53). The ground truth is plotted as blue lines, and the results from the iterative solution obtained by LPNets are shown as red lines. Right: mean-square discrepancy (in all components) between the ground truth and the results obtained by the LPNets.
We will treat that particular system with coincidence of centres of buoyancy and gravity and show the applicability of LPNets approach there as well.

The Hamiltonian for Kirchhoff's equations consists of the kinetic energy of rotational motion and the energy of the translational motion with velocity \(\mathbf{v}\) and mass tensor \(\mathbb{M}\):

\[H(\mathbf{\Pi},\mathbf{p})=\frac{1}{2}\mathbf{\Pi}\cdot\mathbb{T}^{-1}\mathbf{\Pi}+\frac{1} {2}\mathbf{p}\cdot\mathbb{M}^{-1}\mathbf{p}\,. \tag{60}\]

The Poisson bracket for the underwater vehicle is expressed as the Lie-Poisson bracket for the group of rotations and translations \(SE(3)\):

\[\{F,H\}=-\mathbf{\Pi}\cdot\left(\frac{\partial F}{\partial\mathbf{\Pi}}\times\frac{ \partial H}{\partial\mathbf{\Pi}}\right)-\mathbf{p}\cdot\left(\frac{\partial F}{ \partial\mathbf{\Pi}}\times\frac{\partial H}{\partial\mathbf{p}}-\frac{\partial H }{\partial\mathbf{\Pi}}\times\frac{\partial F}{\partial\mathbf{p}}\right)\,. \tag{61}\]

This bracket has a specific form coming from the geometry of semidirect product groups, see (Holm et al., 1998, 2009; Leonard, 1997). Kirchhoff's equations of motion for the underwater vehicle are

\[\dot{\mathbf{H}} =-\frac{\partial H}{\partial\mathbf{\Pi}}\times\mathbf{\Pi}-\frac{ \partial H}{\partial\mathbf{p}}\times\mathbf{p} \tag{62}\] \[\dot{\mathbf{p}} =-\frac{\partial H}{\partial\mathbf{\Pi}}\times\mathbf{p}\,.\]

Equations (62) have two Casimirs: \(C_{1}=\|\mathbf{p}\|^{2}\) and \(C_{2}=\mathbf{\Pi}\cdot\mathbf{p}\). In addition, the total energy given by (60) is also conserved.

_LPNets for Kirchhoff's equations._ Bearing in mind that now we have two momenta \((\mathbf{\Pi},\mathbf{p})\), we shall take the following Hamiltonians: \(H_{1}=\mathbf{A}\cdot\mathbf{\Pi}\) and \(H_{2}=\mathbf{b}\cdot\mathbf{p}\). Equations of motion (62) reduce to

\[\begin{split} H_{1}&=\mathbf{A}\cdot\mathbf{\Pi}:\quad \dot{\mathbf{\Pi}}=-\mathbf{A}\times\mathbf{\Pi}\,,\quad\dot{\mathbf{p}}=-\mathbf{A} \times\mathbf{p}\\ H_{2}&=\mathbf{b}\cdot\mathbf{p}:\quad\dot{\mathbf{\Pi} }=-\mathbf{b}\times\mathbf{p}\,,\quad\dot{\mathbf{p}}=\mathbf{0}\,.\end{split} \tag{63}\]

The first motion is the simultaneous rotation of the vectors \(\mathbf{\Pi}\) and \(\mathbf{p}\) about the same axis \(\mathbf{A}\), by the same amount, with a given angular velocity. This is the transformation \((\mathbf{\Pi},\mathbf{p})\to(\mathbb{R}(\mathbf{A},\theta)\mathbf{\Pi},\mathbb{R}( \mathbf{A},\theta)\mathbf{p})\), where \(\mathbb{R}(\mathbf{A},\theta)\) is the rotation matrix about the axis \(\mathbf{A}\) by the angle \(\theta\). The second motion creates the transformation \((\mathbf{\Pi},\mathbf{p})\to(\mathbf{\Pi}-\mathbf{b}\times\mathbf{p},\mathbf{p})\). These transformation describe the coadjoint action of \(SE(3)\) on an element \(\mathfrak{se}(3)^{*}\), see A.

Based on (63), we thus suggest LPNets for \(SE(3)\)-based Lie-Poisson equations (62).
#### Data preparation.

1. Select the training data couples \((\mathbf{y}^{0},\mathbf{y}^{f})\), where \(\mathbf{y}:=(\boldsymbol{\Pi},\mathbf{p})\).
2. Find the rotation axis \(\mathbf{A}_{j}\) and angles \(\theta_{j}\) that takes \(\mathbf{p}^{0}_{j}\) to \(\mathbf{p}^{f}_{j}\), and the corresponding rotation matrix \(\mathbb{R}(\mathbf{A}_{j},\theta_{j})\). This could be accomplished by finding rotation angles, for example, the Euler or Tait angles, introducing the rotation mapping \(\mathbf{p}^{0}_{j}\) into \(\mathbf{p}^{f}_{j}\). Simultaneous rotation about either \(\mathbf{p}^{0}_{j}\) or \(\mathbf{p}^{f}_{j}\) does not change the end result and must be discarded. While it is possible to proceed in this manner for higher-dimensional groups, for the three-dimensional groups we can utilize a shortcut based on the cross-product of the two vectors, as we have done in the case of a rigid body. We compute \[\mathbf{A}_{j}=\frac{1}{h}\frac{\mathbf{p}^{0}_{j}\times\mathbf{p}^{f}_{j}}{ \|\mathbf{p}^{0}_{j}\|\|\mathbf{p}^{f}_{j}\|}\,.\] (64) The vector \(\mathbf{A}_{j}\) contains all the information necessary for the first step of the algorithm, namely simultaneous rotation. In order to reconstruct the vector \(\mathbf{A}_{j}\) on each time step, we only need the components normal to the vector \(\mathbf{p}^{0}_{j}\). These components can be found by defining two vectors spanning the plane normal to \(\mathbf{p}^{0}_{j}\) in the following way. Take a fixed vector, for example, \(\mathbf{e}_{1}=(1,0,0)^{T}\), and for each \(\mathbf{p}^{0}_{j}\) define \[\boldsymbol{\xi}^{1}_{j}=\frac{\mathbf{p}^{0}_{j}\times\mathbf{e}_{1}}{\| \mathbf{p}^{0}_{j}\times\mathbf{e}_{1}\|}\,,\quad\boldsymbol{\xi}^{2}_{j}= \frac{\mathbf{p}^{f}_{j}\times\boldsymbol{\xi}_{1}}{\|\mathbf{p}^{0}_{j} \times\mathbf{e}_{1}\|}\,.\] (65) Since \((\boldsymbol{\xi}_{1},\boldsymbol{\xi}_{2})\) are unit length and orthogonal to each other and also to \(\mathbf{p}^{0}_{j}\), the vector \(\mathbf{A}_{j}\) defined by (64) can be uniquely reconstructed by \[\mathbf{A}_{j}=\sum_{a=1}^{2}a^{a}_{j}\boldsymbol{\xi}^{a}_{j}\,,\quad a^{a}_ {j}:=\mathbf{A}_{j}\cdot\boldsymbol{\xi}^{a}_{j}\,.\] (66) Variables \(a^{1,2}_{j}\) are the first part of the input for any pair of data points.
3. Find the vector \(\mathbf{b}_{j}\) (or, more precisely, two of its components normal to \(\mathbf{p}_{f}\)), such that \(\boldsymbol{\Pi}^{f}_{j}-(\mathbb{R}(\mathbf{A}_{j},\theta_{j})\boldsymbol{ \Pi}_{0}-\mathbf{b}_{j}\times\mathbf{p}^{f}_{j})\) vanishes. This is accomplished by the following calculation. Define two vectors \(\mathbf{E}_{1,2}\) as \[\mathbf{E}_{1}=\frac{\mathbf{p}^{0}_{j}\times\mathbf{p}^{f}_{j}}{\|\mathbf{p}^ {0}_{j}\times\mathbf{p}^{f}_{j}\|}\,,\quad\mathbf{E}_{2}=\frac{\mathbf{E}_{1} \times\mathbf{p}^{f}_{j}}{\|\mathbf{E}_{1}\times\mathbf{p}^{f}_{j}\|}\,.\] (67)
Clearly, \(\mathbf{E}_{1,2}\) are orthogonal to each other and also to the vector \(\mathbf{p}_{j}^{f}\). Note that \(\mathbf{E}_{1}\) is simply the normalized version of vector \(\mathbf{A}_{j}\) given by (64).
4. We only need to equalize the components of \(\mathbf{\Pi}_{j}^{f}-\mathbf{\Pi}_{j}^{0}\) which are normal to \(\mathbf{p}_{j}^{f}\) on each time step. We thus define the coefficients \(\widetilde{b}^{1,2}\)6 according to Footnote 6: We used tildes to distinguish the coefficients of expansion \(\widetilde{b}^{1,2}\) in the basis of moving vectors \((\mathbf{E}_{1},\mathbf{E}_{2})\) from the coefficients \(b^{1,2}\) in the fixed frame basis \((\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3})\)

\[\widetilde{b}^{1,2}=(\mathbf{\Pi}^{f}-\mathbf{\Pi}^{0})\cdot\mathbf{E}_{1,2} \,,\quad\mathbf{b}=\widetilde{b}^{1}\mathbf{E}_{1}+\widetilde{b}^{2}\mathbf{E }_{2} \tag{68}\]

One can see that by construction, using this algorithm exactly conserves the Casimirs \(|\mathbf{p}|^{2}\) and \(\mathbf{\Pi}\cdot\mathbf{p}\) on every time step.
5. Each pair of mappings \((\mathbf{\Pi}_{j}^{0},\mathbf{p}_{j}^{0})\to(\mathbf{\Pi}_{j}^{f},\mathbf{p}_ {j}^{f})\) (six coordinates) is parameterized by four variables: \((a_{j}^{1,2},\widetilde{b}_{j}^{1,2})\). These four variables are the outputs of the neural networks, whereas \((\mathbf{\Pi}_{j}^{0},\mathbf{p}_{j}^{0})\) are the inputs.

Neural network trainingSimilar to the procedure for the rigid body, we generate 50 trajectories starting in the neighborhood of \(\mathbf{\Pi}^{0}=(1,1,1)\) and \(\mathbf{p}^{0}=(-1,1,2)\). The initial points for trajectories are randomly distributed in the phase space with a uniform distribution in a cube of the size \(2a=0.2\) in every direction. Each trajectory has a 1000 points (not counting the initial point) separated by the time interval of \(h=0.1\), with \(50,000\) data points total used for learning.

EvaluationEvaluation of trajectories using Neural network mapping \((\mathbf{\Pi},\mathbf{p})\) to the variables \((\mathbf{A},\widetilde{b}^{1,2})\) proceeds in next steps:

1. Suppose the initial conditions \((\mathbf{\Pi}^{0},\mathbf{p}^{0})\) are given. Using the Neural Network, estimate the values \((a^{1,2},\widetilde{b}^{1,2})\) corresponding to the inputs.
2. For each starting point of the momentum \(\mathbf{p}^{0}\) and the pair \(a^{1,2}\) evaluated by the neural network, build two vectors \((\boldsymbol{\xi}_{1},\boldsymbol{\xi}_{2})\perp\mathbf{p}^{0}\) according to (65) and reconstruct the axis of rotation \(\mathbf{A}\) according to (64).
3. Compute the three-dimensional rotation matrix \(\mathbb{R}(\mathbf{n}_{\mathbf{A}},\phi)\) about the axis \(\mathbf{n}_{\mathbf{A}}=\mathbf{A}/\|\mathbf{A}\|\) and the angle \(\phi=\operatorname{arc}\,\sin\Big{(}\|\mathbf{A}\|h/\|\mathbf{p}_{0}\|^{2} \Big{)}\).
4. Compute the new linear momentum \(\mathbf{p}^{f}\) and intermediate angular momentum \(\boldsymbol{\Pi}^{f}\) according to \[\mathbf{p}^{f}=\mathbb{R}(\mathbf{n}_{\mathbf{A}},\phi)\mathbf{p}^{0}\,,\quad \boldsymbol{\Pi}^{*}=\mathbb{R}(\mathbf{n}_{\mathbf{A}},\phi)\boldsymbol{\Pi}^{0 }\,.\] (69)
5. Compute the axes \(\mathbf{E}_{1,2}\) according to (67) and update the angular momentum according to \[\boldsymbol{\Pi}_{f}=\boldsymbol{\Pi}_{*}-\widetilde{b}_{1}\mathbf{E}_{1}- \widetilde{b}_{2}\mathbf{E}_{2}\,.\] (70)
6. Reset the final values \((\boldsymbol{\Pi}^{f},\mathbf{p}^{f})\) to the initial conditions and repeat Step 1.

One can see that the algorithm formulated above conserves the Casimirs \(|\mathbf{p}|^{2}\) and \(\boldsymbol{\Pi}\cdot\mathbf{p}\) exactly (_i.e._, to machine precision) on every time step during the trajectory prediction phase.

_Simulation results._

**Data generation** We present the simulation results for the LPNets for the Kirchhoff's equations. The data were obtained using 500 trajectories of 100 time steps with the step \(h=0.1\). All trajectories started in the neighborhood of \(\boldsymbol{\bar{\Pi}}^{0}=(1,1,1)^{T}\) and \(\boldsymbol{\bar{p}}^{0}=(-1,1,2)^{T}\). The initial conditions are uniformly distributed in a cube of the size \(2a=0.2\) in every direction of the phase space. The mass matrix \(\mathbb{M}\) and tensor of inertia \(\mathbb{I}\) in (62) are taken to be \(\mathbb{M}=\text{diag}(1,2,1)\) and \(\mathbb{I}=\text{diag}(1,2,1)\).

**Network parameters** The \(50,000\) data pairs are used to produce \(10,000\) points \((\boldsymbol{\Pi}_{j},\mathbf{p}_{j})\) as input and the parameters \((a^{1},a^{2},\tilde{b}^{1},\tilde{b}^{2})\) as outputs. A fully connected network with 6 inputs, 4 outputs and 6 hidden layers with 64 neurons each, is initialized to describe the mapping. The total number of trainable parameters for the network is 21,508. All activation functions are taken to be sigmoid. The network is trained using Adam algorithm with the step of 0.001, using mean square error as the loss function. From the data, 80% of data are used for training and 20% for validation. The network is optimized using Adams algorithm with the learning rate of \(0^{-3}\) decaying exponentially to \(10^{-5}\) for \(500,000\) epochs. After the optimization procedure, he loss and validation loss reaching the values of \(1.65\cdot 10^{-4}\) and \(1.35\cdot 10^{-3}\), respectively.

**Evaluation of a trajectory** A trajectory is reproduced using the LPNets algorithm as described above, using trained network, starting with the initial conditions \(\boldsymbol{\bar{\Pi}}^{0}=(1,1,1)^{T}\) and \(\boldsymbol{\bar{p}}^{0}=(-1,1,2)^{T}\) until \(t=10\). As usual, the ground truth
solution is produced by a high accuracy BDF algorithm with relative tolerance of \(10^{-13}\) and absolute tolerance of \(10^{-14}\). The comparison between the ground truth and the trajectories obtained by LPNets are presented on Figure 10.

On Figure 11, left panel, we show the conservation of the Hamiltonian (top), the first Casimir \(|\mathbf{p}|^{2}\) (middle) and the second \(\mathbf{\Pi}\cdot\mathbf{p}\) (bottom). Notice that the Casimirs are conserved with higher accuracy than the ground truth, although the ground truth is already conserved to about \(10^{-10}\).

One may wonder whether better results can be obtained by using alternative ways of simulating trajectories or more effective implementations of neural networks. We have to keep in mind that Kirchhoff's equations (62) are chaotic (Holmes et al., 1998; Leonard, 1997). We measure the rate of the divergence of nearby trajectories with the same values of the Casimirs, also known as the first Lyapunov exponent, to be about \(\lambda=0.250\) in units 1/time, starting with the same initial conditions as the simulated trajectory. The _minimum_ growth of errors expected by any numerical scheme is thus proportional to \(e^{\lambda t}\)(Ott, 2002). On

Figure 10: Results of simulation of equations (62) and the corresponding solution of using LPNets procedure. Upper panels: three components of the angular momenta \(\mathbf{\Pi}\); lower panels: three components of linear momenta \(\mathbf{p}\).
the right panel of Figure 11 we present the semilogarithmic plot of the growth of errors, and show that it is corresponds to the growth of errors expected from the chaotic properties of the system.

## 6 G-LPNets: Derivation and test cases

### General theory

Up until now, we have developed a method of computations using transformations that are easily computable, but the coefficients of these transformations must be obtained through the action of a Neural Network. We can extend these results and compute the transformation modules with parameters that are learned through an optimization procedure, producing modules preserving the Poisson structure of the bracket. We proceed as follows.

For Lie-Poisson systems, Section 4 showed that the flow generated by a test Hamiltonian, which is linear in momenta \(H(\mu)=\langle\alpha,\mu\rangle\), can be computed analytically. Instead of taking the test Hamiltonian to be a linear function of momenta, here we take the Hamiltonian to be \(\widetilde{H}=\varphi(H)\), where \(\varphi\) is some nonlinear scalar function of the variable \(H\), which is also a constant of motion7. Instead of equa

Figure 11: Left: Conservation of the Hamiltonian \(H\) (top) and two of the Casimirs \(\left|\mathbf{p}\right|^{2}\) (middle) and \(\mathbf{\Pi}\cdot\mathbf{p}\) (bottom), comparing the results of LPNets (red) and ground truth (blue). The Hamiltonian is conserved to the relative accuracy of less than 0.3%. Notice that LPNets conserves the Casimir exactly (to machine precision) and thus exceeds the ground truth in the conservation of Casimirs. Right: The discrepancy between the results of LPNets and the ground truth. Dashed red line: growth of error expected from the Lyapunov exponent in the semilog scale. The growth of error follows the Lyapunov’s exponent; thus, the neural network is performing as well as could be expected for a chaotic system.
tions (22), we obtain:

Footnote 7: The \(\varphi\)-function \(\varphi(H)\) is defined by the following formula:

\[\dot{\mu}_{a}=\varphi^{\prime}(H)C^{d}_{ab}\alpha^{b}\mu_{d}:=\varphi^{\prime}(H )\mathbb{M}(\alpha)^{d}_{a}\mu_{d}=\varphi^{\prime}(H)\mathbb{N}(\mu)_{ab} \alpha^{b}\,,\quad H=\langle\alpha,\mu_{0}\rangle. \tag{71}\]

Noticing that \(\varphi^{\prime}(H)\) is constant under the flow, generated by the Hamiltonian \(\widetilde{H}=\varphi(H)\), we conclude that the transformation generated by the flow of \(\varphi(H)\) is simply obtained by a time scaling of (23) and can be written as

\[\widetilde{\mathbb{T}}(\alpha,t)\mu_{0}=\mathbb{T}(\alpha,\varphi^{\prime}(H)t )\mu_{0}=e^{\mathbb{M}(\alpha)\varphi^{\prime}(H)t}\mu_{0}\,,\quad H=\langle \alpha,\mu_{0}\rangle\, \tag{72}\]

where \(\mathbb{T}(\alpha,t)\) are defined as in (23). The prefactor \(\varphi^{\prime}(H)=\varphi^{\prime}(\langle\alpha,\mu_{0}\rangle)\) can be sought as a function of the variables using some kind of scalar activation function. Thus, in the G-LPNets framework, we are looking for transformations that are compositions of moduli \(\mathbb{T}_{s}(a_{s},\alpha_{s},b_{s},t)\), \(s=1,\ldots M\) having the functional form

\[\widetilde{\mathbb{T}}_{s}(a_{s},\alpha_{s},b_{s},t)\mu_{0}=e^{\mathbb{M}( \alpha_{s})q_{s}t}\mu_{0}\,,\quad q_{s}=a_{s}\sigma\left(\langle\alpha_{s},\mu _{0}\rangle\right)+b_{s}\,, \tag{73}\]

where \(\sigma\) is the activation function, which we will choose to be the sigmoid, and \((a_{s},\alpha_{s},b_{s})\) are the parameters8to be determined. The transformations (72) are generalized moduli preserving the Poisson bracket and the Casimirs to machine precision.

Footnote 8: In the expression (72), \(\alpha_{s}\) refers to the \(s\) vector in the collection of \(M\) vectors \((\alpha_{1},\ldots,\alpha_{M})\) in the Lie algebra \(\mathfrak{g}\). In the example of the applications to the rigid body motion, for simplicity, we will put \(\alpha_{s}\) to be proportional to the basis vectors of the Lie algebra with some scalar coefficients, which we will also call \(\alpha_{s}\).

The method of G-LPNets selects the order of transformations \(\widetilde{\mathbb{T}}_{s}\) and minimizes a mean-square loss function. For example, for the data comprising \(N\) points specifying the beginning \(\mu_{i}^{0}\) and end \(\mu_{i}^{f}\) values of momenta on each time interval, and assuming that the composition takes \(M\) exactly equal time-substeps, the loss function can be taken to be the mean square error (MSE)

\[L=\frac{1}{N}\sum_{i=1}^{N}\left\|\widetilde{\mathbb{T}}_{M}(a_{M},\alpha_{M}, b_{M},h/M)\circ\ldots\circ\widetilde{\mathbb{T}}_{M}(a_{1},\alpha_{1},b_{1},h/M) \mu_{i}^{0}-\mu_{i}^{f}\right\|^{2}. \tag{74}\]

The use of the mean square loss is advantageous since the derivatives of \(L\) with respect to parameters \(a_{s},\alpha_{s},b_{s}\) can be found analytically. Denoting the collection of parameters with a bar for shortness, for example, \(\bar{a}=(a_{1},a_{2},\ldots)\)_etc_, we naturally have

\[(\bar{a},\bar{\alpha},\bar{b})=\text{arg min}\ L(\bar{a},\bar{\alpha},\bar{b}). \tag{75}\]
There are several points to keep in mind regarding the application of G-LPNets.

1. If there are no Casimirs, it is natural to choose the transformations alternating in a given sequence. For example, for an \(n\)-dimensional system one could choose a repeated application of \(\widetilde{\mathbb{T}}_{1}\), \(\widetilde{\mathbb{T}}_{2}\), \(\widetilde{\mathbb{T}}_{M}\) followed by \(\widetilde{\mathbb{T}}_{1}\) etc. It would be natural (although strictly speaking not necessary) to take the depth of network to be \(M=n\cdot k\), where \(k\) is an integer number.
2. If there are \(d\) independent Casimirs \(C_{j}(\mu)\), \(j=1,\ldots,d\), the transformations producing the motion about the axis parallel to \(\nabla C_{j}\) produce no effective evolution of momenta. One can either choose \(n-d\) transformations on every time step by excluding certain \(\alpha\), or simply apply all transformations in a sequence with the understanding that the resulting \(\alpha\) is not unique. The latter is the approach we will use to describe rigid body rotation below.
3. The advantage of applying all \(\widetilde{\mathbb{T}}_{s}\) in a row is that there will be no accuracy loss if, at some point, the momentum is becoming close to parallel with the given coordinate axis. The disadvantage of applying all transformations in a sequence without excluding any \(\alpha\) lies in the necessity for post-processing for (75), as they are only defined up to corresponding components of gradients that are not in the span of the gradients of Casimir functions \(\nabla_{\mu}C_{j}\), \(j=1,\ldots,d\). In fact, even at the continuous level, the Hamiltonian is defined only up to the Casimirs, since the same dynamics are obtained if an arbitrary Casimir is added to the Hamiltonian. Similarly, all possible solutions of (75) obtained by G-LPNets are in the same equivalence class since they generate the same dynamics in phase space.
4. It is also possible that because parameters are defined only up to a certain vector or vectors, one could encounter vanishing gradients in some directions and a slow down of the convergence to the desired solution due to some numerical artifacts. We have not encountered these artifacts during our solution of the rigid body equations, but we cannot exclude further numerical difficulties when applying G-LPNets for general high-dimensional problems.
5. We expect that G-LPNets will also work for the cases beyond the Lie-Poisson framework, whenever explicit integration of the trajectories with the test Hamiltonians is possible, such as the particle in a magnetic field. The transformations \(\widetilde{\mathbb{T}}_{s}\) are then computed as the generalizations of the Poisson transformations with the unknown time scaling coefficients. As long as
the completeness of these transformations will be achieved, we expect them to provide efficient and accurate data-based computing methods for general Poisson problems.

The most challenging part of applying G-LPNets, in our opinion, is the lack of a general completeness result. It would be nice if the G-LPNets moduli (73) satisfied a completeness result analogous to that of SympNets (Jin et al., 2020). Right now, it seems unlikely that a general result may be proved valid for all Lie-Poisson brackets. Without specifying more information about the particular Lie group, progress in that area may be limited. However, the silver lining here is that for each particular problem, the symmetry of the problem is bound to be known _a priori_, where the exact value of the Hamiltonian may or may not be known. Thus, if one focuses on particular problem at hand, a completeness result of transformations leading to G-LPNets could be feasible to achieve. It is also possible that some of the methods of analysis performed in Jin et al. (2020) for proving completeness of the modules in symplectic space would be applicable in the more general setting of a particular Lie-Poisson bracket.

### Applications to rigid body dynamics

To show the potential power of G-LPNets, we treat the equations of a rigid body, following up on our discussion in Sec. 5.1. Now, instead of learning from trajectories that stay close to the desired area, we aim to learn the whole dynamics and simply choose 50 initial points uniformly distributed in a cube in momentum space \(\mathbf{\Pi}\), \(-2\leq\Pi_{a}\leq 2\), \(a=1,2,3\). Each trajectory is simulated with a high precision ODE solver and output is provided every \(h=0.1\), providing 20 data pairs, with the total of 1000 data pairs. All parameters of the system are exactly as in Sec. 5.1.

Our comparison is the reconstruction of the dynamics of a rigid body which was already done in (Bajars, 2023). Note, however, that (Bajars, 2023) takes all data pairs on _the same_ Casimir surface. In our opinion, since the Casimir surface depends on the initial conditions, a physical system, such as a satellite, could be observed on different Casimir surfaces due to the fact that the thrusters or other external forces have moved it from one Casimir surface to another outside the observation time. The Hamiltonian and physical parameters of the satellite are assumed to be the same, so it makes sense that the ground truth is generated with the same Hamiltonian, but different Casimirs.

After data points are generated, a G-LPnet with 6 transformations (18 parameters) is generated. The transformations \(\widetilde{\mathbb{T}}_{i}\) are rotations about the coordinate axes
in the _fixed_ frame \(\mathbf{e}_{1,2,3}\) by the angle \(\phi_{i}\) (_i.e._, \(e_{1}=(1,0,0)^{T}\)_etc._) The rotations proceed in the sequence

\[\widetilde{\mathbb{T}}_{1}\rightarrow\widetilde{\mathbb{T}}_{2}\rightarrow \widetilde{\mathbb{T}}_{3}\rightarrow\widetilde{\mathbb{T}}_{1}\rightarrow\ldots \tag{76}\]

At the given step \(s\), the rotation angle \(\phi_{s}\) for the given value of the momentum \(\mathbf{\Pi}^{0}\) is given by

\[\phi_{s}=a_{s}\sigma(\mathbf{A}_{s}\cdot\mathbf{\Pi}^{0})+b_{s}\,. \tag{77}\]

The optimization finds the values of \((\bar{\alpha},\bar{\alpha},\bar{b})\) minimizing the MSE loss function. The gradients of MSE functions with respect to parameters are computed analytically. We tried gradient descent methods and discovered that while they work satisfactorily, they do require quite a substantial number of epochs to converge, as was already observed in (Bajars, 2023; Jin et al., 2022). To make computation more efficient, we implemented an optimization procedure based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Fletcher, 2000) into SciPy: www.scipy.org. We let the BFGS algorithm run for 3500 iterations achieving the value of the Loss of less than \(5\cdot 10^{-10}\), a procedure which is several orders of magnitude more efficient than the standard gradient descent-based algorithm.

Using the values of parameters found by the optimization procedure, 10 long solution with \(10,000\) data points each with the time step of \(h=0.1\) (max time \(t_{max}=1000\)) were generated and compared with the ground truth solution obtained by a high precision ODE method. These long-term solutions were generated to compare with the results in (Bajars, 2023). The initial conditions \(\mathbf{\Pi}_{0}\) for these solutions were taken to be random numbers chosen from a random distribution from a cube \(|\Pi_{a}|\leq 1.25\), \(a=1,2,3\). The value of the cube was somewhat smaller than the initial conditions so the solutions are guaranteed to remain within the area where the data were available.

In the left panel of Figure 12 we show a part of one of the 10 trajectories extending to only about \(t=500\) for clarity. This particular trajectory starts with the initial condition \(\mathbf{\Pi}_{0}=(1.011,1.178,0.585)\), the results for other trajectories are similar. The values of momenta are nearly indistinguishable from ground truth for very long time. In the right panel of this Figure, we show all 10 trajectories in the phase space, plotted for all times \(0\leq t\leq 1000\) (10000 iterations). The results from the simulations and the ground truth coincide perfectly. This is due to the fact that both the Casimirs and the Hamiltonians are preserved with very high accuracy, as Figure 13 shows. The Casimir is, again, conserved to machine precision on each step and thus substantially exceeds the accuracy of the calculation for the ground truth solution (still very high at about \(10^{-9}\)). The relative error in energy is of the order of \(0.1\%\) over all long-term solution.
Comparison with previous literature.The case of learning the whole dynamics of a rigid body was considered in (Bajars, 2023). In that paper, the dynamics was considered only on a _single_ Casimir's surface \(|\mathbf{\Pi}|=1\) with 300 data points. The solutions generated by the neural network were also taken on that Casimir's surface. In contrast, we presented the dynamics with initial data taken from a volume of \(\mathbf{\Pi}\) (a cube), and initial conditions for G-LPNets are also taken without any restriction of the Casimir surface. Our relative error of solutions are of the same order as the results presented in (Bajars, 2023). The conservation of Casimir and energy is reported to be of the order of 1.5% in (Bajars, 2023). In our case, the relative error in energy is somewhat better, with the max value of error around 0.4%. The Casimir \(|\mathbf{\Pi}|^{2}\) is conserved with the machine precision on every time step, so after \(10^{5}\) steps a typical error in Casimir is expected to be of the order of \(10^{-11}\div 10^{-10}\). While the relative error of energy in methods presented in (Bajars, 2023) can no doubt be improved to reach the accuracy achieved here, we are not aware of any method preserving the value of the Casimir to the same precision as ours.

## 7 Conclusions

We have derived a novel method of learning the evolution of a general Lie-Poisson system that is capable of predicting the dynamics in phase space with

Figure 12: Left: Results of G-LPNets applied to the motion of a rigid body (red) versus Ground truth (blue) for the individual momenta. Right: Parametric plot of the momenta in the phase space for 10 solutions, chosen to start from uniformly distributed random points in the phase space in the cube \(|\Pi_{a}|\leq 1.25\) in the phase space. The results of G-LPNets are indistinguishable from the ground truth solution.
high precision. Our method learns the system by applying exact Poisson maps for certain test Hamiltonians, which are exactly solvable for any Lie-Poisson bracket. The resulting maps preserve the Poisson bracket under the evolution and also preserve all Casimirs with machine precision. These methods are also applicable to systems beyond Lie-Poisson, such as a particle in magnetic field, as long as the corresponding equations for test Hamiltonians are exactly solvable.

We derive two types of networks. The first one is the Local Lie-Poisson Neural Networks (LPNets) which derive the local Poisson maps using test Hamiltonians that are linear in momenta. The parameters for these Poisson mappings resulting from the test Hamiltonians are then learned using standard methods of data-based learning with Artificial Neural Network (ANN). The advantage of this method is that the local completeness of the mappings is achieved automatically, since they represent exponential maps on a Lie group. An additional advantage of the method is the ability to use all the modern technology of ANNs for the discovery of the mapping from momenta to the parameters of test Hamiltonians.

An alternative method was derived, called Global LPNets, using an example of nonlinear test Hamiltonians. These nonlinear Hamiltonians are arbitrary functions of the local Hamiltonians in the local LPNets approach. The explicit evolution maps on every time step obtained by these methods can be viewed as

Figure 13: Left panel, top: Relative accuracy for the conservation of the Hamiltonian \(H\) (top) (_i.e._, \(\Delta E/E=(E(t)-E(0))/\langle E\rangle\), where \(\langle E\rangle\) is the mean value of the Hamiltonian (top). Left panel, bottom: the corresponding conservation of the Casimir \(C\), computed as \(\Delta C/C=(C(t)-C(0))/\langle C\rangle\) comparing the results of G-LPNets (red) and ground truth (blue), for all 10 simulations. As usual, G-LPNets conserve the Casimir exactly (to machine precision) and thus substantially exceed the ground truth in the conservation of Casimirs. Right: The discrepancy between the results of G-LPNets and the ground truth for the solution, presented on the left part of Figure 12. Again, the discrepancy comes mostly from time mismatch, whereas the amplitude of oscillations is conserved with high precision due to the corresponding high precision in the conservation of energy.
generalizations of the symplectic modules derived in (Jin et al., 2020) for case of a general Lie-Poisson bracket. We have presented an application of these methods to rigid body motion and showed that these methods demonstrate excellent accuracy and efficiency for long-term computation, and the ability to learn the dynamics in the whole phase space from quite a limited number of observations.

While the G-LPNets seem more computationally effective than the local LPNets, we must caution the reader that a completeness result for mappings for G-LPNets obtained through the test nonlinear Hamiltonians is still missing. We believe it is probably highly unlikely that such a result can exist for a general Lie-Poisson system. The completeness result is likely to depend on the structure of the actual Lie group and the corresponding Lie-Poisson bracket. This is an interesting topic that we intend to address in the future.

Another interesting topic is the question of the extension of this method to more general systems. In order to achieve that goal, the equations generated by the Poisson bracket for the test Hamiltonians must be exactly solvable, as it was in the case of the particle in the magnetic field. It will be interesting to compute the conditions on the Poisson bracket for such integrability to occur, which we also plan to undertake in the future. Of particular interest are constrained systems, especially systems with nonholonomic constraints. The Lie-Poisson equations in this case become the Lie-Poisson-d'Alembert's equation, where the right-hand side of (A.3) contains extra terms enforcing vanishing of momenta projections to certain subspaces defined by the constraints, such as in the Suslov problem (Bloch, 2003). Apart from their importance in classical mechanics, these methods were recently found to be important for the variational discretizations of fluids (Gawlik and Gay-Balmaz, 2020), possibly including irreversible thermodynamics processes (Gawlik and Gay-Balmaz, 2022). Extension of data-based computing to nonholonomic Lie-Poisson systems may thus play an important role in the applications of the method of this paper to continuum mechanics, including irreversible processes.

## 8 Acknowledgements

We are grateful to Anthony Bloch, Pavel Bochev, Stephen Bond, Anthony Gruber, Melvin Leok, Tomoki Ohsawa, Tanya Schmah, Andrew Sinclair, Nathaniel Trask and Dmitry Zenkov for fruitful and engaging discussions. SH acknowledges support and experience provided by the internship in ATCO's transformation team and productive exchange with the team members. SH and VP were partially supported by the NSERC Discovery grant.
This article has been co-authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns right, title and interest in and to the article and is responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan [https://www.energy.gov/downloads/doe-public-access-plan](https://www.energy.gov/downloads/doe-public-access-plan).
