NONSEPARABLE SYMPLECTIC NEURAL NETWORKS
Year: 2021

Abstract:
Predicting the behaviors of Hamiltonian systems has been drawing increasing attention in scientific machine learning. However, the vast majority of the literature was focused on predicting separable Hamiltonian systems with their kinematic and potential energy terms being explicitly decoupled while building data-driven paradigms to predict nonseparable Hamiltonian systems that are ubiquitous in fluid dynamics and quantum mechanics were rarely explored. The main computational challenge lies in the effective embedding of symplectic priors to describe the inherently coupled evolution of position and momentum, which typically exhibits intricate dynamics. To solve the problem, we propose a novel neural network architecture, Nonseparable Symplectic Neural Networks (NSSNNs), to uncover and embed the symplectic structure of a nonseparable Hamiltonian system from limited observation data. The enabling mechanics of our approach is an augmented symplectic time integrator to decouple the position and momentum energy terms and facilitate their evolution. We demonstrated the efficacy and versatility of our method by predicting a wide range of Hamiltonian systems, both separable and nonseparable, including chaotic vortical flows. We showed the unique computational merits of our approach to yield long-term, accurate, and robust predictions for large-scale Hamiltonian systems by rigorously enforcing symplectomorphism.

heading: 1 INTRODUCTION
text: A Hamiltonian dynamic system refers to a formalism for modeling a physical system exhibiting some specific form of energy conservation during its temporal evolution. A typical example is a pendulum whose total energy (referred to as the system’s Hamiltonian) is conserved as a temporally invariant sum of its kinematic energy and potential energy. Mathematically, such energy conservation indicates a specific geometric structure underpinning its time integration, named as a symplectic structure, which further spawns a wide range of numerical time integrators to model Hamiltonian systems. These symplectic time integrators have proven their effectiveness in simulating a variety of energy-conserving dynamics when Hamiltonian expressions are known as a prior. Examples encompass applications in plasma physics (Morrison, 2005), electromagnetics (Li et al., 2019), fluid mechanics (Salmon, 1988), and celestial mechanics (Saari & Xia, 1996), to name a few.
On another front, the emergence of the various machine learning paradigms with their particular focus on uncovering the hidden invariant quantities and their evolutionary structures enable a faithful prediction of Hamiltonian dynamics without knowing its analytical energy expression beforehand. The key mechanics underpinning these learning models lie in a proper embedding of the strong mathematical inductive priors to ensure Hamiltonian conservation in a neural network data flow. Typically, such priors are realized in a variational way or a structured way. For example, in Greydanus et al. (2019), the Hamiltonian conservation is encoded in the loss function. This category of methods does not assume any combinatorial pattern of the energy term and therefore relies on the inherent expressiveness of neural networks to distill the Hamiltonian structure from abundant training datasets (Choudhary et al., 2019). Another category of Hamiltonian networks, which we refer to as structured approaches, implements the conservation law indirectly by embedding a symplectic time integrator (DiPietro et al., 2020; Tong et al., 2020; Chen et al., 2020) or composition of linear, activation, and gradient modules (Jin et al., 2020) into the network architecture.
∗shiying.xiong@dartmouth.edu
.
One of the main limitations of the current structured methods lies in the separable assumption of the Hamiltonian expression. Examples of separable Hamiltonian systems include the pendulum, the Lotka–Volterra (Zhu et al., 2016), the Kepler (Antohe & Gladwell, 2004), and the Hénon–Heiles systems (Zotos, 2015). However, beyond this scope, there exist various nonseparable systems whose Hamiltonian has no explicit expression to decouple the position and momentum energies. Examples include incompressible flows (Suzuki et al., 2007), quantum systems (Bonnabel et al., 2009), rigid body dynamics (Chadaj et al., 2017), charged particle dynamics (Zhang et al., 2016), and nonlinear Schrödinger equation (Brugnano et al., 2018). This nonseparability typically causes chaos and instability, which further complicates the systems’ dynamics. Although SympNet in Jin et al. (2020) can be used to learn and predict nonseparable Hamiltonian systems, multiple matrices of the same order with system dimension are needed in the training process of SympNet, resulting in difficulties in generalizing into high-dimensional large-scale N-body problems which are common in a series of nonseparable Hamiltonian systems, such as quantum multibody problems and vortexparticle dynamics problems. Such chaotic and large-scale nature jointly adds shear difficulties for a conventional machine learning model to deliver faithful predictions.
In this paper, we propose an effective machine learning paradigm to predict nonseparable Hamiltonian systems. We build a novel neural network architecture, named nonseparable symplectic neural networks (NSSNNs), to enable accurate and robust predictions of long-term Hamiltonian dynamics based on short-term observation data. Our proposed method belongs to the category of structured network architectures: it intrinsically embeds the symplectomorphism into the network design to strictly preserve the symplectic evolution and further conserves the unknown, nonseparable Hamiltonian energy. The enabling techniques we adopted in our learning framework consist of an augmented symplectic time integrator to asymptotically “decouple” the position and momentum quantities that were nonseparable in their original form. We also introduce the Lagrangian multiplier in the augmented phase space to improve the system’s numerical stability. Our network design is motivated by ideas originated from physics (Tao, 2016) and optimization (Boyd et al., 2004). The combination of these mathematical observations and numerical paradigms enables a novel neural network architecture that can drastically enhance both the scale and scope of the current predictions.
We show a motivational example in Figure 1 by comparing our approach with a traditional HNN method (Greydanus et al., 2019) regarding their structural designs and predicting abilities. We refer the readers to Section 6 for a detailed discussion. As shown in Figure 1, the vortices evolved using NSSNN are separated nicely as the ground truth, while the vortices merge together using HNN due to the failure of conserving the symplectic structure of a nonseparable system. The conservative capability of NSSNN springs from our design of the auxiliary variables (red x and y) which converts
the original nonseparable system into a higher dimensional quasi-separable system where we can adopt a symplectic integrator.

