Unifying physical systems’ inductive biases in neural ODE using dynamics constraints
Authors: Yi Heng Lim
Muhammad Firmansyah Kasim
Venue: Accepted by TMLR ()

Abstract: Conservation of energy is at the core of many physical phenomena and dynamical systems. There have been a significant number of works in the past few years aimed at predicting the trajectory of motion of dynamical systems using neural networks while adhering to the law of conservation of energy. Most of these works are inspired by classical mechanics such as Hamiltonian and Lagrangian mechanics as well as Neural Ordinary Differential Equations. While these works have been shown to work well in specific domains respectively, there is a lack of a unifying method that is more generally applicable without requiring significant changes to the neural network architectures. In this work, we aim to address this issue by providing a simple method that could be applied to not just energy-conserving systems, but also dissipative systems, by including a different inductive bias in different cases in the form of a regularisation term in the loss function. The proposed method does not require changing the neural network architecture and could form the basis to validate a novel idea, therefore showing promises to accelerate research in this direction.

Published in Transactions on Machine Learning Research (04/2023)
Unifying physical systems' inductive biases in neural ODE using dynamics constraints
Yi Heng Lim Machine Discovery Ltd. Ouford, United Kingdom
g.thenggmmachiine-dicoovvv.y...
Muhammad Firmansyah Kasim
outhors caattors@mmaccinneddcvveeerr...
Machine Discovery Ltd. Onford, United Kingdom
Reviewed on OpenReview: /ttps::/openreees net/ forum? id= 20Ab4976aa
Abstract
Conservation of energy is at the of many physical phenomena and dynamical systems. There have been significant number of works in the past few years aimed at predicting the trajectory of motion of dynamical systems using neural networks while adhering to the law of conservation of energy. Most of these works are inspired by classical mechanics such as Hamiltonian and Lagrangian mechanics as well as Neural Ordinary Differential Equations. While these works have been shown to work well specific domains respectively, there lack of unifying method that is more generally applicable without requiring significant changes to the neural network architectures. In this work, we aim to address this issue by providing simple method that could be applied to not just energy- energy-conserving systems, but also dissipative systems, by including a different inductive bias in different cases in the form of a regularisation term in the loss function. The proposed method does not require changing the neural network architecture and could form the basis to validate a novel idea, therefore showing promises accelerate rescarch this direction.
1 Introduction
Learning systems' dynamics with neural networks from data is one of the promising applications of deep neural networks in physical sciences. It has been shown that given enough data, neural networks can uncover the hidden dynamics of physical systems (Chen et al. 2018; Greydanus et al. 2019 Cranmer et al. 2020
Chen et al. 2019 and can be used simulate the systems with new conditions (Greydanus Sosanya 2022
The most straightforward way in learning system's dynamics is by training neural network to learn the system's dynamics given the states at the time, s(t), i.c.
$ ds dt
f(s(t))
(1)
Once the neural network is trained, the system's 's behaviour can be simulated by solving the equation above using ordinary differential equations solver (ODE). This method is typically called neural ordinary differential equations (NODE) (Chen et al. 2018) in the literature.
Despite its simplicity, NODE has a difficulty incorporating inductive biases of systems into training. For example, system is known to conserve energy, the NODE approach above typically results growing or shrinking of the energy
long period of time. A popular workaround is to learn the

---

Published in Transactions on Machine Learning Research (04/2023)
Hamiltonian /rryydamu et al. 2019 Toth et al., 2019 Sanchez-Gonzaler et al., 2019) or Lagrangian CCranmer et al. 2020) of the energy-conserving system with neural network, then get the dynamics as the derivatives of the learned quantity. Although this approach works in learning energy-conserving systems, it has been shown that to incorporate a specific inductive bias, it takes a significant effort to design a special noural network architecture or to oblain stales in special coordinate systems.
Another example is incorporating a bias if the system is known to be dissipative or losing energy. Like in the energy-conserving case, a lot of efforts go into designing special neural network architectures or choosing the optimal weights to ensure stability (Tuor et al. 2020 Xiong et al., 2021 ). Although those methods work well
in getting stable dynamics, changing neural network architectures could potentially reduce the expressiveness of the neural networks.
In this paper, we show how to incorporate various inductive biases in learning systems' dynamics by having various constraints on the dynamics function, including Hamiltonian systems and dissipative systems, without changing the architecture of the neural networks. Where necessary, we also include Invertible Neural Network to transform the coordinate system of the imput states, and enforce the inductive biases in the transformed coordinate system. Our code for this paper is available at
Our paper makes two main contributions in the form of two distinct regularisation terms for two types of systems, namely energy-conserving and dissipative systems:
Derivation of Equation 4 for energy-conserving systems, which is novel regularisation term that incorporates physical constraints into the training of neural networks for such systems.
2. Derivation of Equation 8 for dissipative systems, based on the argument that such systems are asymptotically stable. This is another novel regularisation term that improves the physical consistency of neural network models for dissipative systems.
2 Related works
Motivated by Hamiltonian mechanics, Hamiltonian Neural Network (HNN) was proposed in GGeyydams et al. 2019). HNN takes a set of positions and their corresponding general momenta as input, and models a scalar function called the Hamiltonian. The derivatives of the learned Hamiltonian with respect the inputs are then computed and multiplied by symplectic matrix, giving the learned state derivatives. Since the symplectic matrix is manually introduced, HNN has a nice property of observing symplectic structure through time, thus conserving energy. In practice, however, the applicability of HNN is limited by the fact that it requires input data to be prepared in canonical form, which is usually difficult to obtain. There have been a lot of follow-up works that build on the foundation of HNN and cxpand it to different settings including dissipative systems (Greydanus & Sosanya, 2022 Zhong et al. 2020), generative networks (Toth et al. 2019 and graph networks (Sanchez- Gonzalez et al. 2019), while others try to improve HNN by simplifying it Finzi Gruver al.] 2022).
let al. 2020
To tackle the shortcoming of HNN, another line of work based on Lagrangian mechanics, called Lagrangian Neural Network (LNN), was proposed in (Cranmer et al., 2020). As the name suggests, the Lagrangian is modelled by a neural network, and the expression to compute the second-order state derivatives could be obtained by algebraically rearranging the Euler-Lagrang, equation. With the introduction of Lagrangian prior in the neural network, LNN manages to conserve the total energy of a system without limiting the states to be strictly canonical. Since LNN assumes coordinate system in which states are casily measured in practice, can be applied to a wider range of problems where HNN fails. However, to obtain the second-order state derivatives, the inverse Hessian of the neural network must be computed, making LNN computationally expensive, susceptible to ill-conditioning, and sensitive to poor choice of activation function such as ReLU.
Neural Symplectic Form (NSF) (Chen et al. 2021) exploits the coordinate-free formulation of the Hamiltonian with symplectic 2-form. To speed up the computation, an elegant method to model the symplectic 2-form as an exterior derivative of a paramctorised 1-form was proposod. However, there are few drawbacks in this work. Firstly, without good understanding of exterior calculus and differential geometry, it could be difficult. to understand the motivation behind the work. Secondly, it is expensive to compute the inverse of
2