heading: 2 RELATED WORKS
text: Data-driven physical prediction. Data-driven approaches have been widely applied in physical systems including fluid mechanics (Brunton et al., 2020), wave physics (Hughes et al., 2019), quantum physics (Sellier et al., 2019), thermodynamics (Hernandez et al., 2020), and material science (Teicherta et al., 2019). Among these different physical systems, data-driven fluid receives increasing attention. We refer the readers to Brunton et al. (2020) for a thorough survey of the fundamental machine learning methodologies as well as their uses for understanding, modeling, optimizing, and controlling fluid flows in experiments and simulations based on training data. One of the motivations of our work is to design a versatile learning approach that can predict complex fluid motions. On another front, many pieces of research focus on incorporating physical priors into the learning framework, e.g., by enforcing incompressibility (Mohan et al., 2020), the Galilean invariance (Ling et al., 2016), quasistatic equilibrium (Geng et al., 2020), the Lagrangian invariance (Cranmer et al., 2020), and Hamiltonian conservation (Hernandez et al., 2020; Greydanus et al., 2019; Jin et al., 2020; Zhong et al., 2020). Here, inspired by the idea of embedding physics priors into neural networks, we aim to accelerate the learning process and improve the accuracy of our model.
Neural networks for Hamiltonian systems. Greydanus et al. (2019) introduced Hamiltonian neural networks (HNNs) to conserve the Hamiltonian energy of the system by reformulating the loss function. Inspired by HNN, a series of methods intrinsically embedding a symplectic integrator into the recurrent neural network was proposed, such as SRNN (Chen et al., 2020), TaylorNet (Tong et al., 2020) and SSINN (DiPietro et al., 2020), to solve separable Hamiltonian systems. Combined with graph networks (Sanchez-Gonzalez et al., 2019; Battaglia et al., 2016), these methods were further generalized to large-scale N-body problems induced by interaction force between the particle pairs. Jin et al. (2020) proposed SympNet by directly constructing the symplectic mapping of system variables within neighboring time steps to handle both separable and nonseparable Hamiltonian systems. However, the scale of parameters in SympNet for training N dimensional Hamiltonian system is O(N2), which makes it hard to be generalized to the high dimensional N-body problems. Our NSSNN overcomes these limitations by devising a new Hamiltonian network architecture that is specifically suited for nonseparable systems (see details in Section 5). In addition, the Hamiltonianbased neural networks can be extended to further applications. Toth et al. (2020) developed the Hamiltonian Generative Network (HGN) to learn Hamiltonian dynamics from high-dimensional observations (such as images). Moreover, Zhong et al. (2020) introduced Symplectic ODE-Net (SymODEN), which adds an external control term to the standard Hamiltonian dynamics.

heading: 3 FRAMEWORK

heading: 3.1 AUGMENTED HAMILTONIAN EQUATION
text: We start by considering a Hamiltonian system with N pairs of canonical coordinates (i.e. N generalized positions and N generalized momentum). The time evolution of canonical coordinates is governed by the symplectic gradient of the Hamiltonian (Hand & Finch, 2008). Specifically, the time evolution of the system is governed by Hamilton’s equations as
dq dt = ∂H ∂p , dp dt = −∂H ∂q , (1)
with the initial condition (q,p)|t=t0 = (q0,p0). In a general setting, q = (q1, q2, · · · , qN ) represents the positions and p = (p1, p2, ...pN ) denotes their momentum. Function H = H(q,p) is the Hamiltonian, which corresponds to the total energy of the system. An important feature of Hamilton’s equations is its symplectomorphism (see Appendix B for a detailed overview).
The symplectic structure underpinning our proposed network architecture draws inspirations from the original research of Tao (2016) in computational physics. In Tao (2016), a generic, high-order, explicit and symplectic time integrator was proposed to solve (1) of an arbitrary separable and nonseparable
HamiltonianH. This is implemented by considering an augmented Hamiltonian H(q,p,x,y) := HA +HB + ωHC (2)
with HA = H(q,y), HB = H(x,p), HC = 1
2
( ‖q − x‖22 + ‖p− y‖22 ) (3)
in an extended phase space with symplectic two form dq ∧ dp+ dx ∧ dy, where ω is a constant that controls the binding of the original system and the artificial restraint. Notice that the Hamilton’s equations forH dq dt = ∂H ∂p = ∂H(x,p) ∂p + ω(p− y), dp dt = −∂H ∂q = −∂H(q,y) ∂q − ω(q − x), dx dt = ∂H ∂y = ∂H(q,y) ∂y − ω(p− y),
dy dt = −∂H ∂x = −∂H(x,p) ∂x + ω(q − x),
(4)
with the initial condition (q,p,x,y)|t=t0 = (q0,p0, q0,p0) have the same exact solution as (1) in the sense that (q,p,x,y) = (q,p, q,p). Hence, we can get the solution of (1) by solving (4). Furthermore, it is possible to construct high-order symplectic integrators forH in (4) with explicit updates. Our model aims to learn the dynamical evolution of (q,p) in (1) by embedding (4) into the framework of NeuralODE (Chen et al., 2018). The coefficient ω acts as a regularizer, which stabilizes the numerical results (see Section 4).