---

Published in Transactions on Machine Learning Research (04/2023)
the learnable skew-symmetric matrix which requires ((N3) operations, where 2N is the number of state variables. Thirdly, the training of NSF is unstable and tends to plateau very quickly, see Appendix A.1 To the best of our knowledge, the latter has not been addressed in any work to date.
Invertible Neural Networks (INNs) are neural networks that are elegantly designed such a way that the inverses of the neural networks are always available with just. a forward pass. They have found wide adoption in flow-based generative models and are typically used map data a latent space with a simpler distribution. (Dinh et al.) 2014 2016) propose a simple building block called affine coupling layer whose inverse is easily attainable, and an invertible neural network is obtained by stacking few of these affine coupling layers in sequence. In each affine coupling layer, the inputs are permuted since the INNs are not permutation invariant. (Kingma & Dhariwal 2018) further improves the model by introducing a learnable 1 x 1 invertible convolution to replace the fixed permutation of channel dimension used in (Dinh et al., 2016). In connection with NODE, (Zhi et al.) (201)) leverages INNs map the underlying vector field of an ODE system to a base vector field, and others (Chen al. 2022 in Hamiltonian systems.
Jin et al. 2022) have attempted include INNs
3 Methods
Consider a system with states s E R where its states' dynamics depend on the states at the given time, as written in equation [] The function f from equation [ can be represented with an ordinary neural network. The training can then be done by minimizing the loss function,
C = 9
+ (S).
(2)
The term & is the time-derivative of the observed states s from the experimental data, å is the computed we is the constraint's weight, and C(.) is the constraint applied to the dynamics function, $. The last term, ie. the constraint, is the term that can be adjusted based on the inductive bias of the system.
states' dynamics from equation [
3.1 Hamiltonian systems
Hamiltonian mechanics is framework to describe particles' equations of motion using equation in the phase space of the position and the generalised momenta p. This is in contrast to Newtonian's mechanics where it uses forces to describe the equations of motion. Using the Hamilton's 's equation, H(g, p), a set of equations of motion can obtained by
dq 0H dp 0p"
0H ōq
dt
By defining the state space vector as s (q" " p TT the states' dynamics of Hamiltonian systems can be written (Greydanus al.] 2019: Chen al. 2021
0 ! -I 0
$ (VH(s) where J
where H is the Hamiltonian of the system (lypically the energy in some cases) and VH is the gradient of the Hamiltonian with respect to s. The dynamics in the equation above preserve the quantity H, ie. dH /dt 0, which makes it suitable to simulate systems with constant Hamiltonian. If the dynamics are learned using Hamiltonian neural network (Greydanus et al. 2019), the dynamics will be equal to equation 3 The Jacobian in the Hamiltonian system can be written as 08/0s = J( H H0OOss). As the Hessian of the Hamiltonian (8² H (dsOs) a symmetric matrix, the Jacobian of a Hamiltonian system is a Hamiltonian matrix.
Being a Hamiltonian matrix, the Jacobian of the dynamics in equation 3 has special properties. One of them is having zero trace (V $ 0), which means that there is no source sink in the vector field of &.
3

---

Published in Transactions on Machine Learning Research (04/2023)
Another property is near the local minimum of (s) (if any), the eigenvalues of the Jacobian 08/0s are pure-imaginary which guarantees stability (see appendix A.2 for the derivations).
In order to get the dynamics as in equation 3 by learning the dynamics á directly, it has to be constrained so that the matrix J (0s/0s)) is symmetric. Therefore, we can apply the constraint
08 Us
0s
T
CH (8) = J"
(4)
Øs
where ||-||e is the Frobenius norm to get the matrix J (06/ 0s) as symmetric as possible. Note that J-1 J Theorem 3.1. The dynamics of the states, & E IR, follow the form in equation for function H(s) IR" - if and only if the matria J (a6/// is symmetric.
Proof. If the dynamics & follow equation 3 then the matrix J J- (08/Os) represents the Hessian of the Hamiltonian H, which is symmetric. To prove the converse, note that the matrix J ( - (0z /Os) can a Jacobian matrix of p = (-) where šy and š are the first half and the last half elements s. If the matrix (0p/Os) is a symmetric Jacobian matrix, then the expression below must be evaluated to zero,
Op: Opj 0s, 0s;
Opp (0pr asy
Σ
ds, A ds, = Σ 2.
ds; A ds, = df = 0,
where f= Pids, As df 0, by Poincaré's lemma, there exists an such that da f. In other words, if the matrix (8p/Os) is a symmetric Jacobian matrix, then there exists an a such that its Hessian is equal to the Jacobian matrix, i.e. (0a/0s0s)) (öp/Os). This corresponds to the Hamiltonian H in equation
3.2 Coordinate-transfomna Hamiltonian systems
One of the drawbacks of Hamiltonian systems that they have to use the right coordinates for the states, usually known as the canonical coordinates. In simple cases like simple harmonic motion of mass-spring system, the canonical coordinates easily be found (i.e. they the position and momentum of the mass). However, in more complicated cases, finding the canonical coordinates for Hamiltonian sometimes requires expertisc and analytical trial-and-crror.
Let us denote the canonical states' coordinates as z and the observed states' coordinates as s. As the canonical states' coordinates z are not usually known, we can find them using a learnable invertible transformation,
g(s) and g (z))
then learn the dynamics of with neural network, i.e. å = tz, This way, the dynamics of can be computed by $ = (0z)) The invertible transformation above can be implemented using an invertible neural network {Kiggma & Dhariwal 2018). With the invertible transformation between and the dynamics of s this case can be written as,
$ - G GSGG
(6)
where Hs(s) H and matrix G(s) 0g/Os the Jacobian of the transformation function from to from equation 5 It can also be shown that the dynamics in equation 0 preserve the quantity of H'.
Theorem 3.2. For following the equation b the quantity of H'(s) are constant throughout the time, i.e. dI'' 0.
Proof. Denote the matrix D as D = G- G-'JG-7 -T, where the matrix D is skow-symmetric, because
D" -G-
-D.
4

---

Published in Transactions on Machine Learning Research (04/2023)
The rate of change of H' can be written as (s) u"ss Du, where (H)(s). As the quantity above is a scalar, it must be equal to its transpose. Thus,
Du
[u Du = uu'D'u=-u Du
88 D -D because of its skew-symmetricity. Therefore, (s) =
As this system is similar to the Hamiltomian system in the previous subsection, we can apply the constraint below
CHr($) Cи (2)
where (-) is the constraint from equation 4
3.3 Dissipative systems
Dissipative systems have characteristic that their dynamics are asymptotically stable due to the frictions applied to the systems. In terms of ODE dynamics, this asymptotic stablility can be achieved by having the real parts of all eigenvalues of the Jacobian 08/Os to be negative. This would ensure stability around the stationary point. Therefore, the following constraint can be applied enforce the condition,
05 Ōs
2
CD(()= - Σ max O, Re A
di
where (:'sss are the eigenvalues of the Jacobian 08/0s, as's are the assumed upper bounds of the real parts of the eigenvalues. The upper bounds a, 's are sometimes useful in making sure that the eigenvalues are far from 0, further ensuring the stability.
4 Experiments
In order to see if the methods described above work, we tested them various physical systems.
Task Ideal mass-spring system. The first case is an ideal mass-spring system. The Hamiltonian canonical state of this system is (z, má)" where m is the mass and æ is the displacement of the mass from its equilibrium position. The acceleration given by -kk/m)ë. In our experiments, we set m and the spring constant 1, and the canonical state simply (a. We generated the data by randomly initialising the state from a normal distribution, and added white noise with o 0.1 to the state trajectory. For training, we generated 250 trajectories, with 30 samples in each trajectory equally spaced within t [0,2
Task 2 Ideal rod-pendulum. The second case is ideal rod-pendulum. The Hamiltonian canonical state of this system is (0, mo)T where m is the mass and is the angular displacement of the mass from its equilibrium position. In our experiment, we set m = 1, and the canonical state becomes (r. The angular acceleration is given by 0 = - (g/l) sin 0 with the value of g and set to 3 and respectively. Similar to Task 1, we generated the data by randomly initialising the state from a normal distribution, and added white noise with o 0.1 to the state trajectory. For training, again, we generated 250 trajectories, with 30 samples in each trajectory equally spaced within [0, 2m.. This task is slightly more complicated than Task 1 since it is non-linear dynamical system.
Task Ideal double rod-pendulum. The third case is an ideal double rod-pendulum. We chose a general state of (01, 02 0, 6²)" since the Hamiltonian canonical state this problem is non-trivial to obtain. Again, we set the the two masses m. = m2 = and the two rod lengths data by randomly initialising the state from a normal distribution and obtained the trajectory. With the aforementioned parameters, the two angular accelerations are given by
l2 1, and generated the
0 (sin - sin 202) - 2 (82) + 07) cos(O1 02)
3 - cos cos((2) 20,)
5

---

Published in Transactions on Machine Learning Research (04/2023)
02 2 0in - 02) - cos(0 )03 cos(O, 02))
3 - cos(20) - 282)
Since this a chaotic dynamical system, we omit the addition of noise to the trajectory to avoid introducing irreducible alcatoric uncertainty, and gencratod more data for training: 2000 trajectories, with 300 samples in each trajectory equally spaced within 1 - [0, 2x].
Task 4- Damped single rod-pendulum with redundant states. In this case the pendulum is damped with cocfficient a to model a dissipative system. The angular acceleration is given by b = (g/ll sin 0 ab and a = 0.05. The observed stales are (r, y, å ") where I and y are respectively the horizontal and vertical positions from the pivot, i.e. x == ssin and = -I cos 0. These states are selected instead of (0, d) to illustrate a system where the observed states are redundant. As this is a dissipative system, we use the constraint from equation g for the training.
with m = =
4.1 Training details
For Task 1 and Task 2, we define our method as being a simple Neural Ordinary Differential Equation (NODE), but with constraint 4 added to the loss function. For Task 3, since the input states were not in canonical form, we added a component of Invertible Neural Network (INN), with the FrELA package, based on the invertible architecture defined in KKnggaa & Dhariwal 2018) to transform the state coordinate system in our method, and applied constraint 7 For Task 4, since it is dissipative system, we added constraint g to the loss function.
In the first three lasks, we compared our method with baseline NODE, HNN, LNN and NSF. For Task 4, since is well established that HNN, LNN and NSF would not work in dissipative systems, we only compared our method to baseline NODE. We trained all models in all tasks with Adam optimiser (Kingma & Ba 2014 and the same learning rate of 10-4, We constructed multi-layer perceptrons (MLPs) for all models, with 3 layers of 200 hidden units in each case. INN with 8 blocks was used where needed, with each block having 2 layers of 100 hidden units. We chose softplus as the activation function for all models as it worked well in all cases.
In our experiments, LNN and NSF were difficult to train and training inconsistencies were often incurred by these two models. The training inconsistencies of LNN were pointed out in (Chen al.] 2021) and we showed the inconsistencies of NSF in Appendix A.I]] In those cases, we ran the experiments several times to obtain the best training results that we could compare our method to the best trained models. We found that using a smaller batch size to increase the number of training steps helped in reducing those inconsistencies. We note that this is a weakness of LNN and NSF as the inconsistencies might make training difficult in practical applications.
For our method, we set the cocfficient of the constraint as large as possible without compromising the loss of time derivatives of states. In Task 1, and 4 respectively, we set the coefficient to be 10 104, 10x and 10a. We used batch size of 32 for Task 1, Task 2 and Task 4, and larger batch size of 1280 for Task 3 since there were more data. We ran all experiments for 1000 epochs a single NVIDIA Tesla T4 GPU.
To test, we started from a randomly generated state at time O, and rolled out the subsequent states from time O to 100, and logged the RMSE of the energy deviations. Note that the initial state at testing is different from the training data. We recorded the median, 2.5th percentile and 97.5th percentile since there are outliers in some of the test results.
5 Results
Overall, our method performed consistently well and ranked in the top two methods in all tasks in terms of energy drift from the ground truth. The other methods returned mixed results, performing well in some tasks but not the others. While NSF performed reasonably well in all tasks, we show in Appendix A.1 that training inconsistenccies are often encountered in NSF.
http:///ghhubbcmm//iiiH////
6

---

Published in Transactions on Machine Learning Research (04/2023)
Case
Task 1
Task 2
Task 3
Task 4
NODE
0.086*0.095
00.077
DC
00
0.18 718.10
0.059
HNN
0.0065 -0.000000000 0.0045
+0.013
0.015
1.015+00066
24.44+88.9
N/A
0.012
-11.40
LNN
0.0060 021
0.03220.02
7.92+121141
N/A
00-16
NSF
0.014 0.0093
0.027 0.016
4.29
+5500
-3.05 33.05
N/A
Our method
0.001 1.017+0.000 0,0069
0.02410.063
06008
4.03
+10.52
0.12 -0.012
+0.39
0.0012
0.013
3.16
-0.012
Table 1: RMSE of the energy deviations of 100 test trajectorics. The numbers shown are mcan + standard deviation. Since there are outliers in some of the results, the median, 2.5th and 97.5th of the RMSE in each case are logged. The best result in each case, judging by the 97.5th the RMSE, bolded, and the second best result is highlighted in blue. oo indicates numerical overflow.
Our method performed the best Task and was on par with HNN in Task 2 which logged the best performance. This is consistent with our expectation since the input coordinates are in canonical form. It is worth noting that even though the symplectic structure is enforced inherently in the HNN architecture while our method employed a soft constraint, our method still compared competitively with HNN in the first two tasks. More importantly, as can be seen in Figure [a) and (b), similar to the other known energy-conserving models, the periodic drift in energy of our model oscillates around a constant offset from the x-axis throughout the time span. While this periodic drift is mostly attributed to the state prediction error due to the white noise added to the training data, the constant offsct from the z-axis mcans that the average perceived energy in the system is conserved.
A more interesting point to note is that our method outperformed HNN in Task 1, even though in both models the symplectic structure is preserved. We suspect this is because the inherent symplecticity in HNN acts as a double-edged sword: although it is good at preserving the symplecticity, it reduces the flexibility of the model in some cases especially in the presence of noise.
Task 3 is the most difficult system to model among the first three tasks since the motion double- is chaotic and the input coordinates are not in canonical form. Table [ shows that our method performed the best in terms of energy drift of the predicted system. This is further corroborated by Figure [e) where again, the energy drift is at constant periodic offset from the r-axis.
In Task 4, our method performed significantly better than vanilla NODE, It is interesting to note that while NODE failed spectacularly in the edge cases as indicated by the large RMSE upper bound in terms of energy drift, our method worked reasonably well by just including a simple dissipative bias in the loss function. Figure [d) shows one of the worse predictions that NODE returned, and the corresponding states are shown in Figure and Figure
The states in Task are redundant there are 4 states used while there are only 2 true states. This makes the dynamics of the pendulum lie on a low dimensional manifold in a 4-dimensional space. Without the constraint, the error in learning the dynamics by NODE could send the states outside the manifold where no training data exists, thus making the dynamics unstable. On the other hand, the constraint in equation makes sure the dynamics are stable around some stationary points which typically lie on the manifold. Therefore, a slight error in learning the dynamics our case can be corrected by attracting it back to the stationary point on the manifold. This is an important example of ensuring the stability of a novel system since the true number of states are usually unknown.
6
Pixel pendulum
6.1
Energy-conserving pendulum
Similar to GGryddans et al. 2019: and (Chen et al.) 2021), we tested our model with pixel data. This resembles our model in Task 3 where INN is used to transform the state coordinates, except that in this case
7

---

Published in Transactions on Machine Learning Research (04/2023)
Energs deviation (mass sprting
Iowegy deviation (sinde-penfolom)
Energy deviation (dovble pendalum!
Energy deviation (domped pendulum
140
FIDE
14
0.35
2.0
0.30
NVVV
0.08
0/2
1.0
04004 4444A44
03
00-
00
0.0
6.35
-1.5
(a.)
(b)
(c)
(d)
Figure 1: The energy drift from ground truth in every model for (a) mass-spring (Task 1), (b) single-peddllmm
(Task 2), (c) double-pendulum (Task 3) and (d) damped-pendulne (Task 4). All test cases were taken
from the same distribution as the training data, except for (c) where the test case was sampled from 8 0.5.5 0.5), Vi for better illustration.
(a)
(b)
(c)
(d)
Figure/ 2: The trajectories of the four states of
in Task This sample is picked randomly from one of the edge cases in NODE. (a) shows the trajectory of sin (b) shows the trajectory of cos 0
while (c) and (d) are the trajectories of the horizontal and vertical velocities respectively.
Angle (damped pendulum)
Angular velecity idamped pandulum!
TRUE
2.5
NUUs: HBo
N00e
Morg dynn
-0.
4.0
SIRAEE
ACCEEI
-.5
00
200
60
300
Tims
(a)
(b)
Figure 3: The trajectories of and 0 of daamped-pendulum in Task 4, transformed from Figure 2

---

Published in Transactions on Machine Learning Research (04/2023)
Latent phase pict (Node)
Latent Latentphase pict (H VN)
Latent phase siot (Our method;
031
01
01
0.
-.1
-07-
-03
(a)
(b)
(c)
Figure 4: Phase plot of the learned dynamics in the encoded latent space. Trained (a) vanilla NODE, (b) HNN and (c) our method. The learned dynamics here are the learned states, i.e. angular position and angular velocity of the mass, rolled out in time.
Figure 5: Predictions of the pendulum position with image data from time 0 to 100.
we seek to encode image data in a latent coordinate system. We used the same training hyperparameters as the previous experiments, except that regularisation weight of 500 is used here for our model and 160 videos of 100 frames equally spaced within [o, 100) each were used as the training data. We also found that ReLU activation function worked better for the autoencoder. Figure 5 shows the predictions for = [0, 100] using NODE, HNN, and our method while Figure shows the phase plots in the encoded latent space. The phase plot of NODE, expected, spirals inwards indicating a drift in the symplectic structure, and Figure 5 shows that its predicted pendulum dynamics are barely moving due to the energy drift. On the other hand, HNN and our method performed similarly well with the predicted pendulum dynamics closely resemble the ground truth. The phase plots of both HNN and our method are perfect circles, meaning there is no sink or source in the latent vector field and the symplectic structure is preserved. This is, however, simple proof of concept with a vanilla autoencoder. We note that going forward, we could improve the results by using a better autoencoder or even a pretrained model.
6.2 Damped pendulum
To Lest our model for dissipative systems with image data, we modified the code from OpenAl pendulum gym Brockman et al., 2016) to include a damping factor of 0.05 to the pendulum motion. In this case, we set the to mirror the fact that the number of states in system usually unknown, and added constraint X weighted by 500 the vanilla NODE in our method. In Figure 6 the predictions of NODE noisy due to large errors in latent space. As explained in Task 4, these errors due to the
number of latent states
9

---

Published in Transactions on Machine Learning Research (04/2023)
Ouretted
Figure Predictions of the damped pendulum position with image data from time 100.
predicted states being outside of the 2D manifold. In contrast, with simple constraint, successfully reduced these errors and removed the noise from the pixel predictions.
7 Limitations
Even though our method is simple and proven to work well in different scenarios including both energy- conserving and dissipative systems, we have introduced free hyper-parameter which is the regularisation weight. This weight could be difficult to fine-tune depending on the experiments: setting too low value may not enforce the constraint while setting too high a value may jeopardise the training. This fine-tuning process is further exacerbated if we have more than one regularisation term. As part of the future work, we may explore using Lagrangian multipliers as a potential solution.
8 Broader impact
Unlike classical control theory that has a strong theoretical foundation to explain every phenomenon and trajectory deterministic dynamical system, the results obtained from neural networks could be inexplicable times.
9 Conclusions
In this paper, we proposed simple method enforce inductive bias in a vanilla Neural Ordinary Differential Equation (NODE) to model both energy-conserving and dissipative systems. This was conducted by incorporating regularisation term in the loss function: in energy-conserving systems the Hamiltonian symplectic structure is enforced whereas in dissipative systems the state vector fields are made to spiral towards a sink. Our method unifies different approaches in modelling dynamical systems, and is more generally applicable since it does not. require significant changes to the neural network architectures, nor does it require data specific coordinate systems, thus showing promises to accelerate research in this domain.
10

---

Published in Transactions on Machine Learning Research (04/2023)
References
Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojeiech Zaremba.. Openai gym, 2016.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Yuhan Chen, Takashi Matsubara, and Takaharu Yaguchi. Neural symplectic form: Learning hamiltonian equations on general coordinate systems. Aduances in Newral Information Processing Systems, 34, 2021.
Yuhan Chen, Takashi Matsubara, and Takaharu Yaguchi. Kam theory meets statistical learning theory: Hamiltonian neural networks with non-zero training loss. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pp. 0322-6332, 2022.
Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. In International Conference on Learning Representations, 2019.
Miles Crammer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXiv preprint orXiv:2003.0400 2020.
Laurent Dinh, David Krueger, and Yoshua Bengio. Nice: Non-linear independent components estimation. arXiv preprinl ar Xrx: 140.85566 2014.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp. arXiv preprint 1rXiv:1605.0880 2016.
Mare Finzi, Ke Alexander Wang, and Andrew G Wilson. Simplifying hamiltonian and lagrangian neural networks via explicit constraints. Advances in neural information processing systems, 18:13880-13888 2020.
Sam Greydanus and Andrew Sosanya. Dissipative hamiltonian neural networks: Learning dissipative and conservative dynamics separately. arXiv preprint arXi::200111088, 2022.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in Neural Information Processing Systems, 32, 2019.
Nate Gruver, Marc Finzi, Samuel Stanton, and Andrew Gordon Wilson. Deconstructing the inductive biases of hamiltonian neural networks. arXiv preprint arXiv:2202.04880 04836, 2022.
Pengzhan Jin, Zhen Zhang, Ioannis Kevrekidis, and George Em Karniadakis. Learning poisson systems and trajectorics of autonomous systems via poisson neural networks. IEEE Transactions on Neural Networks and Learning Systems, 2022.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv:1412.6900 2014.
arXiv preprint
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1xl convolutions. Advances in neural information processing systems, 31, 2018.
Alvaro Sanchez-Gonzaler Victor Bapst, Kyle Cranmer, and Peter Battaglia, Hamiltonian graph networks with ode integrators. arXiv preprint 19X9v 2019.
Peter Toth, Danilo Jimenez Rezende, Andrew Jaogle, Sébastien Racanière, Alcksandar Botev, and Irina Higgins. Hamiltonian generative networks. arXiv preprint arXiv:1909.13789, 2019.
Aaron Tuor, Jan Drgona, and Draguna Vrabic. Constrained neural ordinary differential equations with stability guarantees. arXiv preprint Xiix:2004,10883, 2020.
Jie Xiong, Alan Yang, Maxim Raginsky, and Elyse Roscnbaum. Neural networks for transient modeling of circuits. In 2021 ACM/IEEE 3rd Workshop on Machine Learning for CAD (MLCAD), pp. 1-7. IEEE, 2021.
11

---

Published in Transactions on Machine Learning Research (04/2023)
Weiming Zhi, Tin Lai, Lionel Ott, Edwin V Bonilla, and Fabio Ramos. Learning odes via diffeomorphisms for fast and robust integration. arXiv preprint 0000:2207.01650, 2021.
Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Dissipative symoden: Encoding hamiltonian dynamics with dissipation and control into deep learning. Xiv preprint arXi::2002.08660 2020.
12

---

Published in Transactions on Machine Learning Research (04/2023)
Training Losses
Validation Losses
Energy ssinge--pendulum
3.0
1.00
2.5
TRUE
CONVERGED
2.5
FAILED
0.50
2.0
0.25
1.5
0.00
1.0
1.0
-0.50
0.5
0.5
-0.75
0.0
-1.00
Epoch
Epoch
Time
(a)
(b)
(c)
Figure 7: (a) Training losses of 10 identical NSF runs, (b) Validation losses of 10 identical runs, and (c)
Comparison of energy conservation between failed and successful NSF runs.
A Appendix
A.1 NSF training inconsistency
To show the training inconsistency often incurred by NSF, we repeated the training of Task 2 with NSF 10 times using the same data and plotted the losses of each training, see Figure [ We used single batch that included all data, and learning rate of 10-3 here. We note that this problem could be mitigated, but not removed, by using smaller batch size, which we did in the main paper. We also note that this inconsistency is common to all tasks. The early plateaus in some of the runs are not due to the expressivity of the neural networks, since the exact same configuration was used in all runs. Even though NSF, when converged, gives relatively good performance in terms of energy conservation, this training inconsistency unacceptable in any practical application.
A.2 Stability of Hamiltonian near a local minimum
The dynamics of Hamiltonian system equation g can be written as
= (s) where = -I !)
0
(9)
The Jacobian matrix of the dynamics 05/ Ös is then
8 (s)) asds
JB where
(10)
ds
where 82 H() /ösds the symmetric Hessian matrix of the Hamiltonian. Near local minimum of Hs),, the Hessian matrix B is positive definite matrix. With positive definitive property of B, there exists the matrix square root and its inverse, i.e., 1/2 and B-2 which are also symmetric matrices. Therefore, the matrix JB is similar to B BI/2
as
BI//B1/2 BJ)
(11)
As J is a skew-symmetric matrix and B is symmetric matrix, the matrix B1/2JB1/2 skew-symmetric matrix.
Being a skew-symmetric matrix, the cigenvalues of 1//2JB1/2 are pure imaginary. Therefore, from the similarity, we can conclude that the eigenvalues of JB are also pure imaginary near the local minimum of (s)).
13

---

Published in Transactions on Machine Learning Research (04/2023)
Energy (single-pendulum)
1.0
TRUE
0.8
RODE+REG
0.6
0.4
0.2
0.0
-0.2
-0.4
Time
Figure 8: Comparison between HNN and our method in Task 2 when the models are rolled out in a longer time span during inference.
018
0.0
0.5
0.88
80-
0.0
26-3
(a)
(b)
(c)
(d)
Figure 9: (a) Training losses of our method for (a) mass-spring (Task 1), (b) single-pendulum (Task 2), (c) double-pendulum (Task 3) and (d) damped--endulum (Task 4). Logarithmic loss is plotted against the number of training steps in each plot.
A.3 Additional plots
We include Figure show the reliability of method during inference when the prediction spans long time horizon.
We also include Figure to show the training curves of our method for Tasks 1 to 4.
14


Reviews:
[{"review_id": "PpQhYFR9H2", "reviewer": "Reviewer_jz52", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is quite easy to understand and I like the simplicity of the proposed method which seems to perform well across the proposed tasks. I can not judge the novelty of the method, here i would like to hear the other reviewers opinions. The method feels like one of the first things to try. \nAlthough I think the mathematical presentation is sufficient and clear, a few of experimental details would be important to show in the paper. Here, there are a few things that need clarification in order for acceptance in my opinion, see requested changes. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose to use a simple theoretically motivated regularisation term when modelling neural ODEs in order to conserve energy of the systems. This works with and without change of the coordinate systems for which the authors propose to learn invertible neural networks. They also provide a regularisation term to learn to model dissipative systems. \n\nOn a set of dynamical systems as well on a task which requires inferring states from pixel data, they test their method against related work and show promising results. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Can you provide training curves of your method on the tasks? \n2. Can you provide details on how the hyperparameters especially the regularization strength influences the behaviour? \n3. It would be interesting to see the influence of the network depth and width on the performance/generalisation. \n3. Why is the regularization strength so high? The loss function seem not correctly normalized? As far as explained you are using the squared-error loss, why not the MSE?\n4. Can you please revisit the section about training instabilities of the NSF method? As far as I can see, the method does not suffer from instabilities (divergence) but just does not converge for some of the runs. I think it is a bit far fetched to say that the method suffers from instabilities. These slow training dynamics could e.g. be easily be due to improper tuning of the network initialization or optimizers. If you really want to make a point on this, please provide more details and tune all relevant hyperparameters. \n5. What kind of network initialisation are you using?\n6. Can you provide pseudo-code for your method?\n7. It would be interesting to get a feeling of the computational overhead of the method / reg term. Are these Jacobians/Hessians expensive to compute?\n8. Can you compare the training times (in sample efficiency) of your with related work?\n", "broader_impact_concerns": "-", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "lrHSidpmw1", "reviewer": "Reviewer_vbRE", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n\nThe proposed method of adding a regularization term for inductive bias is easy to implement.\n\nIn the experiments presented in the paper, the proposed method has good performance, compared to other similar methods.\n\nWeakness: \n\n> The idea of adding regularization to penalize unwanted aspects of learning is a very common trick in the machine learning community. Hence, the novelty is limited. \n\n> The experiments shown in the paper are only very simple physics systems, which have very simple physics solutions. For a good evaluation of the proposed method, I expect it is tested on more complicated systems. \n\n> I have doubts on the effectiveness of the proposed method. I noticed that all the experiments use super large constraint coefficients so that the enforced constraints are strictly obeyed. According to my personal experience, there is always a tradeoff between learning the system (corresponding to the loss without regularization) and enforcing the regularization. When the constraint coefficient is so large as in this paper, the machine learning models usually learn little to nothing meaningful. Hence, I suspect that the good performance in the presented experiments might be because the physics systems are so simple that they are very easy to learn. If this is the case, the proposed method may not work well for complicated systems.\n\n> Another doubt I have is that: when there exist several inductive biases (for example several conservation laws), there should be multiple regularizers/constraints in the loss function. In that case, tuning the coefficients (which are hyper-parameters) becomes very hard, and often not doable. Hence, I doubt the effectiveness of the method in these cases. (I think this should be a very common situation in physics, as there are many conservation laws or other rules).\n\n> [Minor] The paper claims “changing neural network architectures could potentially reduce the expressiveness of the neural network”. I am not convinced, I hope the authors provide some evidence for this claim. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper considers the problem of learning physical systems’ dynamics by machine learning models.  It proposes to use the regularization term in loss function to incorporate inductive biases which are often inferred by physics, such as conservation laws. The regularization term penalizes dynamics that are away from the physical laws. The paper also experimentally evaluated the proposed method on several simple physics systems. The performance of the proposed method outperforms several other related methods in most of the examined systems.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I think the authors need to evaluate the method in more complicated physics systems, and to show it works well in these scenarios. \n\nI hope the authors can show an efficient way of managing the scenarios of multiple inductive biases and constraints.\n\nI would like the authors to also provide evidence for “it has been shown that to incorporate a specific inductive bias it takes a significant effort ……”\n\nProvide evidence for the claim “changing neural network architectures could potentially reduce the expressiveness of the neural network”.\n", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "5ePN3PcOLv", "reviewer": "Reviewer_rNfb", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\n- The authors clearly explain the challenge of including physics-based constraints and inductive biases and give a good overview of existing methods.\n- The chosen tasks are explained well, and seem like appropriate simple test cases.\n\n\n# Concerns\n\n- It is a little unclear exactly what theorems and methods are novel contributions. I would suggest adding an itemized list of contributions at the end of the introduction.\n- I think the assumed physics knowledge in this paper is going to be a challenge for a large portion of the TMLR audience. I realize that space is at a premium, but a very short introductory example in section 3.1 would be helpful--describing the Hamiltonian and state for an ideal mass-spring for example.\n- IIUC the only test metric is squared deviation from ground truth energy, with no consideration given to whether the actual dynamics are correct. This allows for very incorrect solutions to perform well e.g. a solution where the system sits still in tasks 1-3. I would suggest adding some evaluation metrics to reassure the reader that the dynamics are actually correct.\n- In the constraint for task 4, there are parameters $a_i$ that place upper bounds on the real parts of the eigenvalues of the Jacobian. How are these values chosen? Don't they determine how dissipative the system is, so choosing them \"just right\" already gives your network a significant amount of knowledge? How would we go about choosing these for an unknown dissipative system?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Learning the dynamics of a physical system is often difficult because we can't enforce constraints such as conservation of energy. The authors propose adding a regularization term to the loss function that captures physical constraints. They train models on several benchmark problems and compare to several baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "\n- In Eq. 6, I think there is a `-` sign that is meant to be a `=`.\n- For tasks 1-3 the authors clearly explain how the training data was generated. It would be nice if they could add similar descriptions for the test data in tasks 1-3 and for both datasets in task 4. In particular, they should emphasize that the test data are distinct from the training data.\n- The caption for Table 1 says it shows mean $\\pm$ standard deviation, but the uncertainty is not symmetric. Are these actually the 2.5 and 97.5 percentiles?\n", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes adding a regularization term to the loss function that can capture physical constraints such as energy conservation. It evaluates the method on several simple benchmark problems and compares the results with some natural baselines. 

It was noted in the discussion that adding a regularization penalty is such a common method in machine learning that the results here may not interest the community. Of course, sometimes well-known techniques can lead to important insights if applied in new ways. In this case, the form of the regularizer is perhaps not particularly surprising, but nonetheless there is value in demonstrating its efficacy.

It was also noted that the applications were restricted to simple physical systems, and it was unclear whether the technique would generalize to more complex scenarios. The examples investigated here are indeed fairly simple, but the results are compelling. Furthermore,  no matter where the analysis stopped, there would always be more complex systems out there that were not investigated; it is thus reasonable and natural to limit initial investigations to the simplest settings.

All told, while there are ways the manuscript might be improved in order to further its impact, as it stands there is sufficient value to merit publication.
Decision: Accept (poster)