heading: 3.2 NONSEPARABLE HAMILTONIAN NEURAL NETWORK
text: We learn the nonseparable Hamiltonian dynamics (1) by constructing an augmented system (4), from which we can obtain the energy function H(q,p) by training the neural network Hθ(q,p) with parameter θ and calculate the gradient∇Hθ(q,p) by taking the in-graph gradient. For the constructed network Hθ(q,p), we integrate (4) by using the second-order symplectic integrator (Tao, 2016). Specifically, we will have an input layer (q,p,x,y) = (q0,p0, q0,p0) at t = t0 and an output layer (q,p,x,y) = (qn,pn,xn,yn) at t = t0 +ndt. Algorithm 1 Integrate (4) by using the secondorder symplectic integrator Input: q0,p0, t0, t, dt; φδ1, φδ2, and φδ3 in (5); Output: (q̂, p̂, x̂, ŷ) = (qn,pn,xn,yn) 1 (q0,p0,x0,y0) = (q0,p0, q0,p0) n = floor[(t− t0)/dt] for i = 1→ n do 2 (qi,pi,xi,yi) = φ dt/2 1 ◦φ dt/2 2 ◦φdt3 ◦φ dt/2 2 ◦ φ dt/2 1 ◦ (qi−1,pi−1,xi−1,yi−1); 3 end The recursive relations of (qi,pi,xi,yi), i = 1, 2, · · · , n, can be expressed by the algorithm 1 (also see Figure 8 in Appendix A). The input functions φδ1(q,p,x,y), φ δ 2(q,p,x,y), and φ δ 3(q,p,x,y) in algorithm 1 are qp− δ[∂Hθ(q,y)/∂q]x+ δ[∂Hθ(q,y)/∂p] y  , q + δ[∂Hθ(x,p)/∂p]px y − δ[∂Hθ(x,p)/∂q]  , and 1 2  ( q + x p+ y ) +Rδ ( q − x p− y ) ( q + x p+ y ) −Rδ ( q − x p− y )  , (5)
respectively. Here
Rδ := [ cos(2ωδ)I sin(2ωδ)I − sin(2ωδ)I cos(2ωδ)I ] , where I is a identity matrix. (6)
We remark that x and y are just auxiliary variables, which are theoretically equal to q and p. Therefore, we can use the data set of (q,p) to construct the data set containing variables (q,p,x,y). In addition, by constructing the networkHθ, we show that theorem B.1 in Appendix B holds, so the networks φδ1,φ δ 2, and φ δ 3 in (5) preserve the symplectic structure of the system. Suppose that Φ1 and Φ2 are two symplectomorphisms. Then, it is easy to show that their composite map Φ2 ◦ Φ1 is also symplectomorphism due to the chain rule. Thus, the symplectomorphism of algorithm 1 can be guaranteed by the theorems B.1.

heading: 4 TRAINING SETTINGS AND ABLATION TESTS
text: We use 6 linear layers with hidden size 64 to model Hθ, all of which are followed by a Sigmoid activation function except the last one. The derivatives ∂Hθ/∂p, ∂Hθ/∂q, ∂Hθ/∂x, ∂Hθ/∂y are all obtained by automatic differentiation in Pytorch (Paszke et al., 2019). The weights of the linear layers are initialized by Xavier initializaiton (Glorot & Bengio, 2010).
We generate the dataset for training and validation using high-precision numerical solver (Tao, 2016), where the ratio of training and validation datasets is 9 : 1. We set the dataset (qj0,p j 0) as the start input and (qj ,pj) as the target with j = 1, 2, · · · , Ns, and the time span between (qj0,p j 0) and (q
j ,pj) is Ttrain. Feeding (q0,p0) = (q j 0,p j 0), t0 = 0, t = Ttrain, and time step dt in Algorithm 1 to get the predicted variables (q̂j , p̂j , x̂j , ŷj). Accordingly, the loss function is defined as
LNSSNN = 1
Nb Nb∑ j=1 ‖q(j) − q̂(j)‖1 + ‖p(j) − p̂(j)‖1 + ‖q(j) − x̂(j)‖1 + ‖p(j) − ŷ(j)‖1, (7)
where Nb = 512 is the batch size of the training samples. We use the Adam optimizer (Kingma & Ba, 2015) with learning rate 0.05. The learning rate is multiplied by 0.8 for every 10 epoches.
Taking systemH(q, p) = 0.5(q2 + 1)(p2 + 1) as an example, we carry out a series of ablation tests based on our constructed networks. Normally, we set the time span, time step and dateset size as T = 0.01, dt = 0.01 and Ns = 1280.
The choice of ω in (4) is largely flexible since NSSNN is not sensitive to the parameter ω when it is larger than a certain threshold. Figure 2 shows the training and validation losses with different ω in the network trained by clean and noise datasets. Though the convergence rates are slightly different in a small scope, the examples with various ω are able to converge to the same size of training and validation losses. Here, we set ω = 2000, but ω can be smaller than 2000. The only requirement for picking ω is that it has to be larger than O(10), which is detailed in Appendix C.
We pick the L1 loss function to train our network due to its better performance. Figure 3 compares the validation losses with different training loss functions in the network trained by clean and noise datasets. Figure 3(a) shows that either the network trained by L1 or MSE with a clean dataset can converge to a small validation loss, but the network trained by L1 loss converges relatively faster. Figures 3(b) and 3(c) both show that the network trained by L1 with noise dataset can converge to a smaller validation loss. In addition, we already introduced a regularization term in the symplectic integrator embedded in the network; thus, there is no need to add the regularization term in the loss function.
The integral time step in the sympletic integrator is a vital parameter, and the choice of dt largely depends on the time span Ttrain. Figure 4 compares the validation losses generated by various integral time steps dt based on fixed dataset time spans Ttrain = 0.01, 0.1 and 0.2 respectively in the training process. The validation loss converges to a similar degree with various dt based on fixed Ttrain = 0.01 and Ttrain = 0.1 in 4(a) and (b), while it increases significantly as dt increases based
on fixed Ttrain = 0.02 in 4(c). Thus, we should take relatively small dt for the dataset with larger time span Ttrain

heading: 5 COMPARISONS WITH OTHER METHODS

heading: 5.1 METHODOLOGIES
text: We compare our method with other recently proposed methods, such as HNN (Greydanus et al., 2019), NeuralODE (Chen et al., 2018), TaylorNet (Tong et al., 2020), SSINN (DiPietro et al., 2020), SRNN (Chen et al., 2020), and SympNet (Jin et al., 2020). There are several features distinguishing our method from others, as shown in Table 1. HNN first enforces conservative features of a Hamiltonian system by reformulating its loss function, which incurs two main shortcomings. On the one hand, it requires the temporal derivatives of the momentum and the position of the systems to calculate the loss function, which is difficult to obtain from real-world systems. On the other hand, HNN doesn’t strictly preserve the symplectic structure, because its symplectomorphism is realized by its loss function rather than its intrinsic network architecture. NeuralODE successfully bypasses the time derivatives of the datasets by incorporating an integrator solver into the network architecture.
Embedding the Hamiltonian prior into the NeuralODE, a series of methods are proposed, such as SRNN, SSINN, and TaylorNet, to predict the continuous trajectory of system variables; however, presently these methods are only designed to solve separable Hamiltonian systems. Instead of updating the continuous dynamics by integrating the neural networks in NeuralODE, SympNet adopts a symplectomorphism composed of well-designed both linear and non-linear matrices to intrinsically map the system variables within neighboring time steps. However, the parameters scale in the matrix map for training N dimensional Hamiltonian system in SympNet is O(N2), which makes it hard to generalize to the high dimensional N-body problems. For example, in Section 6, we predict the dynamic evolution of 6000 vortex particles, which is challenging for the training process of the SympNet on the level of O(60002).
NSSNN overcomes the weaknesses mentioned above. Under the framework of NeuralODE, NSSNN utilizes continuously-defined dynamics in the neural networks, which gives it the capability to learn the continuous-time evolution of dynamical systems. Based on Tao (2016), NSSNN embeds the symplectic prior into the nonseparable symplectic integrator to ensure the strict symplectomorphism, thereby guaranteeing the property of long-term predictability. In addition, unlike SympNet, NSSNN is highly flexible and can be generalized to high dimensional N-body problems by involving the interaction networks (Sanchez-Gonzalez et al., 2019), which will be further discussed in Section 6.

heading: 5.2 EXPERIMENTS
text: We compare five implementations that learn and predict Hamiltonian systems. The first one is NeuralODE, which trains the system by embedding the network fθ → (dq/dt, dp/dt) into the Runge-Kutta (RK) integrator. The other four, however, achieve the goal by fitting the Hamiltonian Hθ → H based on (1). Specifically, HNN trains the network with the constraints of the Hamiltonian symplectic gradient along with the time derivative of system variables and then embeds the welltrainedHθ into the RK integrator for predicting the system. The third and fourth implementations are ablation tests. One of them is improved HNN (IHNN), which embeds the well-trainedHθ into the nonseparable symplectic integrator (Tao’s integrator) for predicting. The other is to directly embed Hθ into the RK integrator for training, which we call HRK. The fifth method is NSSNN, which embedsHθ into the nonseparable symplectic integrator for training. For fair comparison, we adopt the same network structure (except that the dimension of output layer in NeuralODE is two times larger than that in the other four), the same L1 loss function and same
size of the dataset, and the precision of all integral schemes is second order, and the other parameters keep consistent with the one in Section 4. The time derivative in the dataset for training HNN and IHNN is obtained by the first difference method
dq dt ≈ q(Ttrain)− q(0) Ttrain and dp dt ≈ p(Ttrain)− q(0) Ttrain . (8)
Figure 5 demonstrates the differences between the five methods using a spring systemH = 0.5(q2 + p2) with different time span Ttrain = 0.4, 1 and same time step dt = 0.2. We can see that by introducing the nonseparable symplectic integrator into the prediction of the Hamiltonian system, NSSNN has a stronger long-term predicting ability than all the other methods. In addition, the prediction of HNN and IHNN lies in the dataset with time derivative; consequently, it will lead to a larger error when the given time span Ttrain is large.
Moreover, the datasets obtained by (11) in HNN and IHNN are sensitive to noise. Figure 6 compares the predictions of (q, p) for the system H = 0.5(q2 + 1)(p2 + 1), where the network is trained by the dataset with noise ∼ 0.05U(−1, 1), along with time span Ttrain = 0.2 and time step dt = 0.02. Under the condition with noise, NSSNN still performs well compared with other methods. Also, we compare the convergent error of a series of Hamiltonian systems with differentH trained with noisy data in Appendix D, which generally shows better robustness than HNN does.

heading: 6 MODELING VORTEX DYNAMICS OF MULTI-PARTICLE SYSTEM
text: For two-dimensional vortex particle systems, the dynamical equations of particle positions (xj , yj), j = 1, 2, · · · , Nv with particle strengths Γj can be written in the generalized Hamiltonian form as
Γj dxj dt = −∂H p ∂yj , Γj dyj dt = ∂Hp ∂xj , with Hp = 1 4π Nv∑ j,k=1 ΓjΓk log(|xj − xk|). (9)
By including the given particle strengths Γj in Algorithm 1, we can still adopt the method mentioned above to learn the Hamiltonian in (9) when there are fewer particles. However, considering a system with Nv 2 particles, the cost to collect training data from all Nv particles might be high, and the training process can be time-consuming. Thus, instead of collecting information from all Nv particles to train our model, we only use data collected from two bodies as training data to make predictions of the dynamics of Nv particles.
Specifically, we assume the interactive models between particle pairs with unit particle strengths Γj = 1 are the same, and their corresponding Hamiltonian can be represented as network Ĥθ(xj ,xk), based on which the corresponding Hamiltonian of Nv particles can be written as (Battaglia et al., 2016; Sanchez-Gonzalez et al., 2019)
Hpθ = Nv∑ i,j=1 ΓjΓkĤθ(xj ,xk). (10)
We embed (10) into the symplectic integrator that includes Γj to obtain the final network architecture.
The setup of the multi-particle problem is similar to the previous problems. The training time span is Ttrain = 0.01 while the prediction period can be up to Tpredict = 40. We use 2048 clean data samples to train our model. The training process takes about 100 epochs for the loss to converge. In Figure 7, we use our trained model to predict the dynamics of 6000-particle systems, including Taylor and Leapfrog vortices. We generate results of Taylor vortex and Leapfrop vortex using NSSNN and HNN and compare them with the ground truth. Vortex elements are used with corresponding initial vorticity conditions of Taylor vortex and Leapfrop vortex (Qu et al., 2019). The difficulty of the numerical modeling of these two systems lies in the separation of different dynamical vortices instead of having them merging into a bigger structure. In both cases, the vortices evolved using NSSNN are separated nicely as the ground truth shows, while the vortices merge together using HNN.

heading: 7 LIMITATIONS
text: The network with the embedded integrator is often more time-consuming to train than the one based on the dataset with time derivative. For example, the ratio of training time of the methods HNN and NSSNN is 1 : 3 when dt = Ttrain, and the training time of the recurrent networks further increases with the decreasing of dt. Although a smaller dt often has higher discretization accuracy, there is a tradeoff between training cost and predicting accuracy. Additionally, a smaller dt may potentially cause gradient explosion. In this case, we may want to use the adjoint method instead. Another limitation lies in the assumption that the symplectic structure is conserved. In real-world systems, there could be dissipation that makes this assumption unsatisfied.

heading: 8 CONCLUSIONS
text: We incorporate a classic ideal that maps a nonseparable system to a higher dimensional space making it quasi-separable to construct symplectic networks. With the intrinsic symplectic structure, NSSNN possesses many benefits compared with other methods. In particular, NSSNN is the first method that can learn the vortex dynamical system, and accurately predict the evolution of complex vortex structures, such as Taylor and Leapfrog vortices. NSSNN, based on the first principle of learning complex systems, has potential applications in fields of physics, astronomy, and weather forecast, etc. We will further explore the possibilities of neural networks with inherent structure-preserving ability in fields like 3D vortex dynamics and quantum turbulence. In addition, we will also work on general applications of NSSNN with datasets based on images or other real scenes through automatically identifying coordinate variables of Hamiltonian systems based on neural networks.

heading: ACKNOWLEDGMENTS
text: This project is supported in part by Neukom Institute CompX Faculty Grant, Burke Research Initiation Award, and ByteDance Gift Donation. Yunjin Tong is supported by the Dartmouth Women in Science Project (WISP), Undergraduate Advising and Research Program (UGAR), and Neukom Scholars Program.

heading: A NETWORK ARCHITECTURE
text: Figure 8(a) shows the forward pass of NSSNN is composed of a forward pass through a differentiable symplectic integrator as well as a backpropagation step through the model. Figure 8(b) plots the schematic diagram of NSSNN. For the constructed network Hθ(q,p), we integrate (4) by using the second-order symplectic integrator (Tao, 2016). Specifically, The input layer of the integrator is (q,p,x,y) = (q0,p0, q0,p0) at t = t0 and the output layer is (q,p,x,y) = (qn,pn,xn,yn) at t = t0 + ndt. The recursive relations of (qi,pi,xi,yi), i = 1, 2, · · · , n, are expressed by the algorithm 1.

heading: B SYMPLECTOMORPHISMS
text: One of the most important features of the time evolution of Hamilton’s equations is that it is a symplectomorphism, representing a transformation of phase space that is volume-preserving. In the setting of canonical coordinates, symplectomorphism means the transformation of the phase flow of a Hamiltonian system conserves the symplectic two-form
dq ∧ dp ≡ N∑ j=1 (dqj ∧ dpj) , (11)
where ∧ denotes the wedge product of two differential forms. The rules of wedge products can be found in Lee (2010). In the two-dimensional case, (11) can be understood as the area element of the surface. In this case, the symplectomorphism can be interpreted as the area element of the surface is constant. As proved below, our constructed network structure intrinsically preserves Hamiltonian structure.
Theorem B.1. For a given δ, the mapping φδ1, φδ2, and φδ3 in (5) are symplectomorphisms.
Proof. Let
(tqj , t p j , t x j , t y j ) = φ δ j(q,p,x,y), j = 1, 2, 3. (12)
.
From the first equation of (5), we have
dtq1 ∧ dt p 1 + dt x 1 ∧ dt y 1 =dq ∧ d [ p− δ ∂Hθ(q,y)
∂q
] + d [ x+ δ
∂Hθ(q,y) ∂p
] ∧ dy
=dq ∧ dp+ dx ∧ dy + δ [ ∂Hθ(q,y) ∂q∂y − ∂Hθ(q,y) ∂y∂q ] dq ∧ dy
=dq ∧ dp+ dx ∧ dy.
(13)
Similarly, we can prove that dtq2 ∧ dt p 2 + dt x 2 ∧ dt y 2 = dq ∧ dp+ dx∧ dy. In addition, from the third equation of (5), we can directly deduce that dtq3 ∧ dt p 3 + dt x 3 ∧ dt y 3 = dq ∧ dp+ dx ∧ dy.
Suppose that Φ1 and Φ2 are two symplectomorphisms. Then, it is easy to show that their composite map Φ2 ◦ Φ1 is also symplectomorphism due to the chain rule. Thus, the symplectomorphism of algorithm 1 can be guaranteed by the theorem B.1.
C DETERMINING COEFFICIENT ω
To further elucidation, the HamiltonianHA+HB without the binding, i.e.,H with ω = 0, in extended phase space (q,p,x,y) may not be integrable, even if H(q,p) is integrable in the original phase space (q,p). However,HC is integrable. Thus, as ω increases, a larger proportion in the phase space for H corresponds to regular behaviors (Kolmogorov, 1954). For H(q, p) = (q2 + 1)(p2 + 1)/2, shown in Fig. 9, we compare the trajectories starting from [q(0), p(0), x(0), y(0)] = (−3, 0,−3, 0) calculated by the symplectic integrator (Tao, 2016) with different ω, where the calculation accuracy is second order accuracy and the time interval is 0.001. As Figs. 9(a), (b), (c), and (d) shown, the chaotic region in phase space is significantly decreasing until forming a stable limit cycle. We define = ‖(q, p)− (x, y)‖2 as the calculation error of this system, shown in Figs. (e), (f), (g), and (h) that the error is decreasing with ω increasing, which fits the quantitative results of phase trajectory well.

heading: D OTHER EXPERIMENTS
text: We consider the pendulum, the Lotka–Volterra, the Spring, the Hénon–Heiles, the Tao’s example (Tao, 2016), the Fourier form of nonlinear Schrödinger and the vortex particle systems in our implementation. The Hamiltonian energies of these systems (except vortex particle system) are summarized as follows:
Pendulum system: H(q, p) = 3(1−cos(q))+p2. Lotka–Volterra system: H(q, p) = p−ep+2q−eq . Spring system: H(q, p) = q2 + p2. Hénon–Heiles system: H(q1, q2, p1, p2) = (p21 + p22)/2 + (q21 + q22) + (q 2 1q2 − q32/3)/2. Tao’s example (Tao, 2016): H(q, p) = (q2 + 1)(p2 + 1)/2. Fourier form of
nonlinear Schrödinger equation: H(q1, q2, p1, p2) = [ (q21 + p 2 1) 2 + (q22 + p 2 2) 2 ] /4−(q21q22 +p21p22− q21p 2 2 − p21q22 + 4q1q2p1p2).
The network is trained by the dataset with noise∼ 0.1U(−1, 1). The training time span, integral time step, and validation time span are 0.01, 0.01, and 0.1, respectively. Table 2 compares the Hamiltonian deviation H = ‖H(qtruth,ptruth)−H(qpredict,ppredict)‖2/‖H(qtruth,ptruth)‖2 and the prediction error p = ‖qtruth − qpredict‖1 + ‖ptruth − ppredict‖1. It is clearly from the Table 2 that NSSNN either outperforms or has similar performances as NeuralODE and HNN do.

Review summaries:
The paper extends the symplectic family of network architectures towards modeling nonseparable Hamiltonian dynamic systems. More specifically, the paper implements a  symplectic integration schema (from Tao (2016)) for solving arbitrary nonseparable (and separable) Hamiltonian systems within a symplectic neural network architecture. The results from several modeling tasks show that the proposed Nonseparable Symplectic NNs (NSSNNs) are more robust and accurate than vanilla HNNs and NeuralODEs when applied to nonseparable Hamiltonian systems. 

This paper describes a deep learning approach for predicting Hamiltonian systems. The original paper enforces conservation in the loss function. Several of the follow-up papers embed a symplectic integrator instead, but these couldn't handle non-separable systems. This paper can both handle non-separable systems and use a symplectic integrator to enforce conservation. They demonstrate their system on quite a few examples and show lower error (both in prediction and the deviation in the Hamiltonian) than a NeuralODE or the original HNN paper. The final example, in Figure 5, shows a compelling visual improvement.